IterativeSolvers.jl icon indicating copy to clipboard operation
IterativeSolvers.jl copied to clipboard

use deterministic RNG initial state

Open stevengj opened this issue 4 years ago • 7 comments

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 avatar Oct 30 '21 12:10 stevengj

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

stevengj avatar Oct 30 '21 20:10 stevengj

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...

stevengj avatar Oct 30 '21 20:10 stevengj

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

mschauer avatar Nov 01 '21 11:11 mschauer

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?)

stevengj avatar Nov 01 '21 13:11 stevengj

Ping @haampie

mschauer avatar Dec 14 '21 08:12 mschauer

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

mschauer avatar Jan 17 '22 17:01 mschauer

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?

mschauer avatar Sep 01 '22 14:09 mschauer