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

Submitted by:  Dan Sebald <sebald>
Submitted on:  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

Add a New Comment(Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( 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>
Project 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>
Project 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>
Project Administrator
Tue 24 Apr 2018 05:39:39 PM UTC, comment #30:

The warning can be turned off with:

So the suggestion of

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"

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

John W. Eaton <jwe>
Project 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>
Project 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:

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:

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:

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>
Project 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

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?

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

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:


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

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:

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


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:

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

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.

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:

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:

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

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>
Project Administrator
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

and on Ubuntu 18.04 you have

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

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

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

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:

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

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>
Project 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

Mike Miller <mtmiller>
Project Administrator
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>
Project Administrator
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>
Project Administrator
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>
Project Administrator
Thu 19 Apr 2018 09:27:54 PM UTC, comment #8:

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


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:

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>
Project Administrator
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:

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:

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:

The NOT SYMMETRIC case has

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>
Project Administrator
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>
Project Administrator
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

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:

Mike Miller <mtmiller>
Project Administrator
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:

and a more detailed breakdown

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:

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)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project 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.4