ellipse, combine A0 and A1 integral, sort !sanity_check, parallel ptrapz

[wrs28]

Summary of changes

Circle contour extended to ellipse (useful when, say, you are interested in complex eigenvalues close to the real axis).
More useful still would be the ability to define a polygon, which could define a strip along the real or imaginary axes, but that is not done here.

I've combined the quadratures for A0 and A1, which reduces the number of calls to the linear operator by half, leading to significant speedup.

I added a the sorting step for the case when sanity_check=false. This means that the eval/evecs will always be sorted by their distance from σ, regardless of the sanity_check status.

Lastly, I extended the quadrature methods by a parallel version of :ptrapz, called :ptrapz_parallel.

I changed some of the function documentation to reflect this.

---

Timing statistics due to the changes

Setup:

addprocs(7);
using BenchmarkTools, LinearAlgebra, SparseArrays, NonlinearEigenproblems
N=1000;
As = [sparse(1.0I,N,N), spdiagm(0=>1.0*(1:N))];
fns = [identity, one];
nep = SPMF_NEP(As,fns);
@benchmark contour_beyn(nep; neigs=4, σ=-3.6, radius=2, quad_method=:ptrapz)

The results before the changes are

BenchmarkTools.Trial: 
  memory estimate:  895.74 MiB
  allocs estimate:  108568
  --------------
  minimum time:     433.630 ms (16.68% GC)
  median time:      443.142 ms (17.24% GC)
  mean time:        443.792 ms (17.11% GC)
  maximum time:     455.777 ms (17.56% GC)
  --------------
  samples:          12
  evals/sample:     1

The results for the serial code after the changes are

BenchmarkTools.Trial: 
  memory estimate:  563.58 MiB
  allocs estimate:  73531
  --------------
  minimum time:     251.076 ms (19.99% GC)
  median time:      258.093 ms (21.74% GC)
  mean time:        258.392 ms (21.69% GC)
  maximum time:     270.839 ms (23.05% GC)
  --------------
  samples:          20
  evals/sample:     1

The results for the parallel code after changes (there is no pre-change counterpart) (use quad_method=:ptrapz_parallel)

BenchmarkTools.Trial: 
  memory estimate:  4.05 MiB
  allocs estimate:  2108
  --------------
  minimum time:     52.615 ms (0.00% GC)
  median time:      59.058 ms (0.00% GC)
  mean time:        63.584 ms (0.53% GC)
  maximum time:     123.447 ms (0.00% GC)
  --------------
  samples:          79
  evals/sample:     1

[jarlebring]

Sweet!

I think sorting will not work so well if we have an ellipse. Do you think it would work to change the distance measure with appropriate weighting?

Like this:

weightednorm= x-> sqrt(radius[2]*real(x-σ)^2+radius[1]*imag(x-σ)^2)
...
sorted_index=map(weightednorm,λ)

(I think I got the weighting wrong...)

[wrs28]

I guess I have to ask what the significance of sorting is if the contour is not a circle. To me it seems that the idea of returning the neigs eigenpairs that are closest to σ inherently assumes that the search is done over a disk.

If the user chooses not to search over a disk, they should not expect to get the closest eigenpairs.
I see no harm in sorting the eigenvalues by distance from σ, but of course they should not (and cannot) be interpreted as being the closest eigenvalues in an absolute sense.

So while I don't see any harm to defining a weighted norm for the ellipse, I don't really see it as especially meaningful, since "closest" is only a meaningful measure if the contour was a circle to begin with.

What do you think?

[jarlebring]

I view the contour as also the region of interest. That's the same as closest to σ for circles but not ellipses (unless we use weighting). Contour integral methods sometimes capture eigvals outside the contour. Suppose we use

σ=0, radius=[1;3], neigs=1

and that algorithm finds eigenvalues

λ=[1.5im 2]

Then the user will get the eigenvalue λ=1.5im, which is not inside the contour, probably not the eigval she wanted. This will matter if we are interested in (close to) real eigenvalues.

Discarding eigvals outside the contour would also be okay solution for me.

[wrs28]

I see, that's a good point, and one that has come up in my own use.

Would it not be better to address the problem head on by throwing away eigenvalues that are outside the contour, maybe throwing a warning and suggesting higher N or lower tol? This could be done for general contours by checking a winding number, but for an ellipse it can be checked analytically (by comparing the distance to r(\theta) ).

[jarlebring]

Sounds good!

[wrs28]

Great, I'll add an inside-contour check this weekend.

[wrs28]

hmm, I hate to nitpick, but I think the logic surrounding sanity_check is still a little off.

Right now, the user requests neigs eigenpairs, checks the contour for k eigenpairs, and finds p "non-vanishing" singular values which is supposed to be associated with p eigenpairs. So long as neigs≤p<k everything makes sense.

The actual Beyn algorithm as I understand it constructs a matrix B which is p×p, and computes the eigenvalues and eigenvectors from that. But as written, we construct a B which is k×k, and then subsequently check which eigenvectors don't work via errmeasure, and then we subsequently sort and check for being in the interior of the contour.

I think it makes more sense to explicitly construct a p×p matrix B, and then take those results and either pass them to the user or sanity-check them. That way, if the user chooses to forgo the sanity check, they still get eigenpairs that are valid within the algorithm. Whereas now, if there is no sanity check, they get the eigenvectors and eigenvalues of the full k×k B matrix, which even in principle contains fictitious eigenpairs.

[jarlebring]

You're totally right that about Beyn's method and k vs p (seems easy to fix), and indeed it requires reworking the sanity_check code (less obvious).

[wrs28]

I think I fixed it with the most recent commit. The k->p was easy to fix, the sanity check I think is also a matter of taking k->p. I did not look into what errmeasure does, but I assume it just checks the norm of A*v-λv in this case. If that's true, then the sanity check still does what it supposed to, but we don't leave for it the work of filtering out the k-p fictitious eigenpairs.

The commits I just made seem to check out with the example we worked on previously for σ=-9, radius=3.999. This has two eigenpairs very close to the contour (5 and 13), but technically outside of it. Interestingly, in this case, with the default N, the outside eigenpairs are still basically correct.

The most recent commit does not throw away the outside eigenpairs per se, it just moves them to the end of the line and throws a warning. If all neigs are inside the contour, it warns that some were detected outside the contour but not returned, and inaccuracy might be possible. If 1:neigs includes eigenpairs that are outside the contour, a slightly different warning is thrown indicating that the last few eigenpairs are outside the contour, and that the results may be inaccurate.

This gives sufficient warning to the user, guides them on how to improve accuracy, but leaves it to them to decide what to do with the extra eigenpairs. I did this because I have found that knowing what the "wrong" eigenvalue is can help diagnose how to fix it, if it indeed needs fixing.

How does that sound?

[jarlebring]

Looks very good. In particular that you don't discard wrt to boundary but instead put them at the end. This is robust wrt numerical errors that make eigvals jump over boundary.

Is it possible to somehow also avoid throwing a warning if sanity_check=false? If it is actually a real use-case, this warning can get very annoying. It would be less annoying if it could be silenced. Earlier, sanity_check=false implied no warnings.

I will have a look at this PR again tomorrow, but it looks ready to be merged.

[wrs28]

good idea about the sanity_check. If the user goes out of their way to turn it off, they are explicitly choosing not to check for accuracy, and such warning would indeed be unnecessary.

I've just made the change.