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

 
 

bug #57185: eigs does not consistently order eigenvalues

Submitted by:  Muhali <muhali>
Submitted on:  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

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)

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

    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 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.5