bugGNU Octave - Bugs: bug #53700, eigs test failure related to...

 
 

bug #53700: eigs test failure related to ARPACK generating real NaN rather than complex NaN+1i*NaN

Submitter:  Dan Sebald <sebald>
Submitted:  Thu 19 Apr 2018 04:50:06 PM UTC
   
 
Category:  Test Suite Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 26 Apr 2018 12:39:28 PM UTC, comment #33: 

Marco: I pushed the changeset from comment #32 to stable.

Can this report be closed now, or is more work needed for the 4.4 release?

John W. Eaton <jwe>
Group administrator
Thu 26 Apr 2018 06:49:17 AM UTC, comment #32: 

I would stay with if (isreal (d))... since we already did it. I attach a changeset.

(file #44034)

Marco Caliari <caliari>
Group Member
Wed 25 Apr 2018 08:40:55 PM UTC, comment #31: 

Is more work on this problem needed for 4.4?

John W. Eaton <jwe>
Group administrator
Tue 24 Apr 2018 05:39:39 PM UTC, comment #30: 

The warning can be turned off with:


warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");


So the suggestion of


+%! d = eigs (A, floor((length(A)-1)/2), "sm", opts);


won't work.  How about floor(length(A)/4) rather than floor((length(A)-1)/2)?  That would put the NEV somewhere in the middle of the range where the Krylov subspace isn't adequate.

It's up to Marco to make the choice.  Check for real/imaginary or pick an NEV somewhere that generates lots of 1.0+-iZi and NaN+iNaN.  Probably doesn't matter too much.  We're seeing so much variation in the convergence behavior for all the different libraries and compilers out there.

Dan Sebald <sebald>
Tue 24 Apr 2018 05:27:52 PM UTC, comment #29: 

After the patch on the stable branch, I see the following behavior"


octave:1> A = magic (100) / 10 + eye (100);
octave:2> opts.v0 = (1:100)';
octave:3> opts.maxit = 10;
octave:4> for i=9:20; d = eigs (A, i, "sm", opts); printf ("%d: %s", i, disp(d(end))); end
warning: eigs: Only 2 of the 9 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
9:  NaN + NaNi
warning: eigs: Only 2 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
10:  NaN + NaNi
warning: eigs: Only 4 of the 11 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
11:  NaN + NaNi
warning: eigs: Only 4 of the 12 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
12:  NaN + NaNi
warning: eigs: Only 4 of the 13 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
13:  NaN + NaNi
warning: eigs: Only 6 of the 14 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
14:  NaN + NaNi
warning: eigs: Only 7 of the 15 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
15:  NaN + NaNi
warning: eigs: Only 7 of the 16 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
16:  NaN + NaNi
warning: eigs: Only 9 of the 17 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
17:  NaN + NaNi
warning: eigs: Only 13 of the 18 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
18:  NaN + NaNi
warning: eigs: Only 11 of the 19 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
19:  NaN + NaNi
warning: eigs: Only 15 of the 20 requested eigenvalues converged
warning: called from
    eigs at line 285 column 18
20:  NaN + NaNi


I begin to see non-NaN results at i=49:


49:  1.0000e+00 - 9.9142e-13i
50:  1.0000e+00 + 4.9814e-13i




John W. Eaton <jwe>
Group administrator
Tue 24 Apr 2018 03:28:20 PM UTC, comment #28: 

Do the two patches pushed today for bug #53719 also resolve this report?

Rik <rik5>
Group administrator
Mon 23 Apr 2018 04:31:04 PM UTC, comment #27: 

OK, got it.

Note that "lm" ends up being a different case; it doesn't necessarily run the same ARPACK function or the same implementation of that function.  So let's not pursue that idea.

Also, an observation here is that on Mint 18.2 I find that ARPACK issues error -14 when selecting NEV to be 8 or less.  On Ubuntu 18.04 (the system that produces one convergent real 1.0 eigenvalue) it fails for NEV equal to 5 or less:


octave:16> d = eigs (A, 5, "sm", opts);
error: __eigs__: eigs: error -14 in dneupd
error: called from
    eigs at line 265 column 18
octave:16> d = eigs (A, 6, "sm", opts);


So, lots of library/compiler variation going  on.  (I don't even see a -14 error listed in the documentation I cited previously.)

Nonetheless, I suspect that the real/complex issue is more likely to occur when NEV is just above that transition where it doesn't produce an error.  A more assured way of always getting a complex result would be to push the NEV higher:


octave:22> for i=9:50; d = eigs (A, i, "sm", opts); printf("%s%s", printf("%d: ",i), disp(d(end))); end
9:  NaN + NaNi
10:  NaN
11:  NaN + NaNi
12:  NaN + NaNi
13:  NaN + NaNi
14:  NaN + NaNi
15:  NaN + NaNi
16:  NaN + NaNi
17:  NaN + NaNi
18:  NaN + NaNi
19:  NaN + NaNi
20:  NaN + NaNi
21:  NaN + NaNi
22:  NaN + NaNi
23:  NaN + NaNi
24:  NaN + NaNi
25:  NaN + NaNi
26:  NaN + NaNi
27:  NaN + NaNi
28:  NaN + NaNi
29:  NaN + NaNi
30:  NaN + NaNi
31:  NaN + NaNi
32:  NaN + NaNi
33:  NaN + NaNi
34:  NaN + NaNi
35:  NaN + NaNi
36:  NaN + NaNi
37:  NaN + NaNi
38:  NaN + NaNi
39:  NaN + NaNi
40:  NaN + NaNi
41:  NaN + NaNi
42:  NaN + NaNi
43:  NaN + NaNi
44:  NaN + NaNi
45:  NaN + NaNi
46:  NaN + NaNi
47:  NaN + NaNi
48:  NaN + NaNi
49:  NaN + NaNi
50:  1.0000e+00 + 6.6006e-13i


That is, the NEV=49 case has so many converged eigenvalues of 1.0 that there is bound to be complex conjugate pairs, and also plenty of NaN + NaNi as well.  So the following would also be an easy way of fixing this test:


--- a/scripts/sparse/eigs.m     Thu Apr 19 01:58:44 2018 -0400
+++ b/scripts/sparse/eigs.m     Mon Apr 23 11:28:32 2018 -0500
@@ -1409,8 +1409,8 @@
 %! opts.v0 = (1:100)';
 %! opts.maxit = 10;
 %! warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
-%! d = eigs (A, 10, "sm", opts);
-%! assert (d(10), NaN+1i*NaN);
+%! d = eigs (A, floor((length(A)-1)/2), "sm", opts);
+%! assert (d(end), NaN+1i*NaN);
 %!testif HAVE_ARPACK
 %! A = toeplitz ([0, 1, zeros(1, 8)], [0, -1, zeros(1, 8)]);
 %! A = kron (A, eye (10)) + kron (eye (10), A);


Dan Sebald <sebald>
Mon 23 Apr 2018 02:13:49 PM UTC, comment #26: 

@Dan: for the "lm" version, I have real NaNs both with 9 and 10.
NEV to 50 is p to 100.

Marco Caliari <caliari>
Group Member
Mon 23 Apr 2018 08:10:10 AM UTC, comment #25: 

Oh, so the reason for the less robust system for NEV=10 is that the Krylov subspace dimensionality is not even near the cardinality of the largest eigenvalue group, which is 97 in this case.  I.e., 97 eigenvalues of 1.0.  That's what the algorithm is struggling with.  On the other hand, had we chosen a system for which there were more variation in the eigenvalues, e.g., 1.0 1.0 3.0 3.0 3.5 2-2j 2+2j, etc., and used a Krylov subspace of 10 then we'd get good convergence.

Going the other direction, increasing NEV to 50 (or is it p=50?) makes the Krylov subspace nearer to the size of the largest eigenvalue group in this case.

Thanks for the explanation.

Would there be any other type of more assured test?  Say, if instead of looking for "sm" eigenvalues, we looked for "lm".  The largest eigenvalues have a group cardinality of 1 or 2 and therefore those few really large eigenvalues should converge nicely on all systems, compilers and linear algebra libraries with a choice of NEV=10.  I'm seeing


octave:32> d = eigs (A, 10, "lm", opts);
octave:33> d
d =

   50006.00000
    2887.60700
   -2885.60700
       1.00000
           NaN
           NaN
           NaN
           NaN
           NaN
           NaN


But that too ends up being not complex, so I don't know if that gets us anywhere different.

How about choosing an odd number for NEV?  Will that force there to be a conjugate pair somehow?


octave:34> d = eigs (A, 9, "lm", opts);
octave:35> d
d =

   50006.00000 +     0.00000i
    2887.60700 +     0.00000i
   -2885.60700 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
           NaN +         NaNi
           NaN +         NaNi
           NaN +         NaNi
           NaN +         NaNi


What do others see for such tests?  Is it just as arbitrary?

Dan Sebald <sebald>
Mon 23 Apr 2018 06:59:06 AM UTC, comment #24: 

@Dan: I try to make a cumulative answer to your deep analysis. The main "problem" in ARPACK is a proper choice of nev (number of desired eigenvalues, k in eigs) and ncv (dimension of the Krylov subspace of approximate eigenpairs (p in eigs). As written at the end of section 2.3.3 of Arpack User's Guide, "it is best to avoid setting nev in such a way that will split clusters of eigenvalues". The choice 10 is then to avoid (but it was made in order to produce a failure!). Later "as a rule of thumb, ncv >= 2*nev is reasonable". ncv=2*nev is the choice made in eigs. Therefore, requiring 10 eigenvalues correspond to using 20 as Krylov subspace. Requiring 50 eugenvalues correspond to using 100 as Krylov subspace, which is the dimension of the original matrix: therefore, you are pretty sure that you find all the eigenvalues. You can force a new ncv by setting opts.p="something larger than 20" and keeping k=10. Solution will improve.
Z = B * X (that is Z = X, in our case) is done in dnaup2, which is called be dnaupd.
Matlab can find only few eigenvalues, too.

Therefore, although it is not possible to exclude a hidden bug in the chain eigs.m, _eigs_.cc, eigs-base.cc, arpack, "Why when we request less eigenvalues do we have a less robust situation?" can be explained as above and it is not a clear symptom of a bug.

Marco Caliari <caliari>
Group Member
Mon 23 Apr 2018 03:46:28 AM UTC, comment #23: 

I got 4 with atlas and 2 with ref blas and openblas.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Mon 23 Apr 2018 03:27:32 AM UTC, comment #22: 

@Dmitri: Could you look at the results for the failing case for all those versions and see if the number of converged eigenvalues is the same?  The code is:


A = magic (100) / 10 + eye (100);
opts.v0 = (1:100)';
opts.maxit = 10;
d = eigs (A, 10, "sm", opts);


I've experimented with things like

opts.maxit = 100;
opts.tol = 1e-13;

to see at what point I begin to see all the eigenvalues converging.  I'm surprised that I have to set the tolerance to 1e-13 to 1e-12 to get all ten eigenvalues converging.  That's four orders of magnitude above the floating point resolution.

Again, picking the number of eigenvalues to be large seems to make convergence better, so it has me wondering if there is some error in the loop.  I've experimented a bit, but can't uncover anything.  The dnaupd function:

http://www.caam.rice.edu/software/ARPACK/UG/node137.html


      F77_FUNC (dnaupd, DNAUPD)
        (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
         F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
         k, tol, presid, p, v, n, iparam,
         ipntr, workd, workl, lwork, tmp_info
         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));


wants the programmer to supply the following when ido=1:


c          IDO =  1: compute  Y = OP * Z  and Z = B * X where
c                    IPNTR(1) is the pointer into WORKD for X,
c                    IPNTR(2) is the pointer into WORKD for Y,
c                    IPNTR(3) is the pointer into WORKD for Z.


Currently, Y = OP * Z is done using a L U decomposition approach that appears to be a sparse routine with P[] and Q[] indexing.  The P Q pairs are pretty much 1,1 2,2 3,3 etc. but along the way the indices get out of sync by 1, like 98 takes the place of 6.  I really don't know the sparse code well enough to say why it is choosing the indices in that way.

I also noticed that currently there is no explicit code that sets Z = B * X where in this case B is the identity matrix.  So, I added a little hunk of code that will do that, and it didn't seem to affect the result any.  I suppose the DNAUPD is setting vector Z = X by default.

Otherwise, all the code associated with EigsRealNonSymmetricMatrixShift appears to be reasonable in terms of how the library function should be used.  If anything, my concern would be that the L-U sparse decomposition carries out Y = OP * Z precisely.

I think Marco's suggestion is the best at this point.  I'm sure I won't be the only person who's system produces a single convergent eigenvalue that doesn't come in a conjugate pair.

Dan Sebald <sebald>
Fri 20 Apr 2018 06:54:37 PM UTC, comment #21: 

it works for me wit hany blas:


LD_PRELOAD=/usr/lib64/atlas/libtatlas.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.33795 seconds.

LD_PRELOAD=/usr/lib64/atlas/libsatlas.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.33742 seconds.

LD_PRELOAD=/usr/lib64/libopenblaso.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.68253 seconds.
octave:2> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.50371 seconds.

LD_PRELOAD=/usr/lib64/libopenblasp.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.68722 seconds.


LD_PRELOAD=/usr/lib64/libopenblas.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.42171 seconds.


LD_PRELOAD=/usr/lib64/libblas.so ./run-octave  -f -q
octave:1> tic; test eigs; toc
PASSES 179 out of 179 tests
Elapsed time is 1.40842 seconds.


Notice the times are different wrt blas version -- atlas seems faster than ref/openblas. So it is possible that actual solutions are different. (I think the size of the problem is too small to judge the efficiency of particular blas implementation.)

Dmitri.
--

 

Dmitri A. Sergatskov <dasergatskov>
Fri 20 Apr 2018 06:40:37 PM UTC, comment #20: 

I see there is a difference in algorithms when running

d = eigs (A, 10, "sm", opts);

versus

d = eigs (A, 100, "sm", opts);

The former doesn't use __eigs__.  Rather it uses the, presumably, non-sparse eig() routine.  So that is one difference, but more investigation shows the situation hasn't changed.

d = eigs (A, 100, "sm", opts);
d = eigs (A, 99, "sm", opts);
d = eigs (A, 98, "sm", opts);

use the eig() routine.  For any k <= 97, the routine that is used is __eigs__().  In eigs.m, here is the hunk of code that makes that decision:


        if (rows (a) - k < 3)
          call_eig = true;
        endif


I suppose that is a crude guess as to whether something might be sparse.  So, let's try


octave:42> d = eigs (A, 97, "sm", opts);
octave:44> size(d(d==d))
ans =

   97    1

octave:45> d(97)
ans =  1.0000e+00 + 2.1304e-13i


The algorithm behaves well when k is large (those last three eigenvalues are the non 1.0 values).

Where does the algorithm seem to break?  At 50, which happens to be half the number of columns or eigenvalues, and that probably isn't a coincidence either.


octave:63> d = eigs (A, 50, "sm", opts);
octave:64> d = eigs (A, 49, "sm", opts);
warning: eigs: Only 43 of the 49 requested eigenvalues converged
warning: called from
    eigs at line 268 column 18


and it so happens that the number of iterations through that do{} loop jumps from 100 to 243 before the ido=99 comes back when transitioning from k=50 to k=49.

I've got a feeling there is some type of bug in Octave's do{} loop.  In the following do loop there are several comments about exactly where some important information lies in Eig ARPACK result:


  do
    {
std::cerr << "counterval " << counterval++ << "\n";
      F77_INT tmp_info = octave::to_f77_int (info);

      // On exit, ip(4) <= k + 1 is the number of converged eigenvalues.
      // See dnaupd2.
      F77_FUNC (dnaupd, DNAUPD)
        (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
         F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2),
         k, tol, presid, p, v, n, iparam,
         ipntr, workd, workl, lwork, tmp_info
         F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));
      // k is not changed

std::cerr << "ido = " << ido << "\n";
      info = tmp_info;

      if (disp > 0 && ! octave::math::isnan (workl[iptr (5)-1]))
        {
          if (iter++)
            {
              os << "Iteration " << iter - 1 <<
                 ": a few Ritz values of the " << p << "-by-" <<
                 p << " matrix\n";
              if (ido == 99) // convergence
                {
                  os << "    " << workl[iptr(5)+k] << "\n";
                  for (F77_INT i = 0; i < k; i++)
                    os << "    " << workl[iptr(5)+i-1] << "\n";
                }
              else
                {
                  // the wanted Ritz estimates are at the end
                  for (F77_INT i = p - k - 1; i < p; i++)
                    os << "    " << workl[iptr(5)+i-1] << "\n";
                }
            }

          // This is a kludge, as ARPACK doesn't give its
          // iteration pointer.  But as workl[iptr(5)-1] is
          // an output value updated at each iteration, setting
          // a value in this array to NaN and testing for it
          // is a way of obtaining the iteration counter.
          if (ido != 99)
            workl[iptr(5)-1] = octave::numeric_limits<double>::NaN ();
        }


Could there be something strange going on with the location of some of these results when k in the ARPACK is below N_cols/2 as opposed to greater than N_cols/2?  And that information is being used to incorrectly "prime" the next pass of the loop?

I don't have time to look into this, but it sure seems to me that the Eigs ARPACK routines shouldn't be failing for this example, i.e., there should be all convergent eigenvalues and no NaNs.

Dan Sebald <sebald>
Fri 20 Apr 2018 05:23:24 PM UTC, comment #19: 

At first I thought there was something wrong about the system setup, i.e., this large matrix with numbers much larger than eye() caused an instability somewhere.  However, experimenting a bit, I don't quite understand what is happening with the number of eigenvalues.  Going back to the original test:


octave:2> A = magic (100) / 10 + eye (100);
octave:3> opts.v0 = (1:100)';
octave:4> opts.maxit = 10;
octave:5> d = eigs (A, 10, "sm", opts);
warning: eigs: Only 1 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
octave:6> d(d==d)
ans =  1.00000


However, change the requested number of eigenvalues from 10 to 100, i.e., all the eigenvalues and surprisingly there is no convergence problem:


octave:7> d = eigs (A, 100, "sm", opts);
octave:8> d(d==d)
ans =

       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 -     0.00000i
       1.00000 +     0.00000i
       1.00000 +     0.00000i
   -2885.60700 +     0.00000i
    2887.60700 +     0.00000i
   50006.00000 +     0.00000i


Doing sum(A), every value is 50006, and it probably isn't a coincidence that the last eigenvalue in the list is 50006.  Seems reasonable.

So I'm confused.  Why when we request less eigenvalues do we have a less robust situation?  Typically when one puts a restriction on an algorithm it leads to improved robustness...say if we know some eigenvalues are problematic we might choose sigma and the number of eigenvalues in a way that squeezes those problematic eigenvalues out of the picture.

Dan Sebald <sebald>
Fri 20 Apr 2018 05:08:04 PM UTC, comment #18: 

I think Marco's fix is the right one. The purpose of the test is to verify that the function was asked for 10 eigenvalues, a vector of length 10 is returned, and the last one is either NaN or (NaN,NaN). It doesn't matter whether 1 or 5 valid eigenvalues were returned, as long as it's less than 10, and it doesn't matter what their values are.

Mike Miller <mtmiller>
Group Member
Fri 20 Apr 2018 05:04:20 PM UTC, comment #17: 

@Marco: I don't see any problems with checking that the eigenvalues are within a certain tolerance.  When it comes to applications, users often have to do such things in order to interpret results.

Dan Sebald <sebald>
Fri 20 Apr 2018 04:39:17 PM UTC, comment #16: 

@Mike: I didn't see your post about Ubuntu 18.04 until now.  Your result is different than what you reported from a different system, so technically it isn't reproduced.  On the older system you have


>> imag(d)
ans =

   0.0000e+00
  -8.6837e-13
   8.6837e-13


and on Ubuntu 18.04 you have


>> d(idx)
ans =

   1.00000 + 0.00000i
   1.00000 - 0.00000i
   1.00000 + 0.00000i
   1.00000 - 0.00000i


So, I suspect that the algorithm is finding eigenvalues by iteratively solving some equation that could result in a negative residual somewhere that results in some imaginary pair, and if said residual is positive the result is a single real value.  Just guessing, but these numerical algorithms are prone to this sort of thing with a large space and lots of iterations.

In other words, my setup could be slightly different in terms of compilers and libraries I installed that gives me a real 1.0 eigenvalue in the first location rather than a complex 1.0 +- 0.0 pair.

When I think of this test, I don't know what it is supposed to reveal.  It's sort of testing that realm of the algorithm where there is numerical breakdown, i.e., it's analyzing the limitation of the algorithm and how it falls apart for a larger, over-parameterized system.  Just exactly how the algorithm starts to fail is highly specific to the numerical methods--very chaotic in the sense that one tiny perturbation will produce some completely different way it fails.  And to illustrate, let's change the allowed iterations on my system


octave:25> A = magic (100) / 10 + eye (100);
octave:26> opts.v0 = (1:100)';
octave:27> opts.maxit = 100;
octave:28> d = eigs (A, 10, "sm", opts);
warning: eigs: Only 5 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 268 column 18
octave:29> d(d==d)
ans =

   1.00000 + 0.00000i
   1.00000 - 0.00000i
   1.00000 + 0.00000i
   1.00000 - 0.00000i
   1.00000 + 0.00000i


The number of iterations is almost like a random stopping criteria.

Wouldn't it be more useful to have a test that is acting in a more robust region of the algorithm?

Dan Sebald <sebald>
Fri 20 Apr 2018 08:45:23 AM UTC, comment #15: 

@Dan: on my system, d(1) is


9.999999999991880e-01 + 9.837437573694459e-13i


so, it is not possible to assert (d(1), 1). If we introduce a tolerance, then other bug reports will pop up sayng that the tolerance is not met. I would stay with if (isreal (d)) ...

Marco Caliari <caliari>
Group Member
Fri 20 Apr 2018 08:29:01 AM UTC, comment #14: 

@Marco: Good suggestion.  I reached a similar conclusion because I think there isn't anything wrong with setting up the routines.  It is simply that the library function is behaving slightly different in terms of convergence:


  F77_INT k0 = k;  // original number of eigenvalues required
  F77_FUNC (dneupd, DNEUPD)
    (rvec77, F77_CONST_CHAR_ARG2 ("A", 1), sel, dr, di, z, n, sigmar,
     sigmai, workev,  F77_CONST_CHAR_ARG2 (&bmat, 1), n,
     F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
     ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)
     F77_CHAR_ARG_LEN(2));
  // On exit, if (and only if) rvec == true, k may have been increased by one
  // and be equal to ip(4), see dngets.

