bugGNU Octave - Bugs: bug #57185, eigs does not consistently order...

 
 

bug #57185: eigs does not consistently order eigenvalues

Submitter:  Muhali <muhali>
Submitted:  Thu 07 Nov 2019 11:19:07 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
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

Wed 04 Dec 2019 08:31:52 AM UTC, comment #17: 

As said, ARPACK does not guarantee any particular order. In eigs-base.cc there are some attempts to sort the eigenvalues. But they were workarounds in order to reproduce the old matlab behavior. Now, the eigenvalues in matlab are sorted (and Matlab does not use ARPACK anymore). I think we should remove from eigs-base.cc those old attempts and sort the results in eigs.m.

Marco Caliari <caliari>
Group Member
Fri 15 Nov 2019 02:04:39 AM UTC, comment #16: 

Seems like another addition for Matlab compatibility will be to map new options for the type of eigenvalues to return to their old two-letter codes based on symmetry and real/complex nature of input.


sigma                      Description                                    sigma (R2017a and earlier)
----------------------------------------------------------------------------------------------------

scalar (real or complex)   The eigenvalues closest to the number sigma.   No change

'largestabs' (default)     Largest magnitude.                             'lm'

'smallestabs'              Smallest magnitude. Same as sigma = 0.         'sm'

'largestreal'              Largest real.                                  'lr', 'la'

'smallestreal'             Smallest real.                                 'sr', 'sa'

'bothendsreal'             Both ends, with k/2 values with largest and    'be'
                           smallest real part respectively
                           (one more from high end if k is odd).



Rik <rik5>
Group administrator
Fri 15 Nov 2019 12:14:25 AM UTC, comment #15: 

I have found a test case where the first two values returned are not in the correct order.  The required variables are attached in the file 57185.var.  The m-file code is in tst_57185.m and is quite simple.


load ("57185.var")

[V,d] = eigs (tstmat, 5, 'lr', opts);
diag (d)


The result for me is


ans =

   0.7613
   0.6569
   1.0751
   1.4364
   2.1913


The order of the first two eigenvalues should be swapped.

One question is whether this is a valid way to call eigs.  I think it is.

Also, note that if only the eigenvalues are requested then Octave orders things correctly


d = eigs (tstmat, 5, 'lr', opts);


This might offer a clue to what is going on.

The second question is whether we should address this in the C++ code of _eigs_ or in eigs.m?

In the case above, simple code to sort and re-order the eigenvalues and eigenvectors is


[ds, idx] = sort (diag (d));
d = diag (ds);
V = V(:, idx);



Rik <rik5>
Group administrator
Fri 15 Nov 2019 12:03:26 AM UTC, comment #14: 

Adding David Bateman to CC list since he is deeply familiar with the eigs function and the interface to ARPACK.

Rik <rik5>
Group administrator
Thu 14 Nov 2019 10:38:05 PM UTC, comment #13: 

This is the behavior we need to match for Matlab compatibility.


Value of sigma   Output sorting
----------------------------------------------

'largestabs'     Descending order by magnitude

'largestreal'    Descending order by real part

'largestimag'    Descending order by imaginary part

'smallestabs'    Ascending order by magnitude

'smallestreal'   Ascending order by real part
'bothendsreal'   Ascending order by real part

'smallestimag'   Ascending order by imaginary part

'bothendsimag'   Descending order by absolute value of imaginary part


Rik <rik5>
Group administrator
Sat 09 Nov 2019 07:39:44 PM UTC, comment #12: 

The issue of the ordering of singular values depends on the underlying function eigs().  That function does not guarantee a particular ordering.  However, Matlab has changed behavior relatively recently (R2017b).  They now state


Behavior changed in R2017b

    Changes to sorting order of output

    eigs now sorts the output according to the value of sigma. For example, the command eigs(A,k,'largestabs') produces k eigenvalues sorted in descending order by magnitude.

    Previously, the sorting order of the output produced by eigs was not guaranteed.

    Reproducibility

    Calling eigs multiple times in succession now produces the same result. Set 'StartVector' to a random vector to change this behavior.

    Display

    A display value of 2 no longer returns timing information. Instead, eigs treats a value of 2 the same as a value of 1. Also, the messages shown by the 'Display' option have changed. The new messages show the residual in each iteration, instead of the Ritz values.


