Search code examples
octavenumerical-methodseigenvalue

Octave eigs() function bugged?


Running Octave 6.3.0 for Windows. I need to get the smallest eigenvalue of some matrix.eigs(A,1,"sm") is supposed to do that, but I often get wrong results with singular matrices.

eigs(A) (which returns all the the first 6 eigenvalues/vectors) is correct (at least at the machine precision):

>> A = [[1 1 1];[1 1 1];[1 1 1]]
A =

   1   1   1
   1   1   1
   1   1   1

>> [v lambda flag] = eigs(A)
v =

   0.5774  -0.3094  -0.7556
   0.5774  -0.4996   0.6458
   0.5774   0.8091   0.1098

lambda =

Diagonal Matrix

   3.0000e+00            0            0
            0  -4.5198e-16            0
            0            0  -1.5831e-17

flag = 0

But eigs(A,1,"sm") is not:

>> [v lambda flag] = eigs(A,1,"sm")
warning: eigs: 'A - sigma*B' is singular, indicating sigma is exactly an eigenvalue so convergence is not guaranteed
warning: called from
    eigs at line 298 column 20

warning: matrix singular to machine precision
warning: called from
    eigs at line 298 column 20

warning: matrix singular to machine precision
warning: called from
    eigs at line 298 column 20

warning: matrix singular to machine precision
warning: called from
    eigs at line 298 column 20

warning: matrix singular to machine precision
warning: called from
    eigs at line 298 column 20

v =

  -0.7554
   0.2745
   0.5950

lambda = 0.4322
flag = 0

Not only the returned eigenvalue is wrong, but the returned flag is zero, indicating that every went right in the function...

Is it a wrong usage of eigs() (but from the doc I can't see what is wrong) or a bug?

EDIT: if not a bug, at least a design issue... No problem either when requesting the 2 smallest values instead of the smallest value alone.

>> eigs(A,2,"sm")
ans =

  -1.7700e-17
  -5.8485e-16

EDIT 2: the eigs() function in Matlab online just runs fine and return the correct results (at the machine precision)

>> A=ones(3)

A =

     1     1     1
     1     1     1
     1     1     1

>> [v lambda flag] = eigs(A,1,"smallestabs")

v =

   -0.7556
    0.6458
    0.1098


lambda =

  -1.5831e-17


flag =

     0

Solution

  • After more tests and investigations I think I can answer that yes, Octave eigs() has some flaw.

    eigs(A,1,"sm") likely uses the inverse power iteration method, that is repeatedly solving y=A\x, then x=y, starting with an arbitrary x vector. Obviously there's a problem if A is singular. However:

    • Matlab eigs() runs fine in such case, and returns the correct eigenvalue (at the machine precision). I don't know what it does, maybe adding a tiny value on the diagonal if the matrix is detected as singular, but it does something better (or at least different) than Octave.
    • If for some (good or bad) reason Octave's algorithm cannot handle a singular matrix, then this should be reflected in the 3rd return argument ("flag"). Instead, it is always zero as if everything went OK.

    eigs(A,1,"sm") is actually equivalent to eigs(A,1,0), and the more general syntax is eigs(A,1,sigma), which means "find the closest eigenvalue to sigma, and the associated eigenvector". For this, the inverse power iteration method is applied with the matrix A-sigma*I. Problem: if sigma is already an exact eigenvalue this matrix is singular by definition. Octave eigs() fails in this case, while Matlab eigs() succeeds. It's kind of weird to have a failure when one knows in advance the exact eigenvalue, or sets it by chance. So the right thing to do in Octave is to test if (A-sigma.I) is singular, and if yes add a tiny value to sigma: eigs(A,1,sigma+eps*norm(A)). Matlab eigs() probably does something like that.