std::cerr << "ip(4) = " <<  ip(4) << "\n";
  if (! rvec && ip(4) > k)
    k = ip(4);


ip(4) is 1, not 3 in my case.  It could be that convergence of the first one is solid, but the other two are tenuous and dependent on the compiler.  Looking at what Mike listed:


>> imag(d)
ans =

   0.0000e+00
  -8.6837e-13
   8.6837e-13
          NaN
          NaN
          NaN
          NaN
          NaN
          NaN
          NaN


if only that first eigenvalue is examined, then the followup code  that checks for the imaginary component is not going to set "have_cplx_eig = true;".  It's only the second and third cases where have_cplx_eig is set true in Mike's case.

In addition to the isreal() check, how about testing that d(1) == 1.0 or (1.0,0.0) in addition to the check on d(10)?  It seems that's the reliable eigenvalue.

Dan Sebald <sebald>
Fri 20 Apr 2018 06:50:33 AM UTC, comment #13: 

There is already a test in eigs of type


%! if (isreal (d))
%!   assert (d(4), NaN);
%! else
%!   assert (d(4), NaN +1i*NaN);
%! endif


Look for bug #49523, starting from comment #41. It was a blas-dependent result. So, the easy fix is to introduce if-else-endif in this test, too.

@Dmitri: can you check with different blas versions (ref, openblas, atlas)?

