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.