JuliaLinearAlgebra / IterativeSolvers.jl

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

accept user-specified rng for random-number generation #303

Open stevengj opened 3 years ago

stevengj commented 3 years ago

Rather than calling rand(...), it would be good to call rand(rng, ...) with an optionally user-specified rng::AbstractRNG keyword argument in algorithms that require pseudo-random initialization, e.g. lobpcg (https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/ae01dfe228138d0db9164eb343252150f6c0bfc6/src/lobpcg.jl#L872), bicgstabl (https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/ae01dfe228138d0db9164eb343252150f6c0bfc6/src/bicgstabl.jl#L38), and idrs (https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/ae01dfe228138d0db9164eb343252150f6c0bfc6/src/bicgstabl.jl#L38).

Rationale: if deterministic results are desired (as is often the case), it is better for the user to be able to specify an RNG than setting the global Random.seed!.

Moreover, arguably the default rng should be a new deterministic RNG object created with a fixed seed every time the iterative solver is initialized, so that the algorithm defaults to being deterministic.

(Note: "deterministic" ≠ "reproducible" … e.g. changing the Julia version could change the results. But it is useful to be able to run the same code with the same software on the same hardware multiple times and get exactly the same results as much as possible.)

antoine-levitt commented 3 years ago

Or simply accept rshadow as a kwarg and let the user do the initialization if needed?

stevengj commented 3 years ago

Different iterative algorithms initialize different pseudorandom internal state. It seems like it's asking a lot of the caller to know what that internal state means for different solvers. Also, it's nice to have a uniform API that works with all of the algorithms that employ pseudorandom state.