Marco Caliari <caliari>
Group Member
Fri 20 Apr 2018 06:22:13 AM UTC, comment #12: 

I can't reproduce on Ubuntu 18.04. Building with OpenBLAS, all tests pass, and I get


>> A = magic (100) / 10 + eye (100);
>> opts.v0 = (1:100)';
>> opts.maxit = 10;
>> d = eigs (A, 10, "sm", opts);
warning: eigs: Only 4 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
>> idx = find (d == d);
>> d(idx)
ans =

   1.00000 + 0.00000i
   1.00000 - 0.00000i
   1.00000 + 0.00000i
   1.00000 - 0.00000i

>> imag (d(idx))
ans =

   2.2485e-13
  -2.2485e-13
   1.1523e-12
  -1.1523e-12


Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 10:35:24 PM UTC, comment #11: 

Looking at the code for ARPACK, I think the 'rvec' parameters are input-only, so the changes you have made should be correct. But if they don't make any difference in your case, then needs more investigation.

I don't have a Ubuntu 18.04 to test on at the moment, but I'm sure I will at some point, so I can at least try to reproduce this.

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 10:26:16 PM UTC, comment #10: 

Sorry, I didn't see the attachment until now. Yes, you have modified all of the functions as I was suggesting. I'm not sure if rvec is an in/out parameter or if it's input-only. Either way, this was only a guess, since I have seen other problems with arrays of C++ bools being passed in to Fortran routines.

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 10:04:53 PM UTC, comment #9: 

