eigs(speye(21)) does not work

[ViralBShah]

eigs does not work on identity matrices of size 21 and higher for the default 6 eigenvalues.

julia> eigs(speye(21), nev=6)
ERROR: ARPACKException(3)
 in aupd_wrapper at linalg/arpack.jl:47
 in eigs at linalg/arnoldi.jl:35

julia> eigs(speye(21), nev=9)
ERROR: ARPACKException(3)
 in aupd_wrapper at linalg/arpack.jl:47
 in eigs at linalg/arnoldi.jl:35

julia> eigs(speye(21), nev=10)
([1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0],
21x10 Array{Float64,2}:
  0.252111    -0.246736     0.0192645    0.234598    …   0.130942     0.100895   -0.0290075 
  0.118993     0.190852    -0.393942    -0.232889       -0.002121     0.0280335  -0.213672  
 -0.00306955  -0.157261    -0.0341907   -0.182208       -0.110704    -0.158746   -0.419363  
  0.0751385    0.167007     0.0305194    0.0977773      -0.00449555  -0.318633   -0.0253858 
 -0.0943686   -0.0731143    0.148546     0.330696        0.0214212   -0.213419   -0.0260602 
  0.140182     0.00836135   0.00474463  -0.387105    …  -0.0425574   -0.111232   -0.135756  
 -0.256888    -0.354937    -0.482451     0.0503225       0.626875    -0.119688   -0.0743748 
  0.115359     0.203424    -0.0082043    0.00844204      0.0697025   -0.31223    -0.384733  
 -0.0716209   -0.346583     0.529378    -0.307466        0.325135     0.238119   -0.100524  
  0.132531    -0.230142    -0.0277612   -0.0306749      -0.174249     0.271229    0.244908  
  0.149707     0.438754     0.286695     0.125643    …   0.42474     -0.149599   -0.0382037 
 -0.0220212    0.102343     0.155978    -0.232871       -0.186025     0.0696285  -0.211983  
 -0.163692    -0.167343    -0.0155797    0.171629       -0.261443     0.203264   -0.480922  
 -0.0956       0.0534561   -0.296926    -0.143689       -0.116822     0.0899562   0.00621333
  0.63577     -0.422864     0.0268107   -0.0998203      -0.0911458   -0.403187    0.106614  
  0.41942      0.15776     -0.034414     0.363393    …   0.0734304    0.475078   -0.205165  
  0.11205     -0.0650867   -0.113777     0.0755098      -0.0801624   -0.0160499  -0.143335  
 -0.0523192    0.0961701    0.0374504   -0.0913321       0.0863463    0.0442572   0.211717  
 -0.242706    -0.0330528    0.0144394    0.226972       -0.31918     -0.267941    0.241889  
  0.0206762   -0.18865     -0.0265206    0.393166       -0.061266    -0.0446505  -0.141417  
 -0.269261    -0.128465     0.307113     0.0916609   …  -0.0646332   -0.146322   -0.243291  )

[TXKurt]

Here is the ARPACK comment on return value 3:

c = 3: No shifts could be applied during a cycle of the
c Implicitly restarted Arnoldi iteration. One possibility
c is to increase the size of NCV relative to NEV.
c See remark 4 below.

Here is "remark 4"

c 4. At present there is no a-priori analysis to guide the selection
c of NCV relative to NEV. The only formal requrement is that NCV > NEV.
c However, it is recommended that NCV .ge. 2_NEV. If many problems of
c the same type are to be solved, one should experiment with increasing
c NCV while keeping NEV fixed for a given test problem. This will
c usually decrease the required number of OP_x operations but it
c also increases the work and storage required to maintain the orthogonal
c basis vectors. The optimal "cross-over" with respect to CPU time
c is problem dependent and must be determined empirically.
c

Presently, ncv is set in the first line of aupd_wrapper() in arpack.jl:

ncv = min(max(nev*2, 20), n)

Maybe the caller should have access to setting the parameter ncv...

[ViralBShah]

Yes I thought about the caller having access to ncv, but this really should work out of the box.

[JeffBezanson]

Do we know how to fix this?

[ViralBShah]

With commit 96dcf38 this now works when ncv == n. Clearly not a solution, since ncv x n is the size of the workspace allocated.

[ViralBShah]

@alanedelman Any insights here? I wonder if there is an error with our ARPACK interface or if this has something to do with the Arnoldi factorization.

[TXKurt]

This command works in octave and they are also using ARPACK. Maybe look though the octave interface?

http://octave.sourceforge.net/doxygen/html/d3/d5d/eigs-base_8cc_source.html

At first glance, it looks like the ncv is chosen in the same way:

00815 if (p < 0)
00816 {
00817 p = k * 2;
00818
00819 if (p < 20)
00820 p = 20;
00821
00822 if (p > n - 1)
00823 p = n - 1 ;
00824 }

Sorry I don't have time to look at this any further..
Kurt

[ViralBShah]

I did take a quick look at octave, but could not find the relevant code. Thanks for pointing to it. I wonder what we are doing differently in the problem setup.

[TXKurt]

Look at the link to the octave source above, lines 938-1002
They are only quitting when info < 0, so maybe the info==3 that we are getting is just some information, and we can continue... ?

[ViralBShah]

That does not seem to be the right thing though.

c          Error flag on output.
c          =  0: Normal exit.
c          =  1: Maximum number of iterations taken.
c                All possible eigenvalues of OP has been found. IPARAM(5)
c                returns the number of wanted converged Ritz values.
c          =  2: No longer an informational error. Deprecated starting
c                with release 2 of ARPACK.
c          =  3: No shifts could be applied during a cycle of the
c                Implicitly restarted Arnoldi iteration. One possibility
c                is to increase the size of NCV relative to NEV.
c                See remark 4 below.

[TXKurt]

You're right. I also just tried it and didn't get the right answer (last eigenvalue/vector was wrong)…
-Kurt

Am 20.09.2013 um 06:07 schrieb Viral B. Shah:

That does not seem to be the right thing though.

c Error flag on output.
c = 0: Normal exit.
c = 1: Maximum number of iterations taken.
c All possible eigenvalues of OP has been found. IPARAM(5)
c returns the number of wanted converged Ritz values.
c = 2: No longer an informational error. Deprecated starting
c with release 2 of ARPACK.
c = 3: No shifts could be applied during a cycle of the
c Implicitly restarted Arnoldi iteration. One possibility
c is to increase the size of NCV relative to NEV.
c See remark 4 below.
—
Reply to this email directly or view it on GitHub.

[vtjnash]

tolerance is an output parameter from the first iteration, which is needed in all future uses of the data in arpack. we were resetting it to 0 (if it couldn't find a solution in the first pass), which would certainly be tough to meet at any size.

[ViralBShah]

Thanks. This is awesome detective work.