bugGNU Octave - Bugs: bug #49523, eigs (R2013b and later) uses NaNs,...

 
 

bug #49523: eigs (R2013b and later) uses NaNs, not zeros, for unconverged Ritz values

Submitter:  None
Submitted:  Thu 03 Nov 2016 01:17:22 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Tim Mitchell Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 01 Dec 2017 07:10:06 PM UTC, comment #54: 

I pushed the fix of using isreal() to determine the expected value for the specific failing BIST test.  See http://hg.savannah.gnu.org/hgweb/octave/rev/1b17a5d82e3f.

Hopefully the builbots will be happy now.

Rik <rik5>
Group administrator
Fri 01 Dec 2017 07:22:56 AM UTC, comment #53: 

I agree with Rik's proposal in comment #52: even if we can find a magic combination of v0, tol, maxit working with the existing combinations of blas/lapack/arpack/whatever, it is not sure it will work for any future combination. Here the important thing is to get either NaN or Nan+1i*NaN (and not NaN+1i*0).

To comment #37: -14 in XYeupd is "did not find any eigenvalues to sufficient accuracy". It is just non-convergence with the given parameters and it could be expected in our cases since the maximum number of iterations was strongly decreased. We could add a string to the error message like: "try to increase the tolerance, the maximum number of iterations or the number of requested eigenvalues". The third hint seems strange, but arpack uses internally Krylov spaces (spaces in which to look for eigenpairs) of default dimension opts.p = 2 * number_of_requested_eigenvalues.

Marco Caliari <caliari>
Group Member
Fri 01 Dec 2017 05:36:15 AM UTC, comment #52: 

Maybe Carlo has an idea on how to relax this.  I found that by increasing the maximum number of iterations, opts.maxit, ARPACK would converge and find all the eigenvalues.  That is not what is needed here.

I restored maxit to 1, and then started changing the initialization vector v0.  For some vectors there was rapid convergence and I got the three eigenvalues listed in comment #50 (although they were all real, no imaginary parts).

Maybe there is a specific initial vector which is always far enough from a solution that the answer from comment #46 is obtained.

Otherwise, we could work around it by coding the BIST test like this


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



Rik <rik5>
Group administrator
Fri 01 Dec 2017 05:00:41 AM UTC, comment #51: 

The last representation in comment #50 is the best.  The imaginary component is very small, but non-zero, so Octave can't discard it.

Rik <rik5>
Group administrator
Fri 01 Dec 2017 01:51:17 AM UTC, comment #50: 

This is probably the best (format long e):

d =

   1.69714256958507e-02 + 2.37387114336600e-17i
   1.69714256958507e-02 - 2.37387114336600e-17i
   6.79928617498172e-02 + 0.00000000000000e+00i
                    NaN +                  NaNi

Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 01:49:02 AM UTC, comment #49: 

This is with "format bit"



d =

  0011111110010001011000001111010100011001010100100001100011001000  0011110001111011010111100110111001010010110100110001001111010011i
  0011111110010001011000001111010100011001010100100001100011001000  1011110001111011010111100110111001010010110100110001001111010011i
  0011111110110001011001111111101011101101100100111011000010100111  0000000000000000000000000000000000000000000000000000000000000000i
  0111111111111000000000000000000000000000000000000000000000000000  0111111111111000000000000000000000000000000000000000000000000000i


Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 01:45:05 AM UTC, comment #48: 

octave:8> d = eigs (A, 4, "sm", opts)
warning: eigs: Only 3 of the 4 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
d =

   0.0169714256958507 + 0.0000000000000000i
   0.0169714256958507 - 0.0000000000000000i
   0.0679928617498172 + 0.0000000000000000i
                  NaN +                NaNi

Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 01:43:47 AM UTC, comment #47: 

