use deterministic RNG initial state

Author: stevengj

Project.toml

@@ -11,4 +11,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 RecipesBase = "0.6, 0.7, 0.8, 1.0"
-julia = "1.3"
+julia = "1.4"

docs/src/eigenproblems/lobpcg.md

@@ -13,6 +13,10 @@ lobpcg!
 
 A `LOBPCGIterator` is created to pre-allocate all the memory required by the method using the constructor `LOBPCGIterator(A, B, largest, X, P, C)` where `A` and `B` are the matrices from the generalized eigenvalue problem, `largest` indicates if the problem is a maximum or minimum eigenvalue problem, `X` is the initial eigenbasis, randomly sampled if not input, where `size(X, 2)` is the block size `bs`. `P` is the preconditioner, `nothing` by default, and `C` is the constraints matrix. The desired `k` eigenvalues are found `bs` at a time.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.
 
 ## References
 Implementation is based on [^Knyazev1993] and [^Scipy].

docs/src/linear_systems/bicgstabl.md

@@ -15,11 +15,16 @@ The method is based on the original article [^Sleijpen1993], but does not implem
 
 The `r` and `u` factors are pre-allocated as matrices of size $n \times (l + 1)$, so that BLAS2 methods can be used. Also the random shadow residual is pre-allocated as a vector. Hence the storage costs are approximately $2l + 3$ vectors.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.
+
 !!! tip
     BiCGStabl(l) can be used as an [iterator](@ref Iterators).
 
-[^Sleijpen1993]: 
+[^Sleijpen1993]:
 
-    Sleijpen, Gerard LG, and Diederik R. Fokkema. "BiCGstab(l) for 
-    linear equations involving unsymmetric matrices with complex spectrum." 
+    Sleijpen, Gerard LG, and Diederik R. Fokkema. "BiCGstab(l) for
+    linear equations involving unsymmetric matrices with complex spectrum."
     Electronic Transactions on Numerical Analysis 1.11 (1993): 2000.

docs/src/linear_systems/idrs.md

