Implemented Cullum/Willoughby method for removing spurious eigenvalues.

[wdj]

See attached files for details.

The Lanczos unit tests pass, even with a new, tighter convergence tolerance and use of the percent overshoot option.

spurious_eigenvalues_fix.txt

[dalg24-jenkins]

Can one of the admins verify this patch?

[aprokop]

ok to test

[dalg24]

OK to test

[Rombur]

OK to test

[dalg24]

Make sure to remove the print statements to the standard output and cleanup the test before we merge.

[dalg24]

s/algoithm/algorithm/

[dalg24]

Any reason why you went for 5e-12 instead of 1e-12?

[dalg24]

Readability How about instead of declaring tol2 you introduce something like:

auto const is_strictly_ordered = [] (double eval1, double eval2) { return eval1 + 5e-12 < eval2; };

And then in the body of the for-loop you do:

is_repeated = (i > 0 && !is_strictly_ordered(evals[i-1], evals[i]) ||
              (i < n-1 && !is_strictly_ordered(evals[i], evals[i+1]);
is_marked = (0 == i) || !is_strictly_ordered(evals[i-1], evals[i]);

I am not sure about the name yet but I do see that every time tol2 appears we are checking strict ascending order on eigenvalues and I find it hard to read in the current form (mixing eval1 >= eval2 - eps, eval2 <= eval1 + eps, etc.)

[dalg24]

Unused variable.

[dalg24]

Do not shadow these variables. Clear the vectors before you return or explicitly make a tuple from default constructed vectors.

[dalg24]

Shadowing variables sub_diagonal_aux and evecs_aux .

[dalg24]

Assert it.

[Rombur]

No need it is guaranteed by Lapack

[dalg24]

Please do

ASSERT(std::is_sorted(evals.begin(), evals.end()), "eigenvalues are not in ascending order");

before the next for-loop.

[dalg24]

I feel like it might be more readable with STL algorithms:

auto first = evals.begin();
auto const last = evals.end();
for (int i = 0; i < n; ++i) {
  if (is_repeated[i]) continue;
  first = std::find_if(first, last,
                       [&is_strictly_ordered, &evals, i](double x) {
                         return !is_strictly_ordered(x, evals[i]);
                       });
  if (first != last && !is_strictly_ordered(evals[i], *first)) is_spurious[i] = true;
}

@aprokop @Rombur Any opinion on this?

Disclaimer I have not actually compiled and tested the pseudo code I wrote :)

[wdj]

One thing about the issue of handling the (a-tol, a+tol) interval. In the CW algorithm some attention is given to making the intervals exactly consistent regarding the multiplicity test and the spuriousness test. it is possible but I think doubtful that the inclusion or exclusion of the endpoints of the interval would make much difference. So I think the simplification of the code is probably ok --

[aprokop]

@wdj I will take care of this PR.

[aprokop]

Running for 60 eigenvalues (non-deflated mode, multiplicity 1) produces correct eigenvalues but wrong eigenvectors. When looking for 15 or 30 eigenvectors, everything seems to be fine.
An example of the wrong eigenvector:

-1.008e+00 3.550e+00 -1.950e-04 -6.619e-06 2.551e-05 ...

which is completely off for a diagonal matrix.

failed_60.txt

[aprokop]

Unfortunately, it seems that the new approach also has issues with multiplicity > 1 cases. For example, searching for 20 eigenvectors in a diagonal matrix with multiplicity 2 eigenvalues results in a weird situation where during the second pass of deflation all tridiagonal eigenvalues are declared to be spurious. And, in fact, the eigenvalues of both Lanczos T matrix and the reduced T2 matrix become negative, which is absurd.

failed_10_multiplicity2.txt

[aprokop]

Reading more about Cullum, it seems that instead of computing the T2 eigenvalues, it is more efficient to use Sturm's theorem to just find the number of T2's eigenvalues in intervals surrounding T1's eigenvalues. This is a significantly cheaper calculation.

[aprokop]

Further checking reveals that for multiplicity > 1 the wrong eigenvectors are coming from the first Lanczos solve. When searching for 10 eigenvalues with a single Lanczos solve for a multiplicity 2 diagonal matrix, the found eigenvalues are correct, but eigenvectors are not.

[wdj]

Running for 60 eigenvalues (non-deflated mode, multiplicity 1) produces correct eigenvalues but wrong eigenvectors. When looking for 15 or 30 eigenvectors, everything seems to be fine.
An example of the wrong eigenvector:

-1.008e+00 3.550e+00 -1.950e-04 -6.619e-06 2.551e-05 ...

which is completely off for a diagonal matrix.

failed_60.txt

There are two problems here. First there is a bug such that if max iterations is hit, the code does not recover gracefully but gives garbage. Adding the following line immediately after the lanczos iteration loop fixes this: "it = it < maxit ? it : maxit;". Second, the convergence tolerance is much smaller than the solver can handle. I was able to get convergence for these combinations of tolerances and required iterations: 1e-5/703, 1e-6/739, 1e-7/1791, 1e-8/2651. So 1e-11 is unfortunately far smaller than what can be handled practically. Not sure immediately whether this is due to the eigenproblem itself, the algorithm, or the implementation. Hopefully a smaller tolerance will meet the needs for the mfmg algorithm - ?

[Rombur]

So using the coarser tolerance, are the eigenvectors correct?

[wdj]

Yes, the test passes (eigenvalues, eigenvectors) --

[wdj]

With the change, if max its is reached without convergence, the best approximate eigenvectors are returned, rather than garbage --

[Rombur]

OK very nice! Thanks Wayne. @aprokop can you take a look at this PR

[wdj]

I haven't checked the multiplicity problem yet, this may or may not fix it --

[aprokop]

Looking at it right now.

[aprokop]

Multiplicity 1 works (I used 1e-4 tolerance), tested for up to 100 eigenvalues/eigenvectors. Multiplicity > 1 with deflation breaks down for larger number of eigenvalues (for example, multiplicity 2 works with 40 distinct eigenvalues, but breaks down with 50). The spurious eigenvalues are there, and present an interesting case of having some eigenvalues be correct but with a higher multiplicity than in the original matrix.

This is certainly an improvement, and the non-deflated version of the method can be considered ready.

[wdj]

One problem with the deflated method is that for some reason it is hard to drive the eigenpair residual norm low enough when in a deflation cycle. However I'm having some success solving this problem with the code below, which picks, from a set of Lanczos tridiagonal eigenpairs with (artificial) multiplicity, the eigenvalue from the cluster of repeats with the lowest residual error estimate.

  // Identify repeated eigenvalues.
  // A "marked" eigenvalue here is nonrepeated or best of a set of repeated.

  std::vector<bool> is_repeated(n);
  std::vector<bool> is_marked(n);
  std::vector<int> elt_in_cluster(n);

  int i_marked_in_cluster = -1;
  double sn_i_marked_in_cluster = 0;
  for (int i = 0; i < n; ++i) {
    const bool is_lower_repeated = i>0 && evals[i-1] >= evals[i] - tol2;
    const bool is_upper_repeated = i<n-1 && evals[i+1] <= evals[i] + tol2;
    is_repeated[i] = is_lower_repeated || is_upper_repeated;
    const bool is_cluster_start = ! is_lower_repeated;
    const bool is_cluster_end = ! is_upper_repeated;

    //is_marked[i] = is_cluster_start; // old way

    const double sn = std::abs(evecs_aux[n - 1 + n * i]);

    if (is_cluster_start) {
      elt_in_cluster[i] = 0;
      is_marked[i] = true;
      i_marked_in_cluster = i;
      sn_i_marked_in_cluster = sn;
    } else {
      elt_in_cluster[i] = elt_in_cluster[i-1] + 1;
      if (is_repeated[i] && sn < sn_i_marked_in_cluster) {
        is_marked[i_marked_in_cluster] = false;
        is_marked[i] = true;
        i_marked_in_cluster = i;
        sn_i_marked_in_cluster = sn;
      }
    }

    if (is_cluster_end) {
      i_marked_in_cluster = -1;
    }
  } // i

[wdj]

The other problem is that the deflation is not "perfect" -- a zero eigenvalue is being introduced from the projection -- there is some "noise" from iterating on the projected operator that introduces components in the null space.

I have some ideas for solving this. I need to ask, can it be assumed that the operator is positive semidefinite (or, even better, positive definite)?

[Rombur]

I need to ask, can it be assumed that the operator is positive semidefinite (or, even better, positive definite)?

It is positive semidefinite but zero is a valid eigenvalue.

[wdj]

Got it. If in Lanczos the original operator is shifted up by 1 (times identity), then in that case the spectrum of the original operator will be clearly distinguished from a zero eigenvalue, then after the solve and removal of spurious zero one can shift down by 1 to correct. I'm going to experiment with something like that.