Well that is the problem.  Could you redo the test with openblas (comment #45) and use 'format long' before displaying the variable d?  It looks like Octave couldn't narrow the variable from complex to real because there was a very small imaginary value.


Rik <rik5>
Group administrator
Fri 01 Dec 2017 01:18:45 AM UTC, comment #46: 

With atlas it is the same as on i7 (with any blas):

LD_PRELOAD=/usr/lib64/atlas/libtatlas.so ./run-octave -f -q --no-gui
octave:1> test eigs
PASSES 179 out of 179 tests
octave:2>  A = toeplitz ([0, 1, zeros(1, 8)], [0, -1, zeros(1, 8)]);
octave:3>  A(1, 1) = 1;
octave:4>  A = kron (A, eye (10)) + kron (eye (10), A);
octave:5>  opts.v0 = (1:100)';
octave:6>  opts.maxit = 1;
octave:7> d = eigs (A, 4, "sm", opts)
warning: eigs: Only 2 of the 4 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
d =

   0.016971
   0.016971
        NaN
        NaN


Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 01:13:07 AM UTC, comment #45: 

With openblas:

octave:2> A = toeplitz ([0, 1, zeros(1, 8)], [0, -1, zeros(1, 8)]);
octave:3>  A(1, 1) = 1;
octave:4>  A = kron (A, eye (10)) + kron (eye (10), A);
octave:5>  opts.v0 = (1:100)';
octave:6>  opts.maxit = 1;
octave:7> d = eigs (A, 4, "sm", opts)
warning: eigs: Only 3 of the 4 requested eigenvalues converged
warning: called from
    eigs at line 265 column 18
d =

   0.016971 + 0.000000i
   0.016971 - 0.000000i
   0.067993 + 0.000000i
        NaN +      NaNi

Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 01:09:42 AM UTC, comment #44: 

If you could run interactively and get the value for d(4) that would help.

In the message from assert, the 'O' just stands for observed and the 'E' stands for expected.  The message 'complex != real' explains what is going on.  The two asserts are


assert (d(4), NaN);
assert (imag (d(4)), 0);


and it is failing on the first one with


ASSERT errors for: assert (d (4),NaN)

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


The expected is just NaN, a real number.  The observed was a complex value.  You can see that the second test is designed to ensure that d(4) is a real value by checking that the imaginary part was equal to 0.  Perhaps a better test would be to check directly with isreal.


assert (isreal (d(4)))



Rik <rik5>
Group administrator
Fri 01 Dec 2017 12:41:03 AM UTC, comment #43: 

I meant to say openblas is 0.2.20 (the latest stable release).

Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Fri 01 Dec 2017 12:38:17 AM UTC, comment #42: 

This erros happens with either openblas-0.2.0 (custom compiled for this machine) or reference blas ("standard" Fedora package).
No problem with atlas.
The build bot is AMD FX-8350; my workstation is i7-2600 and does
not show this error (openblas is different , but reference
blas is exactly the same).

I am not sure what to do about it. I can set buildbot to
use atlas.

On related note -- the error seems strange to me. Should not
observed and expected values be numbers rather than letters?

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 11:28:53 PM UTC, comment #41: 

A-a-a-nd I see that fedora buildbot still fails with:
http://buildbot.octave.org:8010/builders/clang-fedora/builds/1015/steps/test/logs/stdio


 ***** testif HAVE_ARPACK
 A = toeplitz ([0, 1, zeros(1, 8)], [0, -1, zeros(1, 8)]);
 A(1, 1) = 1;
 A = kron (A, eye (10)) + kron (eye (10), A);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(4), NaN);
 assert (imag (d(4)), 0);
!!!!! test failed
ASSERT errors for:  assert (d (4),NaN)

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

it could be related to openblas library it uses. I will investigate.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 10:54:36 PM UTC, comment #40: 

I pushed Marco's relaxed BIST test changeset here: http://hg.savannah.gnu.org/hgweb/octave/rev/30bbef7bf919.

Closing report again.

Rik <rik5>
Group administrator
Thu 30 Nov 2017 08:59:08 PM UTC, comment #39: 

I fixed the incorrect error reporting here http://hg.savannah.gnu.org/hgweb/octave/rev/30240b146061.

I will review eigs3.diff and push that shortly as well.

Rik <rik5>
Group administrator
Thu 30 Nov 2017 08:46:43 PM UTC, comment #38: 

Yes, with eigs3.diff eigs passes all tests.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 08:45:52 PM UTC, comment #37: 

