bugGNU Octave - Bugs: bug #47236, output of eigs() with 'sa' option...

 
 

bug #47236: output of eigs() with 'sa' option sorted differently from Matlab

Submitted by:  None
Submitted on:  Tue 23 Feb 2016 03:14:27 PM UTC  
 
Category: Octave FunctionSeverity: 2 - Minor
Priority: 5 - NormalItem Group: Matlab Compatibility
Status: FixedAssigned to: None
Originator Name: Oliver SampsonOriginator Email: -unavailable-
Open/Closed: ClosedRelease: 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 Jump to the original submission)

Tue 14 Nov 2017 05:46:38 AM UTC, comment #18:

This bug was fixed by one of the csets for bug #45153. Marking as fixed and closig report.

Rik <rik5>
Project Administrator
Tue 18 Apr 2017 03:10:42 PM UTC, comment #17:

I prepared a script in order to test all the cases. When sigma is a scalar, the 'sm' mode is used, so it is enough to fix those cases. Moreover, the functional forms should behave like the matrix forms.
The script passes on Matlab R2014b. I fixed the initial vector in order to have reproducible results, but different initial vectors can lead to different orders (you can try it by removing the input opts). Results are different when asking only the eigenvalues or the eigenvectors, too. ARPACK should sort the eigenvalues someway, but for the operator OP, which can be A or inv(A) or inv(A-sigma*I). And Matlab can sort the result, this I do not know. Anyway, I commented some weird results. For instance, for the real symmetric 'lm' case, the eigenvalues are not sorted when the eigenvectors are requested. In the complex cases, it is rather the other way round. In Octave eigs, there are already some ordering algorithms, as it is clear from comment #13. I think they were introduced to (try to) emulate Matlab's behavior. The script does not pass on Octave. With my patch in comment #14, it passes the real cases. The final comment is that Octave uses (or should use) arpack-ng, while it is not clear what Matlab does (maybe a modified version of arpack-ng?).

I see three possibilities:
1) give back the eigenvalues exactly as given by ARPACK, without any reordering.
2) try to fix all the cases in order this script to pass. But, because of the initial vector, because of possibly different ARPACK versions, we will be never sure to be Matlab compatible.
3) sort all the eigenvalues given by ARPACK using the subfunction select already in eigs.m (and therefore remove any attempt to sort eigenvalues inside eigs-base.cc).

I think 3) is the most reliable choice.

