RandomizedPreconditioners

This package contains several randomized preconditioners, which use randomized numerical linear algebra to construct approximate inverses of matrices. These approximate inverses can dramatically speed up iterative linear system solvers.

Given a positive semidefinite matrix A, the Nyström Sketch  ≈ A is constructed by

using RandomizedPreconditioners
 = NystromSketch(A, k, r)

where k and r are parameters with k ≤ r.

We can use Â to construct a preconditioner P ≈ A + μ*I for the system (A + μ*I)x = b, which is solved by conjugate gradients.

If you need P (e.g., IterativeSolvers.jl), use

P = NystromPreconditioner(Â, μ)

If you need P⁻¹ (e.g., Krylov.jl), use

Pinv = NystromPreconditionerInverse(Anys, μ)

These preconditioners can be simply passed into the solvers, for example

using Krylov
x, stats = cg(A+μ*I, b; M=Pinv)

The package LinearSolve.jl defines a convenient common interface to access all the Krylov implementations, which makes testing very easy.

using RandomizedPreconditioners, LinearSolve
 = NystromSketch(A, k, r)
P = NystromPreconditioner(Â, μ)

prob = LinearProblem(A, b)
sol = solve(prob, IterativeSolversJL_CG(), Pl=P)

The sketching algorithms below use a rank r sketching matrix and have complexity O(n²r). The parameter k ≤ r truncates the sketch, which can improve numerical performance. Possible choices include k = r - 10 and k = round(Int, 0.95*r). Sketches allow for faster (approximate) multiplication (* and mul!) and are used to construct preconditioners.

using RandomizedPreconditioners
 = NystromSketch(A, k, r)

using RandomizedPreconditioners
 = EigenSketch(A, k, r)

The Randomized SVD uses the powered randomized rangefinder [2, Alg. 9] with powering parameter q. Small values of q (e.g., 5) seem to perform well. Note that the complexity increases to O(n²rq).

using RandomizedPreconditioners
 = RandomizedSVD(A, k, r; q=10)

We implement two algorithms for a randomized estimate of the maximum eigenvalue for a PSD matrix: the power method and the Lanczos method.

using RandomizedPreconditioners
const RP = RandomizedPreconditioners

λmax_power =  RP.eigmax_power(A)
λmax_lanczos = RP.eigmax_lanczos(A)
λmin_lanczos = RP.eigmin_lanczos(A)

The Lanczos method can estimate the maximum and minimum eigenvalue simultaneously:

λmax, λmin = RP.eig_lanczos(A; eigtype=0)

There are several choices for the random embedding used in the algorithms. By default, this package uses Gaussian embeddings (and Gaussian test matrices), but it also includes the SSFT and the ability to add new test matrices by implementing the TestMatrix interface.

A TestMatrix, Ω, should implement matrix multiplication for itself and its adjoint by implementing the !mul method. See Martinsson and Tropp [2] Section 9 for more.

• Test Matrices
  • TestMatrix type
  • Sparse maps
  • Subsampled trigonometric transform
  • DCT & Hartley option for SSRFT
  • Tensor product maps
• Rangefinders
  • Lanzcos randomized rangefinder
  • Chebyshev randomized rangefinder
  • Incremental rangefinder with updating
  • Subsequent orthogonalization
  • A posteriori error estimation
  • Incremental rangefinder with powering
  • Incremental rangefinder for sparse matrices
• Sketches & Factorizations
  • Powering option / incorporating rangefinder into Nystrom sketch
  • powerURV (w. re-orthonormalization)
  • CPQR decomposition
  • Improve randomized SVD
• Preconditioners
  • Add preconditioner for symmetric systems
  • Preconditioner for least squares
• Performance
  • Avoid redoing computations in adaptive sketch
  • General performance
• Documentation
  • More complete general docs
  • Least squares example (sketch & solve, iterative sketching, sketch & precondition)

[1] Zachary Frangella, Joel A Tropp, and Madeleine Udell. “Randomized Nyström Preconditioning.” In:arXiv preprint arXiv:2110.02820(2021). https://arxiv.org/abs/2110.02820

[2] PG Martinsson and JA Tropp. “Randomized numerical linear algebra: foundations & algorithms (2020).” In: arXiv preprint arXiv:2002.01387. https://arxiv.org/abs/2002.01387