Dan - no, I meant it is following the right path. It's a very large if-else so it was hard to associate your debug prints with line numbers. But it's supposed to be the condition where A is not sparse, not complex, A is not a function handle, sigma is zero, have_sigma is true, symmetric is false. I think it's the case you labeled "     NOT SYMMETRIC".

I think you changed the wrong routine in eigs-base.cc, it should be 'EigsRealNonSymmetricMatrixShift'.

Since it works fine in Debian and in Fedora with the same version of ARPACK, this could be due to some difference like default compiler options used in Ubuntu (which have definitely made a difference in other situations).

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 09:27:54 PM UTC, comment #8: 

FWIW -- on Fedora 28 arpck 3.5.0 -- no problem

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 19 Apr 2018 08:49:07 PM UTC, comment #7: 

You meant "not be doing the right thing"?  I don't see the

std::cerr << "PROMOTED REAL TO COMPLEX\n";

printed out.

Anyway, I see what you mentioned, in the eigs-base.cc routine.  Attached is the changes I've currently made.  Are those F77_LOGICALs what you had in mind?  I.e., basically do the following in all locations:


@@ -627,6 +627,7 @@
   M b(_b);
   F77_INT n = octave::to_f77_int (m.cols ());
   F77_INT mode = 1;
+  F77_LOGICAL rvec77 = rvec;
   bool have_b = ! b.isempty ();
   bool note3 = false;
   char bmat = 'I';
@@ -839,7 +840,7 @@
   double *d = eig_val.fortran_vec ();

   F77_FUNC (dseupd, DSEUPD)
-    (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
+    (rvec77, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
      F77_CONST_CHAR_ARG2 (&bmat, 1), n,
      F77_CONST_CHAR_ARG2 ((typ.c_str ()), 2), k, tol, presid, p, v, n, iparam,
      ipntr, workd, workl, lwork, info2 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1)


Those mods don't change the outcome.

(file #43976)

Dan Sebald <sebald>
Thu 19 Apr 2018 08:02:43 PM UTC, comment #6: 

This case should have fallen into the "promoted real to complex" case, it seems to be doing the right thing for you.

I just noticed that the 'rvec' parameter is passed in as a C++ 'bool' type. But it is not properly converted to a 'F77_LOGICAL' type, it is just thrown on the call stack for each of the ARPACK routines DSEUPD, DNEUPD, ZNEUPD. Fixing that may resolve this for you.

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 07:38:30 PM UTC, comment #5: 

OK, yeah 3.5.0+real-2 here as well (64-bit build).

Here's settings for what could be relevant:


  BLAS libraries:                -lopenblas
[snip]
  Build static libraries:               no
  Build shared libraries:               yes
  Dynamic Linking:                      yes (dlopen)
  64-bit array dims and indexing:       yes
  64-bit BLAS array dims and indexing:  no
  OpenMP SMP multithreading:            yes
  Truncate intermediate FP results:     yes


I don't have a debugger set up at the moment and it is inconvenient to jump back and forth between systems, so let me just print out some stuff as a start...

It appears to me that the tests jump back and forth between complex versions of the library call and real versions of the library call.  The one that fails is real.  I.e., in the following:


  octave_idx_type nconv;
  if (a_is_complex || b_is_complex)
    {
std::cerr << "I AM COMPLEX\n";
      ComplexMatrix eig_vec;
      ComplexColumnVector eig_val;
[snip]
      if (nargout < 2)
        retval(0) = eig_val;
      else
        retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
    }
  else if (sigmai != 0.)
    {
std::cerr << "PROMOTED REAL TO COMPLEX\n";
      // Promote real problem to a complex one.
      ComplexMatrix eig_vec;
      ComplexColumnVector eig_val;
[snip]
     if (nargout < 2)
        retval(0) = eig_val;
      else
        retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
    }
  else
    {
std::cerr << "I AM REAL\n";
      if (symmetric)
        {
std::cerr << "    SYMMETRIC\n";
          Matrix eig_vec;
          ColumnVector eig_val;
[snip]
         if (nargout < 2)
            retval(0) = eig_val;
          else
            retval = ovl (eig_vec, DiagMatrix (eig_val), double (info));
        }
      else
        {
std::cerr << "    NOT SYMMETRIC\n";
          ComplexMatrix eig_vec;
          ComplexColumnVector eig_val;
[snip]
          if (nargout < 2)
            retval(0) = eig_val;
          else
            retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
        }


Notice in the above that the I AM REAL case is the only one which returns a real vector, otherwise the rest return complex.  Here is the output of that failed case:


octave:2> A = magic (100) / 10 + eye (100);
octave:3> opts.v0 = (1:100)';
octave:4> opts.maxit = 10;
octave:5> d = eigs (A, 10, "sm", opts);
ans = eigs 3
CALLED __eigs__
I AM REAL
     NOT SYMMETRIC
warning: eigs: Only 1 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 268 column 18


The NOT SYMMETRIC case has


          ComplexMatrix eig_vec;
          ComplexColumnVector eig_val;


so one would think that is the octave_value that is returned because I don't see anywhere in which there is conditional code that would check the size of the imaginary portion, say, and discard it if the value is negligible.

How Octave is returning a non-complex value has me wondering.  But the interesting thing is that even though that last case (i.e., REAL and NOT SYMMETRIC, lines 517 through 555 of _eigs_.cc) uses a ComplexColumnVector for the return value, the routines of the ARPACK are all of the EigsRealXYZ() variety, not complex.  (In the I AM COMPLEX and PROMOTED COMPLEX scenarios the functions are EigsComplexXYZ().)  Could this be a compiler issue with calling those EigsRealXYZ routines and eig_val being ComplexColumnVector?  That might also explain why you are seeing three eigenvalues of 1 while my build only reports one eigenvalue of 1.

Dan Sebald <sebald>
Thu 19 Apr 2018 05:37:01 PM UTC, comment #4: 

I don't think you need to post on a bunch of other bugs, this is a relatively minor and well-contained issue.

What you can do is post whether your system is 32-bit or 64-bit, whether you are using BLAS or ATLAS or OpenBLAS with Octave, and whether you have configured Octave for 64-bit indexing or not.

Then if you want to debug this, you can step into the '__eigs__' function with gdb and see where the results you are getting diverge from another system where this test does not fail.

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 05:26:28 PM UTC, comment #3: 

The library you installed is (I think) just a wrapper around the real ARPACK library, not a reimplementation, so you should have both installed. The package you really need is libarpack2-dev, the '++' is not needed and not used by Octave.

Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 05:14:19 PM UTC, comment #2: 

Hmmmm, :-|

What library are we looking at?  When I installed all the packages needed to make Octave build, I did a search for ARPACK and picked whatever appeared in the Muon (is it?) package manager.  I see nothing there about 3.x.x series.  Instead, what I've picked is


libarpack++2-dev
  Object-oriented version of the ARPACK package (development)
libarpack++2c2a
  Object-oriented version of the ARPACK package (runtime)


Have I picked the wrong packages?

Dan Sebald <sebald>
Thu 19 Apr 2018 05:00:33 PM UTC, comment #1: 

With ARPACK 3.5.0 on Debian (Debian package version 3.5.0+real-2, same as Ubuntu dev), I get no failures:


>> test eigs
PASSES 179 out of 179 tests
>> A = magic (100) / 10 + eye (100);
>> opts.v0 = (1:100)';
>> opts.maxit = 10;
>> d = eigs (A, 10, "sm", opts);
warning: eigs: Only 3 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
>> assert (d(10), NaN+1i*NaN);
>> d(10)
ans =  NaN + NaNi
>> d
d =

   1.00000 + 0.00000i
   1.00000 - 0.00000i
   1.00000 + 0.00000i
       NaN +     NaNi
       NaN +     NaNi
       NaN +     NaNi
       NaN +     NaNi
       NaN +     NaNi
       NaN +     NaNi
       NaN +     NaNi

>> imag(d)
ans =

   0.0000e+00
  -8.6837e-13
   8.6837e-13
          NaN
          NaN
          NaN
          NaN
          NaN
          NaN
          NaN


Mike Miller <mtmiller>
Group Member
Thu 19 Apr 2018 04:50:06 PM UTC, original submission:  

I have one test failing on a very up-to-date Ubuntu 18.04 with KDE (not that the latter matters, I suppose), and it seems along the lines of this bug report. If it isn't the same, I'll create a separate bug report.

I will isolate just the failing test:


octave:1> test eigs
***** testif HAVE_ARPACK
 A = magic (100) / 10 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 10;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 10, "sm", opts);
 assert (d(10), NaN+1i*NaN);
!!!!! test failed
ASSERT errors for:  assert (d (10),NaN + 1i * NaN)

  Location  |  Observed  |  Expected  |  Reason
     ()           O            E         real != complex


and a more detailed breakdown


octave:2> A = magic (100) / 10 + eye (100);
octave:3> opts.v0 = (1:100)';
octave:4> opts.maxit = 10;
octave:5> d = eigs (A, 10, "sm", opts);
warning: eigs: Only 1 of the 10 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
octave:6> assert (d(10), NaN+1i*NaN);
error: ASSERT errors for:  assert (d (10),NaN + 1i * NaN)

  Location  |  Observed  |  Expected  |  Reason
     ()           O            E         real != complex
octave:6> d(10)
ans =  NaN
[snip]
octave:8> d
d =

   1.00000
       NaN
       NaN
       NaN
       NaN
       NaN
       NaN
       NaN
       NaN
       NaN


So, it appears that there is non-convergence but this is probably expected behavior because the test is disabling the printout of warning messages about non-convergence. It may be the only issue is that ARPACK is generating "real" (float) values where another library is generating "complex" values. Perhaps there is a change in recent versions of ARPACK, don't know. If you want me to investigate further, I can.

Dan Sebald <sebald>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #44034:  eigs.m.diff added by caliari (986B - application/x-tex)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by sebald (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-04-26 caliari StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2018-04-26 caliari StatusIn Progress Patch Submitted
    2018-04-26 caliari Attached File- Added eigs.m.diff, #44034
    2018-04-20 caliari StatusWorks For Me In Progress
    2018-04-20 mtmiller StatusNone Works For Me
    2018-04-19 sebald Attached File- Added annotation_and_F77_LOGICAL.diff, #43976

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code