@@ -14,4 +14,9 @@ idrs!
 The current implementation is based on the [MATLAB version](http://ta.twi.tudelft.nl/nw/users/gijzen/idrs.m) by Van Gijzen and Sonneveld. For background see [^Sonneveld2008], [^VanGijzen2011] and [the IDR(s) webpage](http://ta.twi.tudelft.nl/nw/users/gijzen/IDR.html).
 
 [^Sonneveld2008]: IDR(s): a family of simple and fast algorithms for solving large nonsymmetric linear systems. P. Sonneveld and M. B. van Gijzen SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035--1062, 2008
-[^VanGijzen2011]: Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits Bi-orthogonality Properties. M. B. van Gijzen and P. Sonneveld ACM Trans. Math. Software,, Vol. 38, No. 1, pp. 5:1-5:19, 2011
+[^VanGijzen2011]: Algorithm 913: An Elegant IDR(s) Variant that Efficiently Exploits Bi-orthogonality Properties. M. B. van Gijzen and P. Sonneveld ACM Trans. Math. Software,, Vol. 38, No. 1, pp. 5:1-5:19, 2011
+
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument.

docs/src/svd/svdl.md

@@ -14,4 +14,10 @@ The implementation of thick restarting follows closely that of SLEPc as describe
 
 The singular vectors are computed directly by forming the Ritz vectors from the product of the Lanczos vectors `L.P`/`L.Q` and the singular vectors of `L.B`. Additional accuracy in the singular triples can be obtained using inverse iteration.
 
+A deterministic seed is used for generating pseudo-random initial
+data for the algorithm; this can be controlled by passing a
+different pseudorandom number generator (an [`AbstractRNG`](https://docs.julialang.org/en/v1/stdlib/Random/#Random.AbstractRNG)) via
+the `rng` keyword argument, or by passing an initial `v0` vector
+directly.
+
 [^Hernandez2008]: Vicente Hernández, José E. Román, and Andrés Tomás. "A robust and efficient parallel SVD solver based on restarted Lanczos bidiagonalization." Electronic Transactions on Numerical Analysis 31 (2008): 68-85.

src/bicgstabl.jl

@@ -1,5 +1,5 @@
 export bicgstabl, bicgstabl!, bicgstabl_iterator, bicgstabl_iterator!, BiCGStabIterable
-using Printf
+using Printf, Random
 import Base: iterate
 
 mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, smallMatT <: AbstractMatrix, realT <: Real, scalarT <: Number}
@@ -29,13 +29,14 @@ function bicgstabl_iterator!(x, A, b, l::Int = 2;
                              max_mv_products = size(A, 2),
                              abstol::Real = zero(real(eltype(b))),
                              reltol::Real = sqrt(eps(real(eltype(b)))),
+                             rng::AbstractRNG=MersenneTwister(0),
                              initial_zero = false)
     T = eltype(x)
     n = size(A, 1)
     mv_products = 0
 
     # Large vectors.
-    r_shadow = rand(T, n)
+    r_shadow = rand(rng, T, n)
     rs = Matrix{T}(undef, n, l + 1)
     us = zeros(T, n, l + 1)
 
@@ -156,6 +157,7 @@ bicgstabl(A, b, l = 2; kwargs...) = bicgstabl!(zerox(A, b), A, b, l; initial_zer
 - `max_mv_products::Int = size(A, 2)`: maximum number of matrix vector products.
 For BiCGStab(l) this is a less dubious term than "number of iterations";
 - `Pl = Identity()`: left preconditioner of the method;
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 - `abstol::Real = zero(real(eltype(b)))`,
   `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
   tolerance for the stopping condition

src/idrs.jl

@@ -32,6 +32,7 @@ shadow space.
   is the residual in the `k`th iteration;
 - `maxiter::Int = size(A, 2)`: maximum number of iterations;
 - `log::Bool`: keep track of the residual norm in each iteration;
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 - `verbose::Bool`: print convergence information during the iterations.
 
 # Return values
@@ -80,8 +81,8 @@ end
 end
 
 function idrs_method!(log::ConvergenceHistory, X, A, C::T,
-    s::Number, Pl::precT, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false
-    ) where {T, precT}
+    s::Number, Pl::precT, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false,
+    rng::AbstractRNG=MersenneTwister(0)) where {T, precT}
 
     verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
     R = C - A*X
@@ -102,7 +103,7 @@ function idrs_method!(log::ConvergenceHistory, X, A, C::T,
 
     Z = zero(C)
 
-    P = T[rand!(copy(C)) for k in 1:s]
+    P = T[rand!(rng, copy(C)) for k in 1:s]
     U = T[copy(Z) for k in 1:s]
     G = T[copy(Z) for k in 1:s]
     Q = copy(Z)

src/lobpcg.jl

@@ -787,8 +787,8 @@ Finds the `nev` extremal eigenvalues and their corresponding eigenvectors satisf
 function lobpcg(A, largest::Bool, nev::Int; kwargs...)
     lobpcg(A, nothing, largest, nev; kwargs...)
 end
-function lobpcg(A, B, largest::Bool, nev::Int; kwargs...)
-    lobpcg(A, B, largest, rand(eltype(A), size(A, 1), nev); not_zeros=true, kwargs...)
+function lobpcg(A, B, largest::Bool, nev::Int; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+    lobpcg(A, B, largest, rand(rng, eltype(A), size(A, 1), nev); not_zeros=true, rng=rng, kwargs...)
 end
 
 """
@@ -804,6 +804,7 @@ end
 ## Keywords
 
 - `not_zeros`: default is `false`. If `true`, `X0` will be assumed to not have any all-zeros column.
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 
 - `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
     states; if `false` only `results.trace` will be empty;
@@ -826,6 +827,7 @@ function lobpcg(A, largest::Bool, X0; kwargs...)
 end
 function lobpcg(A, B, largest, X0;
                 not_zeros=false, log=false, P=nothing, maxiter=200,
+                rng::AbstractRNG=MersenneTwister(0),
                 C=nothing, tol::Real=default_tolerance(eltype(X0)))
     X = copy(X0)
     n = size(X, 1)
@@ -835,7 +837,7 @@ function lobpcg(A, B, largest, X0;
 
     iterator = LOBPCGIterator(A, B, largest, X, P, C)
 
-    return lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
+    return lobpcg!(iterator, rng=rng, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)
 end
 
 """
@@ -849,6 +851,7 @@ end
 ## Keywords
 
 - `not_zeros`: default is `false`. If `true`, the initial Ritz vectors will be assumed to not have any all-zeros column.
+- `rng::AbstractRNG`: generator for pseudorandom initialization
 
 - `log::Bool`: default is `false`; if `true`, `results.trace` will store iterations
     states; if `false` only `results.trace` will be empty;
@@ -863,21 +866,22 @@ end
 
 """
 function lobpcg!(iterator::LOBPCGIterator; log=false, maxiter=200, not_zeros=false,
+                 rng::AbstractRNG=MersenneTwister(0),
                  tol::Real=default_tolerance(eltype(iterator.XBlocks.block)))
     X = iterator.XBlocks.block
     iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
     if !not_zeros
         @views for j in 1:size(X,2)
             if all(x -> x==0, X[:, j])
-                @inbounds X[:,j] .= rand.()
+                @inbounds X[:,j] .= rand.(rng)
             end
         end
         iterator.constr!(iterator.XBlocks.block, iterator.tempXBlocks.block)
     end
     n = size(X, 1)
     sizeX = size(X, 2)
     iterator.iteration[] = 1
-    while iterator.iteration[] <= maxiter 
+    while iterator.iteration[] <= maxiter
         state = iterator(tol, log)
         if log
             push!(iterator.trace, state)
@@ -927,6 +931,7 @@ function lobpcg(A, largest::Bool, X0, nev::Int; kwargs...)
 end
 function lobpcg(A, B, largest::Bool, X0, nev::Int;
                 not_zeros=false, log=false, P=nothing, maxiter=200,
+                rng::AbstractRNG=MersenneTwister(0),
                 C=nothing, tol::Real=default_tolerance(eltype(X0)))
     n = size(X0, 1)
     sizeX = size(X0, 2)
@@ -946,13 +951,13 @@ function lobpcg(A, B, largest::Bool, X0, nev::Int;
             cutoff = sizeX-(nev-converged_x)
             update!(iterator.constr!, iterator.XBlocks.block[:, 1:cutoff], iterator.XBlocks.B_block[:, 1:cutoff])
             X[:, 1:sizeX-cutoff] .= X[:, cutoff+1:sizeX]
-            rand!(X[:, cutoff+1:sizeX])
+            rand!(rng, X[:, cutoff+1:sizeX])
             rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
             append!(r, rnext, converged_x, sizeX-cutoff)
             converged_x += sizeX-cutoff
         else
             update!(iterator.constr!, iterator.XBlocks.block, iterator.XBlocks.B_block)
-            rand!(X)
+            rand!(rng, X)
             rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=true)
             append!(r, rnext, converged_x)
             converged_x += sizeX

src/simple.jl

@@ -56,13 +56,14 @@ function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxit
 end
 
 """
-    powm(B; kwargs...) -> λ, x, [history]
+    powm(B; rng, kwargs...) -> λ, x, [history]
 
 See [`powm!`](@ref). Calls `powm!(B, x0; kwargs...)` with
-`x0` initialized as a random, complex unit vector.
+`x0` initialized as a random, complex unit vector using
+`rng::AbstractRNG`.
 """
-function powm(B; kwargs...)
-    x0 = rand(Complex{real(eltype(B))}, size(B, 1))
+function powm(B; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+    x0 = rand(rng, Complex{real(eltype(B))}, size(B, 1))
     rmul!(x0, one(eltype(B)) / norm(x0))
     powm!(B, x0; kwargs...)
 end
@@ -149,13 +150,13 @@ function powm!(B, x;
 end
 
 """
-    invpowm(B; shift = σ, kwargs...) -> λ, x, [history]
+    invpowm(B; shift = σ, rng, kwargs...) -> λ, x, [history]
 
 Find the approximate eigenpair `(λ, x)` of ``A`` near `shift`, where `B`
 is a linear map that has the effect `B * v = inv(A - σI) * v`.
 
 The method calls `powm!(B, x0; inverse = true, shift = σ)` with
-`x0` a random, complex unit vector. See [`powm!`](@ref)
+`x0` a random, complex unit vector using `rng::AbstractRNG`. See [`powm!`](@ref)
 
 # Examples
 
@@ -168,8 +169,8 @@ Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 50, ismutating = true)
 λ, x = invpowm(Fmap, shift = σ, tol = 1e-4, maxiter = 200)
 ```
 """
-function invpowm(B; kwargs...)
-    x0 = rand(Complex{real(eltype(B))}, size(B, 1))
+function invpowm(B; rng::AbstractRNG=MersenneTwister(0), kwargs...)
+    x0 = rand(rng, Complex{real(eltype(B))}, size(B, 1))
     rmul!(x0, one(eltype(B)) / norm(x0))
     invpowm!(B, x0; kwargs...)
 end

src/svdl.jl

@@ -101,6 +101,7 @@ If `log` is set to `true` is given, method will output a tuple `X, L, ch`. Where
 
 - `nsv::Int = 6`: number of singular values requested;
 - `v0 = random unit vector`: starting guess vector in the domain of `A`.
+- `rng::AbstractRNG`: generator for pseudorandom initialization of `v0`
   The length of `q` should be the number of columns in `A`;
 - `k::Int = 2nsv`: maximum number of Lanczos vectors to compute before restarting;
 - `j::Int = nsv`: number of vectors to keep at the end of the restart.
@@ -175,7 +176,7 @@ end
 #########################
 
 function svdl_method!(log::ConvergenceHistory, A, l::Int=min(6, size(A,1)); k::Int=2l,
-    j::Int=l, v0::AbstractVector = Vector{eltype(A)}(randn(size(A, 2))) |> x->rmul!(x, inv(norm(x))),
+    j::Int=l, rng::AbstractRNG=MersenneTwister(0), v0::AbstractVector = Vector{eltype(A)}(randn(rng, size(A, 2))) |> x->rmul!(x, inv(norm(x))),
     maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
     verbose::Bool=false, method::Symbol=:ritz, vecs=:none, dolock::Bool=false)