Re-titling bug report to reflect underlying issue.

Rik <rik5>
Group administrator
Sat 09 Nov 2019 07:35:27 PM UTC, comment #11: 

I checked in a fix for svds returning complex results here: https://hg.savannah.gnu.org/hgweb/octave/rev/2213e82bb9d3.

Rik <rik5>
Group administrator
Sat 09 Nov 2019 03:28:58 PM UTC, comment #10: 

I filed a bug report about eigs() not returning strictly real values for symmetric inputs (bug #57196).  That would fix this, but in the meantime, it is appropriate to apply the patch from comment #9.

Rik <rik5>
Group administrator
Fri 08 Nov 2019 09:41:47 PM UTC, comment #9: 

Attached is a patch for the first issue of complex values.

Maybe we should be detecting whether the input to eigs() is symmetric and, if so, guarantee real results from that function.  But this patch was certainly easy to implement and does the right thing.

(file #47825)

Rik <rik5>
Group administrator
Fri 08 Nov 2019 09:25:57 PM UTC, comment #8: 

For the second question, I think we need to check the code.  In svds.m there is


    if (nargout > 1)
      [V, s, flag] = eigs ([sparse(m,m), b; b', sparse(n,n)],
                           b_k, b_sigma, b_opts);
      s = diag (s);
    else
      s = eigs ([sparse(m,m), b; b', sparse(n,n)], b_k, b_sigma, b_opts);
    endif


so we are calling eigs with a different function prototype in each case.  I seem to recall that if only the singular values are requested we may perform a sort on them to ensure uniform ordering.  But, we can't do that arbitrarily when someone has requested the full decomposition because it would also require changing the U and V matrices to ensure that the decomposition


A = U * S * V'


still holds. 

Rik <rik5>
Group administrator
Fri 08 Nov 2019 09:08:42 PM UTC, comment #7: 

I think we should consider fixing this is two steps.  The first question to resolve is the return of complex singular values rather than just real ones.  The complex values are very small, which suggests the algorithm is effective but there is a small loss of precision owing to the fact that this is imperfct numerical computation.  Sample:


octave:9>  x = randn(10, 5) ;
octave:10> y = randn(10, 5) ;
octave:11> z = complex(x, y) ;
octave:12> s = svds (z)
s =

   5.9713 + 0.0000i
   4.6205 - 0.0000i
   3.7797 + 0.0000i
   2.5995 + 0.0000i
   1.9744 - 0.0000i

octave:13> imag (s)
ans =

   9.0733e-17
  -1.2895e-16
   9.0859e-17
   2.8342e-19
  -7.4995e-17


A simple answer might be to take the equivalent, in C++, of


s = real (s);


before returning the values.

Alternatively, it make be possible to declare s as a real array in C++ so that it can never hold imaginary values.

Rik <rik5>
Group administrator
Fri 08 Nov 2019 12:10:33 PM UTC, comment #6: 

Now this makes sense, same for default hg_id 959a63d1ede6:


x = randn(10, 5) ;
y = randn(10, 5) ;
z = complex(x, y) ;

[~, Dr] = svds(x);
svds_x = [diag(Dr), svds(x), svd(x)]
[~, Dc] = svds(z);
svd_z = [diag(Dc), svds(z), svd(z)]

isreal(Dr)
arrayfun (@isreal, svd_z)


The output


svds_x =

   5.1420   5.1420   5.1420
   3.0993   3.0993   3.0993
   2.7733   2.7733   2.7733
   2.1100   2.1100   2.1100
   1.4318   1.4318   1.4318

svd_z =

   2.6499 + 0.0000i   7.3784 - 0.0000i   7.3784 + 0.0000i
   3.8673 + 0.0000i   5.7718 - 0.0000i   5.7718 + 0.0000i
   4.3697 - 0.0000i   4.3697 - 0.0000i   4.3697 + 0.0000i
   5.7718 - 0.0000i   3.8673 + 0.0000i   3.8673 + 0.0000i
   7.3784 + 0.0000i   2.6499 - 0.0000i   2.6499 + 0.0000i

ans = 1
ans =

  0  0  1
  0  0  1
  0  0  1
  0  0  1
  0  0  1


Thus two problems:

1. The order of the singular values depends on the number of output arguments

2. The singular values are complex values.

Kai Torben Ohlhus <siko1056>
Group Member
Fri 08 Nov 2019 09:35:02 AM UTC, comment #5: 

sorry, the function in question was 'svds' not 'svd'.


x = randn(10, 5) ;
y = randn(10, 5) ;
z = complex(x, y) ;

[~, Dr] = svds(x) ; diag(Dr)
svds(x)
[~, Dc] = svds(z) ; diag(Dc)
svds(z)

isreal(Dr)
isreal(Dc)


Muhali <muhali>
Fri 08 Nov 2019 04:08:04 AM UTC, comment #4: 

Agree to comment #1 and comment #2.  Nothing suspicious using


>> __blas_version__ ()
ans = OpenBLAS (config: OpenBLAS 0.3.6 NO_LAPACKE DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=48)
>> __lapack_version__ ()
ans = Linear Algebra PACKage Version 3.8.0


on openSUSE 15.1 using the flatpak installation or a custom build.

Same suspicion as comment #3: What is your output for the functions above and did you customize your Octave somehow (please explain how you installed Octave)?

Kai Torben Ohlhus <siko1056>
Group Member
Fri 08 Nov 2019 03:48:58 AM UTC, comment #3: 

A broken blas/lapack is a suspect...

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 08 Nov 2019 02:26:10 AM UTC, comment #2: 

Also works for me.  In case it was something related to random numbers, I tried several thousand runs with


for i = 1:1000
  x = randn(10, 5) ;
  y = randn(10, 5) ;
  z = complex(x, y) ;

  [~, Dr] = svd(x) ; diag(Dr);
  svd(x);
  [~, Dc] = svd(z) ; diag(Dc);
  svd(z);

  assert (isreal(Dr))
  assert (isreal(Dc))
endfor


It always passes.

Rik <rik5>
Group administrator
Thu 07 Nov 2019 04:46:11 PM UTC, comment #1: 

I get consistent results with 5.1 and 6.0, so I'm not seeing what you are. Are you saying that Dc is not real on your system?

Example output for me (always helps to show your results):


>> [~, Dc] = svd(z) ; diag(Dc)
ans =

   7.0210
   5.2298
   4.5198
   3.4535
   2.0099

>> svd(z)
ans =

   7.0210
   5.2298
   4.5198
   3.4535
   2.0099

>> isreal(Dc)
ans = 1


Mike Miller <mtmiller>
Group Member
Thu 07 Nov 2019 11:19:07 AM UTC, original submission:  


x = randn(10, 5) ;
y = randn(10, 5) ;
z = complex(x, y) ;

[~, Dr] = svd(x) ; diag(Dr)
svd(x)
[~, Dc] = svd(z) ; diag(Dc)
svd(z)

isreal(Dr)
isreal(Dc)


This shows that in the complex case, the ordering depends on having explicit output (Dz) or not. This is inconsistent and also different from Matlab.

svd should also return real and not complex eigenvalues.

Muhali <muhali>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #47871:  57185.var added by rik5 (2KiB - application/octet-stream)
file #47872:  tst_57185.m added by rik5 (68B - text/x-matlab)
file #47825:  bug57185.real.cset added by rik5 (981B - application/octet-stream)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by muhali (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 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-11-15 rik5 Attached File- Added 57185.var, #47871
        Attached File- Added tst_57185.m, #47872
    2019-11-15 rik5 Carbon-Copy- Added dbateman
    2019-11-09 rik5 Severity4 - Important 3 - Normal
        Item GroupInaccurate Result Matlab Compatibility
        Summarysvds returns inconsistent singular values eigs does not consistently order eigenvalues
    2019-11-09 rik5 Dependencies- Depends on bugs #57196
    2019-11-09 rik5 Release5.1.0 dev
    2019-11-08 rik5 Attached File- Added bug57185.real.cset, #47825
    2019-11-08 siko1056 Severity3 - Normal 4 - Important
        StatusWorks For Me Confirmed
        Summarysvd returns inconsistent eigenvalues svds returns inconsistent singular values
    2019-11-08 siko1056 CategoryNone Octave Function
        Item GroupNone Inaccurate Result
    2019-11-07 mtmiller StatusNone Works For Me

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code