sparse matrix eigs, svds fails

[alanedelman]

eigs and svds run into trouble computing too many eigenvalues/singular values
not sure if this is expected behavior

julia> load("sparse.jl")

julia> load("suitesparse.jl")

julia> load("arpack.jl")

julia> A=speye(4)
4x4 sparse matrix with 4 nonzeros:
[1, 1] = 1.0
[2, 2] = 1.0
[3, 3] = 1.0
[4, 4] = 1.0

julia> svds(A)
error code -3 from ARPACK saupd
in svds at /home/jeffb/julia/extras/arpack.jl:230
in svds at /home/jeffb/julia/extras/arpack.jl:259

julia> eigs(A)
Compute fewer eigenvalues using eigs(A, k)
in eigs at /home/jeffb/julia/extras/arpack.jl:83
in eigs at /home/jeffb/julia/extras/arpack.jl:187

julia> eigs(A,2)
([1.0, 1.0],
4x2 Float64 Array:
-0.0760132 -0.534781
-0.521881 0.360963
0.539126 0.685331
-0.656662 0.337693)

[JeffBezanson]

Likely same issue as #918.

[andreasnoack]

Another thing. The eigenvectors seem to be random. Also, there are some problems with symmetry.

julia> srand(126)

julia> Asp=sprandn(10,10,0.2);AspSym=Asp'Asp;

