JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
394 stars 106 forks source link

use deterministic RNG initial state #304

Open stevengj opened 2 years ago

stevengj commented 2 years ago

Closes #303.

Uses a deterministic RNG of MersenneTwister(seed) for generating initial vectors in iterative algorithms that use random initial vectors, controllable via an rng keyword.

(In some cases, the rng keyword is redundant, because the algorithm alternatively allows you to pass the initial vector directly, but it seemed better to me to have a consistent API that can be used with any of the iterative solvers doing pseudorandom initialization.)

(Bumped Julia requirement to 1.4 to support broadcasting RNG objects.)

stevengj commented 2 years ago

Weird, I can't replicate the test failure locally…

stevengj commented 2 years ago

Okay, now I can replicate a failure. It looks like one issue was that lobpcg! calls itself recursively, but I wasn't passing through the rng keyword, so it was generating the same pseudorandom numbers on each recursive call. But there is still some other failure I need to track down...

mschauer commented 2 years ago

The commit a4db20b had some build's passing, are the errors just from having different seeds and too tights bounds for stochastic tests?

stevengj commented 2 years ago

are the errors just from having different seeds and to tights bounds for stochastic tests?

Looks like it. I switched to master and changed the seed! line in test/lobpcg.jl to Random.seed!(Random.make_seed()), to ensure a different seed on every run, and now I'm getting the same failures about 25% of the time — I include that file 100 times, and it failed 26 times.

Note that one of the failures is not just a tolerance thing, but a robustness problem: the Cholesky factorization here is sometimes failing because the matrix is indefinite (a nearly zero eigenvalue that is slightly indefinite … basically it looks like X is rank-deficient, in which case you should probably do some re-initialization?).

The other failure is on

test/lobpcg.jl:316
  Expression: all(isapprox.(adjoint(X) * (B * X), Matrix{T}(I, 3, 3), atol = 2 * n * tol))

where it looks like the tolerance might be too small.

(Not sure why you are using all(isapprox.(..)) instead of just isapprox on the whole matrix here?)

mschauer commented 2 years ago

Ping @haampie

mschauer commented 2 years ago

@stevengj What to do with the tests? I agree that we can soften a bound for the second failure.

mschauer commented 1 year ago

There is one test passing,

CI / Julia 1.6 - ubuntu-latest - x64 - pull_request (pull_request) which seems to be that the lobpg issue is fixed on some level. But shouldn't it pass on all systems as it is deterministic?