bugGNU Octave - Bugs: bug #53719, eigs delivers wrong result in case...

 
 

bug #53719: eigs delivers wrong result in case of complex valued general eigenvalue problem

Submitter:  Michael Asam <stromi>
Submitted:  Mon 23 Apr 2018 06:47:04 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Michael Asam Open/Closed:  * Closed
Release:  * 4.2.2 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Fri 27 Apr 2018 03:20:05 PM UTC, comment #17: 

Ah, there we go (ARPACK, The Next Generation)...

https://github.com/opencollab/arpack-ng/blob/master/SRC/znaupd.f


c                    In mode 3, the vector B * X is already
c                    available in WORKD(ipntr(3)).  It does not
c                    need to be recomputed in forming OP * X.


Dan Sebald <sebald>
Fri 27 Apr 2018 06:25:14 AM UTC, comment #16: 

@Dan: that is original (old) arpack. Look for arpack-ng.

Marco Caliari <caliari>
Group Member
Thu 26 Apr 2018 04:42:00 PM UTC, comment #15: 
Dan Sebald <sebald>
Thu 26 Apr 2018 10:09:38 AM UTC, comment #14: 

I found it. which has to be set to 'LM' in the fortran code. My fault. Plain arpack "converges" to wrong eigenvalues, too. This is maybe expected, since B is not HPD. I think we should use the idea in section 3.2.2 of the documentation. I will elaborate on it. Of course, after 4.4.

@Dan: some time ago I asked for help into trying to make a review and write a sort of documentation of eigs.m->__eigs__.cc->eigs-base.cc->ARPACK, something similar to https://arxiv.org/pdf/cs/0604006. If you are interested, please privately write me.

I close the bug report.

Marco Caliari <caliari>
Group Member
Thu 26 Apr 2018 09:55:09 AM UTC, comment #13: 

@Dan: which documentation are you referring to? I cannot find it.

Marco Caliari <caliari>
Group Member
Thu 26 Apr 2018 08:53:46 AM UTC, comment #12: 

OK, the invalid number of converged Ritz makes sense.

I'm wondering about the following comment in the code:


c           %-----------------------------------------%
c           | Perform y <-- OP*x = inv[A-sigma*M]*M*x |
c           | M*x has been saved in workd(ipntr(3)).  |
c           | The user only need the linear system    |
c           | solver here that takes workd(ipntr(3))  |
c           | as input, and returns the result to     |
c           | workd(ipntr(2)).                        |
c           %-----------------------------------------%


The documentation says that when IDO=1:


c  Mode 3:  A*x = lambda*M*x, M symmetric semi-definite
c           ===> OP =  inv[A - sigma*M]*M   and  B = M.
c           ===> shift-and-invert mode
c           If OP*x = amu*x, then lambda = sigma + 1/amu.
[snip]
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.


That comment above makes it sound as though Z = M X, but then why would the documentation state Y = OP Z?  That would be Y = inv(A) M M X.  I don't follow that part of things.

Dan Sebald <sebald>
Thu 26 Apr 2018 07:13:17 AM UTC, comment #11: 

@Dan: when there is an error in znaupd/zneupd (like maximum number of iterations reached) the number of converged eigenvalues is not set in ARPACK. The final printing of my fortran file (taken by an existing driver) is meant for converged situations. So, no real problem in my fortran file or ARPACK library.

Marco Caliari <caliari>
Group Member
Thu 26 Apr 2018 06:40:22 AM UTC, comment #10: 

@Dan: sorry for the zip file, which misses debug.h and contains subfolders. I uploaded a new one. I did not notice the problem with the number of converged ritz values: so, there is probably something in my fortran file or in ARPACK itself.
In this example sigma is zero. So inv(A-sigma M) is inv(A) which is a diagonal matrix. Therefore, I divided x by the diagonal of A.

