Arpack-ng bug preventing us upgrading from 3.5 to 3.8

[ViralBShah]

If I have done my analysis right - the failure is introduced in https://github.com/opencollab/arpack-ng/compare/0e7d01d..7fc42e5

I'm not sure if we are doing something wrong in our wrappers or Arpack-ng introduced a bug, but it is somewhere in there. The relevant issue is #118.

[ViralBShah]

@andreasnoack does eyeballing that diff suggest something insightful?

[andreasnoack]

My guess would be opencollab/arpack-ng@d04bdf1 but I'm not really sure.

[ViralBShah]

That's my best bet too. I wonder if we are calling arpack differently than other systems that makes this buggy for us.

[fghoussen]

@ViralBShah: did you try to remove the patch opencollab/arpack-ng#80 (comment)? Does it solve the problem?
@caliarim: could you have a look?

[caliarim]

Both the problem described here and the problem related to "Force the residual vector to be in the range of OP" depend on a bad choice of the initial random vector. Many iterative algorithms have some conditions on the initial values which we usually cannot verify and we hope for the best. Matlab has a similar problem with normest(toeplitz([-2,1,0,0])). So, I have no solution.
Of course, this is a general comment, it could be that for this specific example there is a specific bug, but I have no time to investigate now.

[fghoussen]

depend on a bad choice of the initial random vector

@ViralBShah: how to do initialize your vectors?

[ViralBShah]

@fghoussen We do not initialize the vectors: https://github.com/JuliaLinearAlgebra/Arpack.jl/blob/master/src/libarpack.jl#L109

What should we be doing here? I thought ARPACK would initialize them for us.

[fghoussen]

What should we be doing here?

@caliarim may answer. I never had specific problems with initialization (but a lot on computations! :))

[ViralBShah]

I'll note that in #118, just using a rand initialized v0 is sufficient to not trigger the error any more in the case described.

cc @amontoison as well.

[fghoussen]

just using a rand initialized v0 is sufficient to not trigger the error any more

Would say make sense to me (in difficult cases). AFAIR, this way, one may better span the vectorial space solution.
Hope @caliarim could give better explanations / clarifications. Try to init with rand: does it work better?

[ViralBShah]

Yes, init with rand does work better.

[amontoison]

@ViralBShah
It's relevant to have a random v0.
It should be related to the initial vector of the Arnoldi process.
We use it to generate the Krylov subspace and the associated projection (Upper Hessenberg matrix).
It's the b vector here.
The eigenvalues of the Upper Hessenberg (Ritz values) converge to the dominant eigenvalues of A because for a random vector v0, A^n * v0 converges to the eigenvector associated to the dominant eigenvalue of A.
We expect to have this behaviour when v0 doesn't have a specific structure.

[fghoussen]

BTW, in corner cases where arpack would still fail, you can restart arpack using the previous "not-good-enough" solution as v0: for difficult cases, you may run arpack twice, but, the second time using a better initial guess. Hope this helps!

[ViralBShah]

@amontoison I was of the opinion that if we don't provide one, ARPACK will automatically pick a starting initial vector. Would it make sense to use rand to initialize the v0 by default in Arpack.jl?

@fghoussen Thanks for the suggestion. We can add that to the documentation.

[ViralBShah]

Ok, so with ARPACK-ng 3.8, the problem that is provided in #118 fails reliably on every test, even when I give it a starting v0. With ARPACK 3.5, it seems to work reliably with a random v0.

I have created #145 which makes it easy to try out and debug ARPACK 3.8.

[fghoussen]

Could you export your A/B matrix to ASCII matrix market format?
If so, we'll have more chance to understand if the problem comes from arpack-ng or not.

[ViralBShah]

https://gist.github.com/ViralBShah/ae892fe6df0884185313f3f0e10c03cc

[fghoussen]

OK, I'll try to have a look when possible.
What are the constraints of your solve: complex numbers in double precision? is A sym? nb EV? nb CV? Magnitude? Want / need to shift / invert? Tol needed? Ritz or Schur vectors? File is 1-based? etc

[fghoussen]

Seems I can compute OK the 200 first EV with LM magnitude, no shift, no invert, regular tol / nb it.
Ritz vectors OK, but Schur vector fail.

[ViralBShah]

The actual failure case is when doing d, rest = eigs(L, nev=2, which=:LR), in Julia's Arpack.jl (actual Julia code in #118). I know there will be questions about how things are being set up etc., but the behaviour changed between 3.5 and 3.8.

The file is what MatrixMarket.jl generates, I assume 1-based.

[fghoussen]

OK, I'll have a look when possible.
Is the problem symmetric? Can you regenerate the file formatting lines this way i j (real_ij, imag_ij)?

[ViralBShah]

It is not.

[fghoussen]

Can you regenerate the file formatting lines this way i j (real_ij, imag_ij) (unless I can't use the file)?
Missing (, , and ) prevent me from reading the file correctly