(file #40438)

Marco Caliari <caliari>
Project Member
Wed 12 Apr 2017 06:14:57 PM UTC, comment #16:

@Rik: I think the other cases are fine. Anyway, I'm preparing the script as you suggested.

Marco Caliari <caliari>
Project Member
Wed 12 Apr 2017 04:41:00 PM UTC, comment #15:

@Marco: What about my suggestion in comment #13:

"I think what should happen is a script gets written that does the full combinatorial test of all matrix types against all solver options ('lm', 'sm', 'sa', 'la' etc.). Run that script on Octave and Matlab and see where the differences are."

I see that you took my suggestion and made a special exception for "SA". Are there other cases that we need to add? Without a full test it seems like eventually someone is going to come along and say it doesn't work with option 'li'. We will fix that in a different bug report. And then someone will file a bug that it doesn't work with 'si'. This will go on and on until we have eventually tested all options.

Rik <rik5>
Project Administrator
Wed 29 Mar 2017 06:44:43 PM UTC, comment #14:

See bug #45153. For this specific problem, I basically used the patch by Rik in comment #13.

Marco Caliari <caliari>
Project Member
Sun 28 Feb 2016 01:39:25 AM UTC, comment #13:

It turns out this is very specific to the type of matrix being solved. The code for eigs is eventually located in liboctave/numeric/eigs-base.cc. The solvers are templates based on qualities of the matrix. It turns out that only two of the matrix types are re-ordering the eigenvalues. One of these cases is the real symmetric which is what matrix A is.

The code is around line 927.

If I alter the code to exclude the "SA" option

then the results returned are in the same order as eig.

I'm attaching the quick patch I made (eigs_sa.patch) and the test script tst_eigs.m which depends on the file yaledata.

I think what should happen is a script gets written that does the full combinatorial test of all matrix types against all solver options ('lm', 'sm', 'sa', 'la' etc.). Run that script on Octave and Matlab and see where the differences are.

(file #36496, file #36497)

Rik <rik5>
Project Administrator
Sun 28 Feb 2016 01:24:44 AM UTC, comment #12:

Changing the category to Matlab Compatibility rather than Innacurate Result. The result is correct, just presented differently than Matlab.

Rik <rik5>
Project Administrator
Sat 27 Feb 2016 12:44:00 PM UTC, comment #11:

"In which sense -3.5256 is smaller than -3.9126? "

It is closer to zero. Its abs is smaller.

This all started from a professor's demo, written to run in matlab, which does not work in Octave. I think we should attempt to make eigs work like matlab's eigs so that students can use octave instead of matlab.

Doug Stewart <dastew>
Thu 25 Feb 2016 08:32:35 AM UTC, comment #10:

@ Doug, comment #7

In which sense -3.5256 is smaller than -3.9126? So, the order in Octave is consistent between A and L and inverse w.r.t. to Matlab (in both cases).

Anyway, please take a look to #46683, in particular comment #17.

Marco

Marco Caliari <caliari>
Project Member
Tue 23 Feb 2016 10:21:52 PM UTC, comment #9:

Yes I see what Matlab is saying about sorting eigenvalues and vectors. But if we ask for the 3 smallest then is it not sorted to find the 3 smallest? And if it is already sorted to find them then it should be sorted in the output.
A bigger question why sorted one way for the A matrix and the other way for the L matrix!

Doug Stewart <dastew>
Tue 23 Feb 2016 08:32:13 PM UTC, comment #8:

Thanks for confirming what I observed. So the returned eigenvalues are the correct values, they are just in the wrong sort order. And Matlab explicitly states that the diagonal matrix of eigenvalues may not always be in sorted order.

Mike Miller <mtmiller>
Project Administrator
Tue 23 Feb 2016 08:16:59 PM UTC, comment #7:

d =

Diagonal Matrix

2.5453e-02 0 0
0 2.0710e-02 0
0 0 -3.2499e-15

>> lapsh1


m =

Diagonal Matrix

-3.5256 0 0
0 -3.5358 0
0 0 -3.9126
from m smallest number is top left but with d the smallest number is at bottom right.

here is one more example of d

d =

Diagonal Matrix

2.5453e-02 0 0
0 2.0710e-02 0
0 0 4.4038e-15

"sa"
Smallest Algebraic (valid only for real symmetric problems).

matlab always does it with smallest at top left. That is what I was told. and that is what the OP was expecting and the program was expecting.

Doug Stewart <dastew>
Tue 23 Feb 2016 07:58:08 PM UTC, comment #6:

Can you show what the results are and what they should be? I'm assuming it's the sort order of the d diagonal matrix, but the values are correct?

Note that even Matlab's help for eigs claims that "eigs does not always return sorted eigenvalues and eigenvectors. Use sort to explicitly sort the output eigenvalues and eigenvectors in cases where their order is important." (from https://www.mathworks.com/help/matlab/ref/eigs.html)

Mike Miller <mtmiller>
Project Administrator
Tue 23 Feb 2016 06:20:59 PM UTC, comment #5:

I have narrowed it down some more.

clear
load yaledata;

[v, m] = eigs (A, 3, "sa");
m
L= diag(sum(A))-A;
[v, d] = eigs (L, 3, "sa");
d
gplot (A, v(:,[2 1]))

m shows that it is sorted from small to big
but d is wrong.

why the eigs of L are sorted different than the eigs of A is the problem.

(file #36443)

Doug Stewart <dastew>
Tue 23 Feb 2016 05:36:24 PM UTC, comment #4:

As a test example, I think this will do.

Upload the file tstdata.mat which will have the matrices before and after the call to eigs.

Rik <rik5>
Project Administrator
Tue 23 Feb 2016 04:34:12 PM UTC, comment #3:

Could the original reporter please post sample code that unambiguously reproduces the problem? There are a number of examples in the paper, but they seem to rely on data loaded from files which I don't have access to.

In creating an example, please do not upload any proprietary information, such as sample data from a Matlab installation.

Rik <rik5>
Project Administrator
Tue 23 Feb 2016 04:29:50 PM UTC, comment #2:

I made a sparse matrix and it also worked as expected.
I then went back to the original problem which involved a .mat file from a 3rd source.
When I load the data from the .mat file then the ordering of the eigs with 'sa' is wrong!

How or why the loading of the yaleShieldBig.mat file causes this to happen is not understood by me.

Doug Stewart <dastew>
Tue 23 Feb 2016 03:58:30 PM UTC, comment #1:

I tried a simple function and it works correctly

N=26
d = 10*rand(N,1);
q = triu(bsxfun(@min,d,d.').*rand(N),1);
A = diag(d)+q+q.';

[v d]=eigs(A,3,'sa')

So I think the problem might be with some sparse matrix

Doug Stewart <dastew>
Tue 23 Feb 2016 03:14:27 PM UTC, original submission:

When going through the examples contained in this document: http://www.cs.yale.edu/homes/spielman/561/lect01-15.pdf and in particular towards the end with the dodecahedron example, the output from Octave differs significantly from that of Matlab. In particular, the "sa" option for eigs() forces an ordering on the ouput eigenvalues/vectors which doesn't seem to be observed. When following through the example, the output is signficantly different.

I conferred with others on the IRC channel, and they confirmed that the bug is also in 4.0.

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #40438:  test_all_cases.m added by caliari (5KiB - text/x-objcsrc)
file #36496:  eigs_sa.patch added by rik5 (988B - text/x-diff)
file #36497:  tst_eigs.m added by rik5 (259B - d2l/unknowntype)
file #36443:  yaledata added by dastew (61KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by dastew (Posted a comment)
  • -unavailable- added by None (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 12 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 14 Nov 2017 05:46:38 AM UTCrik5StatusPatch Submitted=>Fixed
      Open/ClosedOpen=>Closed
    Tue 18 Apr 2017 03:10:42 PM UTCcaliariAttached File-=>Added test_all_cases.m, #40438
    Wed 29 Mar 2017 06:44:43 PM UTCcaliariStatusConfirmed=>Patch Submitted
    Sun 28 Feb 2016 01:39:25 AM UTCrik5Attached File-=>Added eigs_sa.patch, #36496
      Attached File-=>Added tst_eigs.m, #36497
    Sun 28 Feb 2016 01:24:44 AM UTCrik5Item GroupInaccurate Result=>Matlab Compatibility
    Tue 23 Feb 2016 08:32:13 PM UTCmtmillerSeverity3 - Normal=>2 - Minor
      StatusNone=>Confirmed
      Release3.8.2=>dev
      Summaryoutput of eigs() with \'sa\' option differs from Matlab=>output of eigs() with 'sa' option sorted differently from Matlab
    Tue 23 Feb 2016 06:20:59 PM UTCdastewAttached File-=>Added yaledata, #36443

    Back to the top


    Powered by Savane 3.1-cleanup1