(file #44033)

Marco Caliari <caliari>
Group Member
Wed 25 Apr 2018 05:08:09 PM UTC, comment #9: 

@Marco: I've looked at the code just a little bit.  I have no concrete answers, but just observations.

First, I agree with your thinking that the FORTRAN (non-sparse construct) and C (sparse construct) implementations should behave pretty much the same in terms of convergences.  So, this is a good test.

For the FORTRAN version, I had to remove the debug header file


zndrv4_53719.f:141: Error: Can't open included file 'debug.h'


I then run the program and get some loop outputs for the residuals, then the final summary:


 _NDRV4
 ======

  Size of the matrix is           10
  The number of Ritz values requested is            4
  The number of Arnoldi vectors generated (NCV) is           10
  What portion of the spectrum: SM
  The number of converged Ritz values is    538187800
  The number of Implicit Arnoldi update iterations taken is            2
  The number of OP*x is           17
  The convergence criterion is    1.1102230246251565E-016


Everything looks reasonable there except the number of converged Ritz values reported is unrealistic.  It should report 0.  log2(538187800) is 29.004, so it's not like this could be -1 in disguise or anything.

The other thing that has me wondering is the IDO value reported keeps repeating with the pattern of 2,2,1,2,2,1,2,2,1,2,2,1...


 IDO =           -1
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2
 IDO =            1
 IDO =            2
 IDO =            2


The above does look to be about 17 repeated patterns (i.e., the number of computed OP*x).  However, I'm wondering why in your code for the IDO=-1 and IDO=1 cases there is the division by something resembling the eigenvalues or Ritz values (is it supposed to be (A-sigma*M)^-1?).  Won't that throw the estimates out of sorts?

From documentation

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

when IDO is 1, the routine wants


c          IDO =  1: compute  Y = OP * Z  and Z = B * X where


and what the code is computing is


y <-- OP*x = inv[A-sigma*M]*M*x |


Otherwise, all the FORTRAN vs. C indexing looks about right.  The mx() routine looks to match what the eigs() example is set up as.

Dan Sebald <sebald>
Wed 25 Apr 2018 10:55:26 AM UTC, comment #8: 

I run the original test (asking for 4 eigenvalues) with a plain fortran file, and arpack does not converge, which is better than silently returning the wrong eigenvalues as octave.
Arpack calls several times znaupd which returns a vector x and ask for an operation y = OP(x) (either y = B*x or y = A\x). I displayed the first element of x and y before and after each of the operations. They coincide, in plain fortran and octave, till the end of the first major iteration (about 30 calls to znaupd). Then, eigenvalues are computed. They coincide, even if they are in a different order. This is already strange, but maybe I missed a sorting operation when Arpack directly displays results. Then the first residual has to be computed. In order to do that, a B*x operation is needed. Well, at this point, x is different in plain fortran and octave. At this point, resid, which is an input/output argument for znaupd, is the same as input, but different as output. It seems that znaupd behaves differently when called by the fortran file or when called by octave. I have only one set of blas/lapack/arpack libraries on my system.

I include a zip file, with the modified eigs-base.cc which produces output, with zndrv4_53719.f (gfortran zndrv4_53719.f -lblas -llapack -larpack), and with the output files, with an entry LOOK HERE to find the divergence point. I still have no idea.

(file #44025)

Marco Caliari <caliari>
Group Member
Tue 24 Apr 2018 05:08:19 PM UTC, comment #7: 

I verified that the relative errors in the original Test_GEVP are now all zero.  Marking as fixed and closing report.

Rik <rik5>
Group administrator
Tue 24 Apr 2018 03:04:40 PM UTC, comment #6: 

Marco:  I imported your changes, made a couple of minor edits to the commit messages, then pushed on stable, and merged with default.  All the tests for eigs pass on my system.

John W. Eaton <jwe>
Group administrator
Tue 24 Apr 2018 09:04:03 AM UTC, comment #5: 

I decided for two patches. The first is about the segmentation fault. Moreover, it fixes some conditions in the subfunction select in eigs.m. It should certainly be applied (eigs1.diff).

The second changes the switch between using ARPACK or eig inside eigs.m. Unfortunately, it was not as easy as I thought in comment #3. In fact, we should compare p and rows (a). By default, p is 2 k, but it can be overwritten by opts.p. I also changed a lot of BIST tests in order to really use ARPACK and not eig. Since the default number of eigenvalues is 6, and the default Krylov dimension is 2 6, I forced eigs to use eig when rows (a) < 13 (it was rows (a) < 9). This second patch (eigs2.diff) needs an accurate review (we should be sure that all tests pass again on different platforms) and does not solve the undetected ill-conditioning of the problem in this bug report.

When eig is used, the output argument flag is set to 0. What about setting it to -0 in order to know when eig was used instead of ARPACK?

(file #44003, file #44004)

Marco Caliari <caliari>
Group Member
Mon 23 Apr 2018 03:38:05 PM UTC, comment #4: 

For question 2:

2) the diff fixes in general the crash when B is not SPD, but maybe a bug in sparse_chol is now hidden (is it correct that fact.P() does not exist when info != 0?)

This is correct.  Downstream code should not rely on any of the outputs from the factorization if the factorization has failed.

P is not, in fact a data member of the sparse_chol class, but actually a function.  The code is


      octave_idx_type P (void) const
      {
#if defined (HAVE_CHOLMOD)
        return (minor_p == static_cast<octave_idx_type>(Lsparse->ncol) ?
                0 : minor_p + 1);
#else
        return 0;
#endif
      }


So P is not calculable if Lsparse hasn't been computed.  There is an optional fourth argument that we could pass to the sparse_chol to force some kind of factorization.  Instead of this,


   octave::math::sparse_chol<SparseMatrix> fact (b, info, false);


call


   octave::math::sparse_chol<SparseMatrix> fact (b, info, false, true);


but I doubt that is a good idea.

I will commit your patch whenever you think it is sufficiently tested.

Rik <rik5>
Group administrator
Mon 23 Apr 2018 01:31:58 PM UTC, comment #3: 

The matrix B is not Hermitian positive definite, as required in the documentation, and this is the main problem here. Looking for small eigenvalues does not require explicit inversion of B by cholesky, and therefore this problem is not detected. Matlab succeeds when asking 5 eigenvalues, but returns an error with 4 (ill-conditioned problem). I believe that Matlab uses eig when asking 5 eigenvalues. This is reasonable: since the default choice is to double the number of eigenvalues for the Krylov space, when the Krylov space has the same dimension of the problem, we should use eig. This is possible by replacing


if (rows (a) - k < 3)


with


if (2 * k >= rows (a))


in eigs.m. It seems that ARPACK does not detect problems in convergence and therefore Octave silently returns wrong eigenvalues (when asking 4 eigenvalues). I will implement the example in fortran, to be 100% sure it is a not-detected-problem in ARPACK.

For the crash, there is a problem in eigs-base:make_cholb. For the sparse versions, when the factorization fails, it seems that fact.P () is not available. The test should be as in the full case


if (info != 0)
  return false;
else


fact.P() (which is the second output of chol) is probably set later, I don't know. I attach a diff (not yet a changeset): I had to change a lot of tests, otherwise call_eig was used. Morever, I fixed some conditions in eigs.m.

To summarize:

1) the diff fixes the problem by calling eig (which is reasonable for this example)
2) the diff fixes in general the crash when B is not SPD, but maybe a bug in sparse_chol is now hidden (is it correct that fact.P() does not exist when info != 0?)
3) the diff does not detect that the original problem is ill-conditioned (this is probably a defect in ARPACK, I will check).


