Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Arpack-ng bug preventing us upgrading from 3.5 to 3.8 #138

Open
ViralBShah opened this issue Nov 17, 2021 · 23 comments
Open

Arpack-ng bug preventing us upgrading from 3.5 to 3.8 #138

ViralBShah opened this issue Nov 17, 2021 · 23 comments

Comments

@ViralBShah
Copy link
Collaborator

ViralBShah commented Nov 17, 2021

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
Copy link
Collaborator Author

@andreasnoack does eyeballing that diff suggest something insightful?

@andreasnoack
Copy link
Member

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

@ViralBShah
Copy link
Collaborator Author

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

@fghoussen
Copy link

@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
Copy link

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
Copy link

depend on a bad choice of the initial random vector

@ViralBShah: how to do initialize your vectors?

@ViralBShah
Copy link
Collaborator Author

ViralBShah commented Oct 21, 2022

@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
Copy link

What should we be doing here?

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

@ViralBShah
Copy link
Collaborator Author

ViralBShah commented Oct 21, 2022

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
Copy link

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
Copy link
Collaborator Author

Yes, init with rand does work better.

@amontoison
Copy link

amontoison commented Oct 22, 2022

@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
Copy link

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
Copy link
Collaborator Author

@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
Copy link
Collaborator Author

ViralBShah commented Oct 22, 2022

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
Copy link

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
Copy link
Collaborator Author

@fghoussen
Copy link

fghoussen commented Mar 8, 2023

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
Copy link

fghoussen commented Mar 8, 2023

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
Copy link
Collaborator Author

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
Copy link

fghoussen commented Mar 8, 2023

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
Copy link
Collaborator Author

It is not.

@fghoussen
Copy link

fghoussen commented Mar 10, 2023

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants