Randomized algorithms

reig(A, l)

Compute the spectral (Eigen) decomposition of A using a randomized algorithm.

Arguments

• A: input matrix.

• l::Int: number of eigenpairs to find.

Output

• ::Base.LinAlg.Eigen: eigen decomposition.

Implementation note

This is a wrapper around eigfact_onepass() which uses the randomized samples found using srft(l).

rsvdfact(A, n, p=0)

Compute partial singular value decomposition of A using a randomized algorithm.

Arguments

A: input matrix.

n::Int: number of singular value/vector pairs to find.

p::Int=0: number of extra vectors to include in computation.

Output

::SVD: singular value decomposition.

Implementation note

This function calls rrange, which uses naive randomized rangefinding to compute a basis for a subspace of dimension n (Algorithm 4.1 of [Halko2011]), followed by svdfact_restricted(), which computes the exact SVD factorization on the restriction of A to this randomly selected subspace (Algorithm 5.1 of [Halko2011]).

Alternatively, you can mix and match your own randomized algorithm using any of the randomized range finding algorithms to find a suitable subspace and feeding the result to one of the routines that computes the SVD restricted to that subspace.

rsvd_fnkz(A, k)

Compute the randomized SVD by iterative refinement from randomly selected columns/rows [Friedland2006].

Arguments

• A: matrix whose SVD is desired;

• k::Int: desired rank of approximation (k ≤ min(m, n)).

Keywords

• l::Int = k: number of columns/rows to sample at each iteration (1 ≤ l ≤ k);

• N::Int = minimum(size(A)): maximum number of iterations;

• ϵ::Real = prod(size(A))*eps(): relative threshold for convergence, as measured by growth of the spectral norm;

• method::Symbol = :eig: problem to solve.

  1. :eig: eigenproblem.

  2. :svd: singular problem.

• verbose::Bool = false: print convergence information at each iteration.

Return value

SVD object of rank ≤ k.

rcond(A, iters=1)

Estimate matrix condition number randomly.

Arguments

• A: matrix whose condition number to estimate. Must be square and

support premultiply (A*⋅) and solve (A\⋅).

• iters::Int = 1: number of power iterations to run.

Keywords

• p::Real = 0.05: probability that estimate fails to hold as an upper bound.

Output

Interval (x, y) which contains κ(A) with probability 1 - p.

Implementation note

[Dixon1983] originally describes this as a computation that can be done by computing the necessary number of power iterations given p and the desired accuracy parameter θ=y/x. However, these bounds were only derived under the assumptions of exact arithmetic. Empirically, iters≥4 has been seen to result in incorrect results in that the computed interval does not contain the true condition number. This implemention therefore makes iters an explicitly user-controllable parameter from which to infer the accuracy parameter and hence the interval containing κ(A). ```

reigmin(A, iters=1)

Estimate minimal eigenvalue randomly.

Arguments

• A: Matrix whose maximal eigenvalue to estimate.

Must be square and support premultiply (A*⋅).

• iters::Int=1: Number of power iterations to run. (Recommended: iters ≤ 3)

Keywords

• p::Real=0.05: Probability that estimate fails to hold as an upper bound.

Output

Interval (x, y) which contains the maximal eigenvalue of A with probability 1 - p.

References

Corollary of Theorem 1 in [Dixon1983]

reigmax(A, iters=1)

Estimate maximal eigenvalue randomly.

Arguments

• A: Matrix whose maximal eigenvalue to estimate.

Must be square and support premultiply (A*⋅).

• iters::Int=1: Number of power iterations to run. (Recommended: iters ≤ 3)

Keywords

• p::Real=0.05: Probability that estimate fails to hold as an upper bound.

Output

Interval (x, y) which contains the maximal eigenvalue of A with probability 1 - p.

References

Corollary of Theorem 1 in [Dixon1983]

rnorm(A, mvps)

Compute a probabilistic upper bound on the norm of a matrix A. ‖A‖ ≤ α √(2/π) maxᵢ ‖Aωᵢ‖ with probability p=α^(-mvps).

Arguments

• A: matrix whose norm to estimate.

• mvps::Int: number of matrix-vector products to compute.

Keywords

• p::Real=0.05: probability of upper bound failing.

Output

Estimate of ‖A‖.

See also rnorms for a different estimator that uses premultiplying by both A and A'.

References

Lemma 4.1 of Halko2011

rnorms(A, iters=1)

Estimate matrix norm randomly using A'A.

Compute a probabilistic upper bound on the norm of a matrix A.

ρ = √(‖(A'A)ʲω‖/‖(A'A)ʲ⁻¹ω‖)

which is an estimate of the spectral norm of A produced by iters steps of the power method starting with normalized ω, is a lower bound on the true norm by a factor

ρ ≤ α ‖A‖

with probability greater than 1 - p, where p = 4\sqrt(n/(iters-1)) α^(-2iters).

Arguments

• A: matrix whose norm to estimate.

• iters::Int = 1: mumber of power iterations to perform.

Keywords

• p::Real = 0.05: probability of upper bound failing.

• At = A': Transpose of A.

Output

Estimate of ‖A‖.

See also rnorm for a different estimator that does not require premultiplying by A'

References

Appendix of [Liberty2007].