(file #43992)

Marco Caliari <caliari>
Group Member
Mon 23 Apr 2018 08:50:47 AM UTC, comment #2: 

It seems a problem in _eigs_.cc. I will investigate.

Marco Caliari <caliari>
Group Member
Mon 23 Apr 2018 08:27:34 AM UTC, comment #1: 

Very timely bug report.  We're looking at similar issues here:

https://savannah.gnu.org/bugs/index.php?53700

(but that case is one which follows the "not sparse" path).

Just noting that these three cases work:

[Evectors Evalues] = eigs(A, B, 10, 'SM');
[Evectors Evalues] = eigs(A, B, 9, 'SM');
[Evectors Evalues] = eigs(A, B, 8, 'SM');

but going to

[Evectors Evalues] = eigs(A, B, 7, 'SM');

throws the eigenvalues out of alignment.  However, I think there is some requirement in the _eigs_.cc and/or eigs-base.cc code that the NEV should satisfy 2 <= NCV-NEV, so perhaps internally those first three above are handled with a different algorithm--don't know.

Dan Sebald <sebald>
Mon 23 Apr 2018 06:47:04 AM UTC, original submission:  

Attached I send a very small octave script showing that eigs returns wrong results in case of a general eigenvalue problem with complex matrices.

In addition octave crashes if we search for the largest eigenvalues ('LM').

I tested it in windows as well in linux. Both with the same behaviour.

Michael Asam <stromi>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #44033:  53719_2.zip added by caliari (20KiB - application/zip)
file #44025:  53719.zip added by caliari (20KiB - application/zip)
file #44003:  eigs1.diff added by caliari (3KiB - application/x-tex)
file #44004:  eigs2.diff added by caliari (11KiB - application/x-tex)
file #43992:  eigs.diff added by caliari (9KiB - application/x-tex)
file #43990:  Test_GEVP.m added by stromi (840B - text/x-objcsrc)

 

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 sebald (Posted a comment)
  • -email is unavailable- added by stromi (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 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-04-26 caliari StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2018-04-26 caliari Open/ClosedClosed Open
    2018-04-26 caliari Attached File- Added 53719_2.zip, #44033
    2018-04-25 caliari Attached File- Added 53719.zip, #44025
        StatusFixed In Progress
    2018-04-24 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2018-04-24 jwe StatusPatch Submitted In Progress
    2018-04-24 caliari Attached File- Added eigs1.diff, #44003
        Attached File- Added eigs2.diff, #44004
        StatusIn Progress Patch Submitted
    2018-04-23 rik5 Summaryeigs delivers wron result in case of complex valued general eigenvalue problem eigs delivers wrong result in case of complex valued general eigenvalue problem
    2018-04-23 caliari Attached File- Added eigs.diff, #43992
    2018-04-23 caliari StatusNone In Progress
    2018-04-23 stromi Attached File- Added Test_GEVP.m, #43990

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code