Implemented Cullum/Willoughby method for removing spurious eigenvalues. by wdj · Pull Request #143 · ORNL-CEES/mfmg

wdj · 2019-03-19T16:49:25Z

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 · 2019-03-19T16:49:27Z

Can one of the admins verify this patch?

aprokop · 2019-03-19T17:21:46Z

ok to test

dalg24 · 2019-03-19T17:50:01Z

OK to test

Rombur · 2019-03-20T19:46:08Z

OK to test

dalg24

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

dalg24 · 2019-03-20T14:43:11Z

+  // an eigenvalue of T2, then it is considered spurious.
+  //
+  // The C/W algorithm has other nuances which may necessitate more
+  // revisions of the simplified algoithm implemented here.


s/algoithm/algorithm/

dalg24 · 2019-03-20T14:44:15Z

+
+  // The following tolerance may need adjustment, however this doesn't
+  // seem critical (cf. C/W vol. 1). Cullum/Willoughby use 1e-12.
+  const double tol2 = 5.e-12;


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

dalg24 · 2019-03-20T15:13:04Z

+    is_repeated[i] = (i>0 && evals[i-1] >= evals[i] - tol2) ||
+                     (i<n-1 && evals[i+1] <= evals[i] + tol2);
+    is_marked[i] = 0==i || evals[i-1] < evals[i] - tol2;
+  } // i


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 · 2019-03-20T15:17:47Z

+    // (only need evals, not evecs, so a different LAPACK call may be adequate)
+
+    std::vector<double> evals2;
+    std::vector<double> evecs2;


Unused variable.

dalg24 · 2019-03-20T16:24:47Z

+    std::cout << std::endl;
+    std::cout << std::endl;
+#endif
+    return std::make_tuple(evals, evecs);


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

dalg24 · 2019-03-20T16:28:26Z

+    evals2.erase(evals2.begin());
+    std::vector<double> sub_diagonal_aux = sub_diagonal;
+    sub_diagonal_aux.erase(sub_diagonal_aux.begin());
+    std::vector<double> evecs_aux(n2 * n2);


Shadowing variables sub_diagonal_aux and evecs_aux .

dalg24 · 2019-03-20T17:35:05Z

+#endif
+
+    // Loop over original eigenvalues to check for spuriousness..
+    // NOTE: assuming here evals and evals2 are in ascending order.


No need it is guaranteed by Lapack

dalg24 · 2019-03-20T17:37:03Z

+  //
+  // The C/W algorithm has other nuances which may necessitate more
+  // revisions of the simplified algoithm implemented here.
+


Please do

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

before the next for-loop.

dalg24 · 2019-03-20T20:16:54Z

+
+      } // n2
+
+    } // n


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 :)

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 · 2019-03-22T14:12:06Z

@wdj I will take care of this PR.

aprokop · 2019-03-22T16:01:45Z

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 · 2019-03-22T19:41:35Z

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 · 2019-03-22T20:52:56Z

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 · 2019-03-22T21:17:28Z

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 · 2019-05-09T02:12:22Z

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 · 2019-05-09T12:36:31Z

So using the coarser tolerance, are the eigenvectors correct?

wdj · 2019-05-09T12:58:04Z

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

wdj · 2019-05-09T13:05:05Z

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

Rombur · 2019-05-09T13:09:14Z

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

wdj · 2019-05-09T13:16:10Z

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

aprokop · 2019-05-09T13:46:15Z

Looking at it right now.

aprokop · 2019-05-09T18:18:08Z

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 · 2019-05-10T18:38:14Z

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 · 2019-05-10T18:49:10Z

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 · 2019-05-10T19:00:59Z

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 · 2019-05-10T19:10:50Z

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.

Implemented Cullum/Willoughby method for removing spurious eigenvalues.

507f6ae

ghost assigned wdj Mar 19, 2019

ghost added the in progress label Mar 19, 2019

Rombur closed this Mar 20, 2019

Rombur reopened this Mar 20, 2019

ghost removed the in progress label Mar 20, 2019

ghost assigned Rombur Mar 20, 2019

ghost added the in progress label Mar 20, 2019

dalg24 reviewed Mar 20, 2019

View reviewed changes

aprokop assigned aprokop and unassigned Rombur Mar 21, 2019

aprokop mentioned this pull request May 16, 2019

Implemented Cullum/Willoughby method for removing spurious eigenvalues. #180

Merged

Conversation

wdj commented Mar 19, 2019

Uh oh!

dalg24-jenkins commented Mar 19, 2019

Uh oh!

aprokop commented Mar 19, 2019

Uh oh!

dalg24 commented Mar 19, 2019

Uh oh!

Rombur commented Mar 20, 2019

Uh oh!

dalg24 left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

aprokop commented Mar 22, 2019

Uh oh!

aprokop commented Mar 22, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

aprokop commented Mar 22, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

aprokop commented Mar 22, 2019

Uh oh!

aprokop commented Mar 22, 2019

Uh oh!

wdj commented May 9, 2019

Uh oh!

Rombur commented May 9, 2019

Uh oh!

wdj commented May 9, 2019

Uh oh!

wdj commented May 9, 2019

Uh oh!

Rombur commented May 9, 2019

Uh oh!

wdj commented May 9, 2019

Uh oh!

aprokop commented May 9, 2019

Uh oh!

aprokop commented May 9, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

wdj commented May 10, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

wdj commented May 10, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Rombur commented May 10, 2019

Uh oh!

wdj commented May 10, 2019

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

aprokop commented Mar 22, 2019 •

edited

Loading

aprokop commented Mar 22, 2019 •

edited

Loading

aprokop commented May 9, 2019 •

edited

Loading

wdj commented May 10, 2019 •

edited

Loading

wdj commented May 10, 2019 •

edited

Loading