Arpack.jl

Arpack.jl copied to clipboard

I am trying to use svds but svds for same matrix is giving different singular vectors after every run. Due to this reason my codes are non-reproducible. I also checked in Matlab, there svds gives exactly same result after every run.

So below I have show a small example for real and complex matrix. For real matrix singular vectors are same but for some of the vectors sign is flipped. For complex matrices singular vectors are completely different after every run of svds.

Please advice if anything can be done to fix this issue.

Note that the singular values are the same, and the factorization is correct. This is what we get from ARPACK, and is not wrong per se. Although, it would be nice to have these factorizations be more reproducible. PRs for post-processing and making the results reproducible across runs are welcome.

ARPACK uses a random starting vector by default (generated with this cgetv0 function, which calls the LAPACK clarnv function). This causes the phase of the singular vectors to be random (or just the sign, for real values).

The resulting singular vectors are still correct (the phase is arbitrary in the SVD). However, if you want a deterministic phase, you can use the v0 keyword parameter to pass your own initial guess to svds

In the long run, it would probably be better for Arpack.jl to generate its own random numbers here. That way:

1. It would be thread-safe — I'm not sure if LAPACK's random numbers are thread-safe?
2. It would allow you to get a reproducible result more easily by setting Julia's random seed.

(However, we would have to replicate some of the logic from cgetv0.)

(In the even longer run, it would probably be better to point people towards a pure-Julia Arnoldi implementation; I'm not sure if the current ones have feature parity with ARPACK yet, though?)

The main reason I updated this package was because I wasn't sure if the pure Julia implementations were on feature parity. Also, it is good to have the ARPACK implementation to test against, compare performance, etc.

I'll take some of these comments and add them to the documentation.

Documented in #116. Leaving the issue open in case someone wants to try their hand at the proposed solution.

Note that Julia now has a thread safe RNG.