@Marco: I was just making a guess.  Since eigs is returning NaN for sum of the eigenvalues I assume that some of the eigenvectors must also have NaN in them and that maybe this was upsetting zneupd.  I suppose the real way to figure out the problem is to go look at the Fortran code and figure out what the error code of -14 means.

I will take care of changing the error messages to match the Fortran call involved.

Rik <rik5>
Group administrator
Thu 30 Nov 2017 08:33:06 PM UTC, comment #36: 

@Rik: not sure I understood your comments. znaupd quits with a value of ido in {-1,1,2,99}. 99 means "done" and info is checked. Of course the error message is wrong. Which vector could have norm NaN before zneupd?

Marco Caliari <caliari>
Group Member
Thu 30 Nov 2017 07:24:11 PM UTC, comment #35: 

@Dmitri: Can you test Marco's patch on your failing system?  I don't want to commit it to Mercurial until it is clear that it solves the problem.

@Marco: Are we sure that there is no way for the code in eigs-base.cc to exit from znaupd with an error code that is not checked before calling zneupd?  Brief documentation for zneupd is at http://www.caam.rice.edu/software/ARPACK/UG/node44.html#2332.

According to that, "On the final return from znaupd (indicated by ido = 99), the error flag info must be checked. If info = 0, then no fatal errors have occurred and it is time to post-process using zneupd  to get eigenvalues of the original problem and the corresponding eigenvectors if desired."

The error reporting looks wrong too.


          if (info < 0)
            (*current_liboctave_error_handler)
              ("eigs: error %d in dsaupd", info);  // line 3371


This looks like it should reference znaupd, not dsaupd.

Also, what happens if the norm of one of the vectors is NaN and it attempts to normalize the lengths in zneupd?  That might throw an error.

Rik <rik5>
Group administrator
Thu 30 Nov 2017 07:07:30 PM UTC, comment #34: 

Please remove eigs2.diff and apply eigs3.diff.

@Rik: in fact, I relaxed the tests.

(file #42536)

Marco Caliari <caliari>
Group Member
Thu 30 Nov 2017 06:40:12 PM UTC, comment #33: 

OK, so with  the patch is only one test is failing, although the
errors are different depending on blas version.

Ref. blas:

LD_PRELOAD=/usr/lib64/libblas.so ./run-octave -f -q --no-gui
octave:1> test ("eigs", "normal", "results_blas_1.txt")

>>>>> /d2/home/dima/src/octave/gcc_def/../scripts/sparse/eigs.m

PASSES 178 out of 179 tests
octave:2> test ("eigs", "normal", "results_blas_2.txt")

>>>>> /d2/home/dima/src/octave/gcc_def/../scripts/sparse/eigs.m

PASSES 178 out of 179 tests
octave:3> test eigs
*** testif HAVE_ARPACK
 A = 1i * magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
eigs: error -14 in zneupd


Same with openblas:

LD_PRELOAD=/usr/lib64/libopenblaso.so ./run-octave -f -q --no-gui
octave:1> test ("eigs", "normal", "results_openblas_1.txt")

>>>>> /d2/home/dima/src/octave/gcc_def/../scripts/sparse/eigs.m

PASSES 178 out of 179 tests
octave:2> test ("eigs", "normal", "results_openblas_2.txt")

>>>>> /d2/home/dima/src/octave/gcc_def/../scripts/sparse/eigs.m

PASSES 178 out of 179 tests
octave:3> test eigs
*** testif HAVE_ARPACK
 A = 1i * magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
eigs: error -14 in zneupd

Same test,  different error with atlas:

 LD_PRELOAD=/usr/lib64/atlas/libtatlas.so ./run-octave -f -q --no-gui
octave:1> test eigs
*** testif HAVE_ARPACK
 A = magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 10;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(4), NaN+1i*NaN);
!!!!! test failed
ASSERT errors for:  assert (d (4),NaN + 1i * NaN)

  Location  |  Observed  |  Expected  |  Reason
     ()      1-1.219e-14i   NaN+NaNi     'NaN' mismatch
octave:2> test ("eigs", "normal", "results_atlas_1.txt")

>>>>> /d2/home/dima/src/octave/gcc_def/../scripts/sparse/eigs.m

PASSES 178 out of 179 tests

Dmitri.
--


Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 05:07:23 PM UTC, comment #32: 

@Marco: Is it possible to pick the runtime parameters such as the initialization vector v0 and the maximum number of iterations is such that at least some of the returned values are NaN?  I know you have initialized both of those items, but I wonder whether we should attempt to optimize those, or maybe we should relax the tests.  For example, this is a failing test


A = toeplitz ([-2, 1, zeros(1, 8)]);
 A = kron (A, eye (10)) + kron (eye (10), A);
 opts.v0 = (1:100)';
 opts.maxit = 3;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "lm", opts);
 assert (d(2:4), [NaN; NaN; NaN]);
!!!!! test failed
ASSERT errors for:  assert (d (2:4),[NaN; NaN; NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)         -7.365        NaN        'NaN' mismatch


What if we relaxed the conditions for the assert so that it was


assert (d(3:4), [NaN; NaN]);


This would still verify that unconverged Ritz values are returned as NaN, but wouldn't require that exactly 3 Ritz values be unconverged.

Rik <rik5>
Group administrator
Thu 30 Nov 2017 04:54:24 PM UTC, comment #31: 

@Dmitri: You should be able to run all of the tests and dump the log into a file with


test ("eigs", "normal", "results.txt")


This will stick the logfile into the filename "results.txt".


Rik <rik5>
Group administrator
Thu 30 Nov 2017 04:14:15 PM UTC, comment #30: 

Before your patch 6 eigs test were failing on Fedora and 5 on debian (jwe buildbot). You can see logs at
 
http://buildbot.octave.org:8010/waterfall

There are still multiple failures (some appear vary depending on
blas version used), but I could not figure out how to get
"test" to report all the failures in interactive mode.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 04:02:31 PM UTC, comment #29: 

@Dmitri: Any other test not working?

Marco Caliari <caliari>
Group Member
Thu 30 Nov 2017 03:59:33 PM UTC, comment #28: 

octave:1> A = 1i * magic (100) + eye (100);
octave:2> opts.v0 = (1:100)';
octave:3> opts.maxit = 7;
octave:4> 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
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

(Same with either ref blas, atlas, openblas)

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 03:52:26 PM UTC, comment #27: 

@Dmitri: can you tell me the result of


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


But is it worth/possible to try to construct a set of examples working (that is, with eigs failing) with any possible combination of blas/lapack/arpack libraries?

Marco Caliari <caliari>
Group Member
Thu 30 Nov 2017 02:57:02 PM UTC, comment #26: 

The patch changed some errors, e.g. with the reference blas
I am getting this now:

test eigs
*** testif HAVE_ARPACK
 A = 1i * magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
eigs: error -14 in zneupd

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 08:06:58 AM UTC, comment #25: 

@Rik, Dmitri: are the results really intermittent or 100% reproducible?
In the first case, no idea, since in the tests I initialized the initial vector which should be the only otherwise random choice.
I tried standard blas, atlas and openblas on a Ubuntu systems, I was able to reproduce some failures. Convergence clearly depends on the blas/lapack libraries (and the very stringent tolerance, which is eps by default). I relaxed the tests and now they succeed with all my libraries. Please try the patch.

(file #42528)

Marco Caliari <caliari>
Group Member
Thu 30 Nov 2017 05:53:33 AM UTC, comment #24: 

I tried arpack 3.4 (from Fedora 26) arpack 3.6 (master source from https://github.com/opencollab/arpack-ng).

Same results.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 05:33:43 AM UTC, comment #23: 

I also get different test failures depending on the blas version used:

LD_PRELOAD=/usr/lib64/libopenblas.so.0 ./run-octave -f -q --no-gui
octave:1> test eigs
*** testif HAVE_ARPACK
 A = toeplitz ([-2, 1, zeros(1, 8)]);
 A = kron (A, eye (10)) + kron (eye (10), A);
 opts.v0 = (1:100)';
 opts.maxit = 3;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "lm", opts);
 assert (d(2:4), [NaN; NaN; NaN]);
!!!!! test failed
ASSERT errors for:  assert (d (2:4),[NaN; NaN; NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)         -7.365        NaN        'NaN' mismatch

LD_PRELOAD=/usr/lib64/atlas/libsatlas.so ./run-octave -f -q --no-gui
octave:1> test eigs
*** testif HAVE_ARPACK
 A = magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
ASSERT errors for:  assert (d (3:4),[NaN + 1i NaN; NaN + 1i NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)      1+1.783e-14i   NaN+NaNi     'NaN' mismatch

LD_PRELOAD=/usr/lib64/libblas.so.3 ./run-octave -f -q --no-gui
octave:1> test eigs
*** testif HAVE_ARPACK
 A = toeplitz ([0, 1, zeros(1, 8)], [0, -1, zeros(1, 8)]);
 A(1, 1) = 1;
 A = kron (A, eye (10)) + kron (eye (10), A);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN; NaN]);
 assert (imag (d(3:4)), [0; 0]);
!!!!! test failed
ASSERT errors for:  assert (d (3:4),[NaN; NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)        0.067993       NaN        'NaN' mismatch

Those are 100% reproducible.
I use
arpack-devel-3.5.0-4.fc27.x86_64
arpack-3.5.0-4.fc27.x86_64
(same as fedora buildbots)

Dmitri.

Dmitri A. Sergatskov <dasergatskov>
Thu 30 Nov 2017 05:22:15 AM UTC, comment #22: 

Have to re-open the bug temporarily.  The BIST tests are creating intermittent failures on the automated buildbot testing system.  Dmitri first reported this.  See log files at http://buildbot.octave.org:8010/waterfall.

Some sample errors are:


***** testif HAVE_ARPACK
 A = toeplitz ([-2, 1, zeros(1, 8)]);
 A = kron (A, eye (10)) + kron (eye (10), A);
 opts.v0 = (1:100)';
 opts.maxit = 3;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "lm", opts);
 assert (d(2:4), [NaN; NaN; NaN]);
!!!!! test failed
ASSERT errors for:  assert (d (2:4),[NaN; NaN; NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)         -7.365        NaN        'NaN' mismatch
***** testif HAVE_ARPACK
 A = toeplitz ([-2, 1, zeros(1, 8)]);
 A = kron (A, eye (10)) + kron (eye (10), A);
 Afun = @(x) A * x;
 opts.v0 = (1:100)';
 opts.maxit = 3;
 opts.issym = true;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (Afun, 100, 4, "sm", opts);
 assert (d(2:4), [NaN; NaN; NaN]);
!!!!! test failed
ASSERT errors for:  assert (d (2:4),[NaN; NaN; NaN])

  Location  |  Observed  |  Expected  |  Reason
    (1)        -0.13578       NaN        'NaN' mismatch
***** testif HAVE_ARPACK
 A = magic (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "lm", opts);
 assert (d(4), NaN);
!!!!! test failed
ASSERT errors for:  assert (d (4),NaN)

  Location  |  Observed  |  Expected  |  Reason
     ()       3.8548e-11      NaN        'NaN' mismatch
***** testif HAVE_ARPACK
 A = magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(3:4), [NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
eigs: error -14 in dneupd
***** testif HAVE_ARPACK
 A = 1i * magic (100) / 100 + eye (100);
 opts.v0 = (1:100)';
 opts.maxit = 1;
 warning ("off", "Octave:eigs:UnconvergedEigenvalues", "local");
 d = eigs (A, 4, "sm", opts);
 assert (d(2:4), [NaN+1i*NaN; NaN+1i*NaN; NaN+1i*NaN]);
!!!!! test failed
eigs: error -14 in zneupd


This is stressing the ARPACK library.  And maybe it indicates the ARPACK library used by the buildbots is bad.

Rik <rik5>
Group administrator
Wed 29 Nov 2017 09:12:18 PM UTC, comment #21: 

I checked in Marco's changeset on the development branch here (http://hg.savannah.gnu.org/hgweb/octave/rev/0c1417a8d9d0).

@Marco: You were right that it isn't a good idea to disable all warnings during the BIST tests.  Instead, I created a new warning ID, "Octave:eigs:UnconvergedEigenvalues", and I temporarily disable just this one warning during the BIST tests.  See cset http://hg.savannah.gnu.org/hgweb/octave/rev/d3d38a443df8.

Marking report as fixed and closing it.

Rik <rik5>
Group administrator
Mon 27 Nov 2017 02:35:13 PM UTC, comment #20: 

Great, thank you so much for addressing the issue I reported!

Anonymous
Mon 27 Nov 2017 02:04:10 PM UTC, comment #19: 

I changed my mind and implemented the NaN/NaN+1i*NaN version. Unconverged eigenvalues and eigenvectors are always in the last positions. I used warning("off","all") in the tests, but maybe there is a better way to suppress the specific warning abount unconverged values.

(file #42497)

Marco Caliari <caliari>
Group Member
Fri 24 Nov 2017 05:34:09 PM UTC, comment #18: 

@Rik: I'm going to implement the NaN + i*0 version for both eigenvalues and eigenvectors. Now that eigs in octave is going to be in a good shape, I discovered that matlab switched to "A Krylov Schur algorithm for large eigenproblems". I found only a simple implementation around (https://github.com/dingxiong/KrylovSchur), so I think that matlab implemented it from scratch.

Marco Caliari <caliari>
Group Member
Thu 23 Nov 2017 06:30:30 PM UTC, comment #17: 

@Marco: I had originally proposed setting just the real part to NaN.  But, I agree with you that it looks nicer if NaN + NaNi is returned for complex unconverged eigenvalues.

I don't think it matters from a performance standpoint if we post-process the list of eigenvalues because the number is going to be small, say 50 eigenvalues, relative to the size of the sparse matrix itself which could have thousands to millions of rows or columns.

I'll accept whatever you feel is best.

Rik <rik5>
Group administrator
Thu 23 Nov 2017 03:57:45 PM UTC, comment #16: 

@Rik: Matlab now gives NaN+0.0000i, if there are other complex eigenvalues. Is it what you proposed in comment #12?

Marco Caliari <caliari>
Group Member
Thu 23 Nov 2017 03:42:04 PM UTC, comment #15: 

@Marco: Any progress?  This bug seems easy enough to knock off now that the behavior (return NaN, or Nan+Nani, for unconverged Ritz values and full NaN vectors for corresponding eigenvectors) has been agreed upon.

Rik <rik5>
Group administrator
Tue 14 Nov 2017 05:53:18 AM UTC, comment #14: 

@Marco: I just checked in 3 of your csets improving the behavior of eigs().  This bug report seems like a natural one to finish off as well.  I think you can code Octave to return Nan+NaNi for unconverged values.

Rik <rik5>
Group administrator
Mon 17 Apr 2017 03:50:10 AM UTC, comment #13: 

@Marco: I'm going to change my advice in comment #12.  It seems to me that eigs is generally only used when you want a few eigenvalues from a large sparse matrix.  For a low number, say N < 50, it really isn't much of a performance hit to iterate over the entire return list to see if there are any true complex eigenvalues or if there are only complex eigenvalues where the imaginary part is zero.  Thus, I think you should go ahead and return NaN if all the values are real or Nan+NaNi if there are complex values.

Rik <rik5>
Group administrator
Wed 12 Apr 2017 04:31:01 PM UTC, comment #12: 

The unconverged Ritz values represent an exceptional case.  It is important to warn about them, but the exact form of the warning is not that critical.  My suggestion was to set only the real portion to NaN.  After that, I would rely on Octave's automatic narrowing to change a complex matrix, with all imaginary parts equal to zero, to a real matrix.  If the matrix can't be narrowed because there are non-zero complex values then you do get


NaN + 0.00000i


but that is okay for me since the NaN stands out as the first thing in the entry, and isnan returns true if either the real or imaginary part is NaN.


Rik <rik5>
Group administrator
Wed 12 Apr 2017 07:37:11 AM UTC, comment #11: 

@Rik: If I set them to Complex (NaN, 0.), the result will be somethinkg like (there are other complex eigenvalues)


ans =

  -0.00000 - 3.32975i
  -0.00000 + 3.32975i
       NaN + 0.00000i


or (all the remaining eigenvalues are real)


ans =

   3.2154
      NaN
      NaN


If I set them to Complex (NaN, NaN), the result will be somethinkg like


ans =

  -0.00000 - 3.32975i
  -0.00000 + 3.32975i
       NaN + NaNi


or


ans =

   3.2154 + 0.00000i
      NaN + NaNi
      NaN + NaNi


If we want


ans =

  -0.00000 - 3.32975i
  -0.00000 + 3.32975i
       NaN + NaNi


or


ans =

   3.2154
      NaN
      NaN


we have to check the imaginary part of the converged eigenvalues. I mean, it is not possible to set "without the complex [imaginary] part".

Marco Caliari <caliari>
Group Member
Tue 11 Apr 2017 03:15:01 PM UTC, comment #10: 

I would just set the unconverged values to NaN without the complex part.  The real part being NaN is enough to indicate that there is something wrong with this particular value.

Rik <rik5>
Group administrator
Tue 11 Apr 2017 07:52:57 AM UTC, comment #9: 

In the nonsymmetric case, the nonconverged Ritz values are set to NaN+NaNi in Matlab. You can try


A = toeplitz ([-2, 1, zeros(1, 8)]) + magic (10) * 1e-10;
A = kron (A, eye (10)) + kron (eye (10), A);
opts.v0 = (1:100)';
opts.maxit = 1;
eigs (A, 8, 'sm', opts)


The converged values are complex with imaginary part zero. If you omit opts, you see that all the converged eigenvalues are double. If we follow Matlab and comment #5, in case of nonconverged Ritz values we should set them and the corresponding eigenvectors to NaN+NaNi. An alternative would be: if the converged Ritz values are double, set the nonconverged and the corresponding eigenvectors to (double) NaN. This would save some memory.

Marco Caliari <caliari>
Group Member
Mon 10 Apr 2017 05:28:53 PM UTC, comment #8: 

@Rik: Interesting.  On R2016b, on my hardware, I consistently see that allocating NaNs is about 10 times slower than just zeros.  Allocating ones is also about 10 times slower than just zeros on my hardware.

In any case, I don't feel strongly either way on how to handle the values of the "non-converged" eigenvectors and I do generally agree with trying to prevent garbage in, garbage out.  I was just pointing out that there was a potential downside to explicitly setting all the entries in the eigenvectors.

Anonymous
Mon 10 Apr 2017 04:46:22 PM UTC, comment #7: 

I don't see any real difference in creating matrices with all zeros or all NaNs.  See below


octave:1> tic; x = zeros (1e3,1e3); toc
Elapsed time is 0.00626612 seconds.
octave:2> tic; x = zeros (1e3,1e3); toc
Elapsed time is 0.00639701 seconds.
octave:3> tic; x = NaN (1e3,1e3); toc
Elapsed time is 0.00226998 seconds.
octave:4> tic; x = NaN (1e3,1e3); toc
Elapsed time is 0.00632405 seconds.
octave:5> tic; x = NaN (1e3,1e3); toc
Elapsed time is 0.00363517 seconds.


I don't think performance is an issue here.  And I have always tried to code Octave to avoid Garbage In / Garbage Out phenomena.  I believe many Octave users are professional scientists first, and coders second.  I don't think they will necessarily realize that they need to code defensively and check their eigenvalues before going on to use the eigenvectors.  For that reason, I would prefer to return NaNs for the eigenvectors to immediately inform the coder that something is amiss.

Rik <rik5>
Group administrator
Mon 10 Apr 2017 04:22:19 PM UTC, comment #6: 

Original submitter here. 

@Marco: The MathWorks has their own private fork of the original ARPACK source with various changes of their own, which probably explains why you are seeing different behaviors regarding the eigenvectors of non-converged Ritz values in Matlab compared to ARPACK and arpack-ng.

@Marco, @Rik: an argument against using nans for eigenvecotors whose Ritz values did not converge is that returning a large number of nans will take longer than just returning zeros (or using whatever values are already set for these eigenvectors).  On my laptop, setting up an array of 1 million nans takes about 10 times longer than 1 million zeros.  Perhaps it's a negligible cost but I personally don't see that it's necessary to use nans.  The user already has to check which, if any, eigenvalues are nans.  At that point, it should be clear to the user which eigenvectors are good and which are to be ignored.  That said, I could see a user using an eigenvector without checking if its Ritz value converged!

In my eigsPlus package, I avoided this issue altogether by only returning the converged eigenvalues (and eigenvectors).  ;-)  I believe I also suggested this to The MathWorks long ago but they didn't want to change that.

Anonymous
Mon 10 Apr 2017 04:01:24 PM UTC, comment #5: 

@Marco: I agree with you.  If one of the eigenvalues is NaN, it seems like the eigenvector for that eigenvalue should also be NaN.

Rik <rik5>
Group administrator
Mon 10 Apr 2017 02:59:48 PM UTC, comment #4: 

I checked that ARPACK gives the zero eigenvector for not converged Ritz values. I cannot understand the result in Matlab. For non-converged Ritz values I would return a NaN eigenvector. And I will propose this change to arpack-ng.

Marco Caliari <caliari>
Group Member
Fri 07 Apr 2017 12:29:47 PM UTC, comment #3: 

For the same example below, if you call


[V, D] = eigs (A, 4, 'sm', opts)


in Matlab, you get D with the entry NaN, but V with four columns of numbers. Should we copy such a weird result?

Marco Caliari <caliari>
Group Member
Fri 07 Apr 2017 08:55:17 AM UTC, comment #2: 

I created a reproducible example for this issue.


A = toeplitz ([-2, 1, zeros(1, 8)]);
A = kron (A, eye (10)) + kron (eye (10), A);
opts.v0 = (1:100)';
opts.maxit = 1;
eigs (A, 4, "sm", opts)
warning: eigs: Only 3 of the 4 requested eigenvalues converged
warning: called from
    eigs at line 266 column 18
ans =

   0.00000
  -0.16203
  -0.39851
  -0.39851


Marco Caliari <caliari>
Group Member
Thu 03 Nov 2016 03:38:05 PM UTC, comment #1: 

This seems like a good idea.  The Octave devs are working really hard right now to get the 4.2.0 release out, but after that this could be taken up on the development branch and made a part of the 4.4 release.

Rik <rik5>
Group administrator
Thu 03 Nov 2016 01:17:22 PM UTC, original submission:  

Since R2013b, Matlab's version of eigs uses NaNs as placeholders for the Ritz values, not zeros, which was the older convention and is what Octave appears to still be using (I tested 4.0.0). 

The reason for the change to NaNs is under the older scheme using zeros, if not all the requested eigenvalues converged, it could be difficult to sort out the converged eigenvalues from the placeholder zeros, especially if 0 should be an eigenvalue of the matrix in question.  Furthermore, if a user is interested in computing the spectral abscissa of a stable matrix, they may try calling

 max(real(eigs(A,6,'LR')))

but this would misleading return 0 in the case when not all requested eigenvalues converge. 

In fact, I am the person who pointed these issues out to The MathWorks in the first place and at the time, I made the suggestion for them to use NaNs instead, after I had created my own modified version of eigs doing just that.  It would be great if Octave could adopt this convention as well.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #42536:  eigs3.diff added by caliari (3KiB - application/x-tex)
file #42528:  eigs2.diff added by caliari (2KiB - application/x-tex)
file #42497:  eigs.diff added by caliari (20KiB - 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 dasergatskov (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by None (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 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-11-30 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2017-11-30 caliari Attached File- Added eigs3.diff, #42536
    2017-11-30 caliari Attached File- Added eigs2.diff, #42528
    2017-11-30 rik5 Open/ClosedClosed Open
    2017-11-30 rik5 StatusFixed In Progress
    2017-11-29 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2017-11-27 caliari StatusConfirmed Patch Submitted
    2017-11-27 caliari Attached File- Added eigs.diff, #42497
    2017-04-10 rik5 Release4.0.0 dev
    2016-11-03 rik5 StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code