julia> full(AspSym-AspSym')
10x10 Float64 Array:
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0      
 0.0  0.0      0.0   0.0       0.00972277  0.13768   0.0  0.0  0.0  -2.07257  
 0.0  0.0      0.0  -0.262696  0.0         0.0       0.0  0.0  0.0   0.0      
 0.0  0.0      0.0   0.0       0.0         0.154933  0.0  0.0  0.0   0.0      
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0108206
 0.0  0.13768  0.0   0.0       0.0284376   0.0       0.0  0.0  0.0  -0.0805007
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0      
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0      
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0473264
 0.0  0.0      0.0   0.0       0.0         0.0       0.0  0.0  0.0   0.0      

julia> full(AspSym)-full(AspSym)'
10x10 Float64 Array:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> eigs(speye(5), 3)
([1.0, 1.0, 1.0],
5x3 Float64 Array:
  0.239634   -0.569204  -0.602892
 -0.64662    -0.27262    0.161997
  0.32331    -0.69591    0.336256
 -0.0425561  -0.207542   0.655278
  0.64662     0.27262    0.260424)

julia> eigs(speye(5), 3)
([1.0, 1.0, 1.0],
5x3 Float64 Array:
  0.243169  -0.839674  -0.336662
 -0.327944   0.279359  -0.389231
  0.260929   0.194977  -0.807363
  0.810981   0.140716   0.256699
 -0.327944  -0.39886    0.132015)

[ViralBShah]

Does it work on matlab?

[ViralBShah]

Oh it is for small matrices. Perhaps just fall to eig/svd for less than 20.

[ViralBShah]

@andreasnoackjensen I do not think I have implemented the symmetric drivers from ARPACK. Should be easy to do, and once the SymmetricMatrix stuff is sorted out, it should be convenient to use as well.

[dmbates]

@JeffBezanson This isn't the same issue as #918. The obscure error code means that the number of columns in A is less than the number of singular values requested, which defaults to 6. I'll commit a fix once I have tested it.

[andreasnoack]

@ViralBShah I think the symmetric driver is implemented. The function tests for symmetry and calls saupd if it is. But the problem is that I there is a bug in at least ref and -.

julia> srand(127)

julia> A=sprandn(5,5,0.2);As=A'A
5x5 sparse matrix with 8 nonzeros:
    [1, 1]  =  0.593517
    [4, 1]  =  -0.404331
    [2, 2]  =  1.73909
    [3, 2]  =  0.611348
    [2, 3]  =  0.611348
    [3, 3]  =  0.214909
    [4, 4]  =  4.20678
    [1, 4]  =  -0.404331

looks symmetric, but

julia> issym(As)
false

and it is not the floating point problem

julia> As-As'
5x5 sparse matrix with 2 nonzeros:
    [1, 4]  =  0.404331
    [1, 4]  =  -0.404331

Hm. Two (1,4) elements. Finally, even though show told me that As[1,4]=-0.404331 I get

julia> As[1,4]
0.0

[dmbates]

Of course it gets a little more complicated. ``nev`, the number of eigenvalues or singular values requested must be less than the size of a square matrix or the number of columns in a rectangular matrix. So for speye(4) the best you can do is to get 3 singular values.

julia> A = speye(4)
4x4 sparse matrix with 4 nonzeros:
    [1, 1]  =  1.0
    [2, 2]  =  1.0
    [3, 3]  =  1.0
    [4, 4]  =  1.0


julia> svds(A, 3)
(
4x3 Float64 Array:
 0.485473   0.700208   0.0775653
 0.417536   0.296277  -0.464294 
 0.644222  -0.375598   0.638083 
 0.418279  -0.529957  -0.609315 ,

[1.0, 1.0, 1.0],
3x4 Float64 Array:
 0.485473    0.417536   0.644222   0.418279
 0.700208    0.296277  -0.375598  -0.529957
 0.0775653  -0.464294   0.638083  -0.609315)

With eigs one of the checks is

           ncv = min(max(nev*2, 20), n)
           if ncv-nev < 2 || ncv > n error("Compute fewer eigenvalues using eigs(A, k)") end

which I think is unmodified from @viralshah 's original. The second condition in the if statement doesn't make sense given the previous line (because of the min, you must have ncv <= n) and I'm not sure where the first condition comes from. As I read the documentation you must have ncv-nev > 1, not 2.

Anyway, if I comment out that if statement then it still works as I have added a condition that nev <= n - 1.

I'll commit and let me know if I have botched things up.

[ViralBShah]

The case from #918 still fails.

julia> eigs(rand(4,4))
error code -3 from ARPACK aupd
 in eigs at /Users/viral/julia/extras/arpack.jl:84
 in eigs at /Users/viral/julia/extras/arpack.jl:189

From dnaupd.f:
-3: NCV-NEV >= 2 and less than or equal to N.

@andreasnoackjensen I guess it was a long time back I implemented this - and haven't looked at it recently. Seems like the issues you have discovered are sparse related, and should be a separate issue.

[ViralBShah]

These are the constraints on NEV and NCV. As a result, NEV should be at most N-2, and NCV should be at most N. If these constraints are not satisfied, I would much rather return an error, and suggest that the user use eig or svd.

Also, how many Arnoldi vectors should we use in general? Currently, the code tries to use 2*NEV Arnoldi vectors, if the problem permits. What should it really do?

 NEV     Integer.  (INPUT)
         Number of eigenvalues of OP to be computed. 0 < NEV < N-1.

NCV     Integer.  (INPUT)
         Number of columns of the matrix V. NCV must satisfy the two
         inequalities 2 <= NCV-NEV and NCV <= N.
         This will indicate how many Arnoldi vectors are generated
         at each iteration.  After the startup phase in which NEV
         Arnoldi vectors are generated, the algorithm generates
         approximately NCV-NEV Arnoldi vectors at each subsequent update
         iteration. Most of the cost in generating each Arnoldi vector is
         in the matrix-vector operation OP*x.

[StefanKarpinski]

On Thu, Nov 22, 2012 at 4:04 AM, Viral B. Shah [email protected]:

These are the constraints on NEV and NCV. As a result, NEV should be at
most N-2, and NCV should be at most N. If these constraints are not
satisfied, I would much rather return an error, and suggest that the user
use eig or svd.

While that's fine for interactive usage, it would be deeply annoying for
something that is deployed to work until some unusual situation happens and
then fail when the solver could have produced a correct result.

[ViralBShah]

This is not just about getting all the eigenvalues for small 4x4 matrices. It also applies to getting all the eigenvalues or singular values for an nxn matrix. Just falling back on eig/svd leads to the use of a different algorithm in a non-transparent way, with performance characteristics that are hard to predict.

[JeffBezanson]

Is it really that bad?

If somebody writes eigs(A,n) where n might be the size of A in some cases, the burden of that performance is on that programmer. Either way, they will have to end up writing

if size_condition(n)
  x = eigs(A,n)
else
  x = eig(A)
  ...
end

either for performance or to avoid the exception.

It might make sense to give an error instead of horribly bad performance that will be hard to track down. But in that case, it should only kick in for large sizes where the performance is a real problem, and work for 4x4.

[ViralBShah]

The issue is not about 4x4, but about providing more than n-2 eigenvalues and eigenvectors - if I have understood the ARPACK documentation correctly. Is it so bad that eigs and svds give only n-2 eigenvalues and eigenvectors?

[ViralBShah]

BTW, the new eig that we are going to have is going to support fewer than n eigenvalues, if I have understood @andreasnoackjensen correctly.

[andreasnoack]

@ViralBShah Correct, but only for symmetric matrices. And the new solver is fast (but is not merged in this form yet so you cannot run the following).

julia> load("arpack")

julia> A=sprandn(1000,1000,0.1);Asym=A'A;

julia> min([@elapsed eigs(Asym) for i = 1:5])
1.3453099727630615

julia> min([@elapsed eigs(full(Asym)) for i = 1:5])
0.34864091873168945

julia> min([@elapsed (Z=Array(Float64, 1000, 6); LAPACK.syevr!('V', 'I', 'U', full(Asym), 0.0, 0.0, 995, 1000, Z, -1.0)) for i = 1:5])
0.2948329448699951

julia> Z=Array(Float64, 1000, 6); vals = LAPACK.syevr!('V', 'I', 'U', full(Asym), 0.0, 0.0, 995, 1000, Z, -1.0);

julia> norm(abs(eigs(Asym)[2]) - abs(Z))
8.550264651753738e-13

julia> norm(eigs(Asym)[1] - vals)
1.5701547587165274e-12

I vote for a unified (in eig and svd) approach with dispatch for returning some eigenvalues/vectors and then get rid of eigs and svds. I guess that the most users don't care which PACK that does the work behind the scenes.

[ViralBShah]

You still need eigs and svds for very large sparse matrices - so can't just get rid of it.

[andreasnoack]

Why can't eig have a method for sparse matrices that does what eigs is now doing?

[ViralBShah]

Without eigs, how would the user decide if he wants to use the Arnoldi methods or the direct factorization? What about all the other options that eigs supports?

http://www.mathworks.in/help/matlab/ref/eigs.html

We would then end up where we are in \, where so much functionality is crammed into it, that it is impossible to tell what is happening. We would have to then have even more options supported.

[andreasnoack]

My idea was to simplify the eigenvalue front end by removing and the s in eigs and overload eig with most of the methods of eigs. A few methods with few arguments would need to go, but no functionality. For dense matrices, if you only want some vectors the Arnoldi methods are used or maybe the new LAPACK subroutine for symmetric matrices and when all are requested LAPACK is called. For sparse matrices, the Arnoldi methods are used with a convention of only providing some of the vectors. Like is now for eigs.

The point is that when using eig you are interested in computing some eigenvalues and vectors but not how they are computed. You just expect them to be computed as efficient as possible and returned in the most obvious way, but of course it should be possible to figure out the algorithm from the documentation. The LAPACK routines are accessible either directly or through the factorisations when you need more control. There could be similar wrappers around ARPACK, but I don't think that eigs is the most obvious name for that even though it is MATLAB's choice.

Maybe front ends for different PACKs should be given different names even though the mathematical problem they solve is the same, but for many problems there have been chosen otherwise. * sometimes calls BLAS but not always. \ sometimes calls LAPACK and sometimes something in SuiteSparse and so forth.

[StefanKarpinski]

+1 for Andreas' proposal. I actually think that would be welcomed by many Matlab users even though there would initially be a little bit of friction. "Wait, what? I don't need to think about sparse vs. dense in this code? Ok."

[ViralBShah]

It is not that eig is for dense and eigs is for sparse in matlab. Both work on dense and sparse. The thing about eigs is that it uses ARNOLDI methods to compute a few eigenvalues and eigenvectors, which may be the largest, smallest, or in the neighbourhood of a particular eigenvalue, etc. This makes the eigs interface do a lot of things.

Hence, the question really is if all of this can be stuck into eig and svd without getting confusing. If we can do it in a sane way, with underlying equivalents of what we have done with Factorizations and stuff, it would certainly be worthwhile.

[alanedelman]

Mathematica allows for a Method option which defaults to AUTOMATIC, but can also be
ARNOLDI and LAPACK.

If you read the online documentation it is written in such a manner as to protect the user from too much
arcane knowledge.
(This may be a mostly correct policy, as long as you can drill down if you want to)

Automatic goes like this:
For eigenvalue computation when the input is an matrix of machine numbers and the number of eigenvalues requested is less than 20% of an Arnoldi method is used. Otherwise, if the input matrix is numerical then a method using LAPACK routines is used.

Any of the commands Eigenvalues, Eigenvectors, and Eigensystem can have an optional
argument such as Method->AUTOMATIC, Method->Arnoldi, Method->LAPACK
They seem to work for the vanilla and the generalized eigen problems.

See http://reference.wolfram.com/legacy/v5_2/Built-inFunctions/AdvancedDocumentation/LinearAlgebra/LinearAlgebraInMathematica/MatrixComputations/EigensystemComputations/AdvancedDocumentationLinearAlgebra3.3.4.html

[ViralBShah]

We will only be able to do something about this after keyword args. Even eigs as it stands is a real pain without keyword args.

[ViralBShah]

9aff106 adds support for keyword arguments for eigs and svds.