MatrixFreeRandomizedLinearAlgebra.jl

Randomized algorithms for the top singular values, eigenvalues, and the trace of a linear operator, without building the matrix.

If you can multiply your operator by a vector (and, for some of these, by its adjoint), you can use it here: dense or sparse matrices, GPU arrays, LinearMaps.jl operators, or any type that defines size and *. The methods are randomized approximations from Halko, Martinsson, and Tropp. They are most useful when the operator is large, a matrix-vector product is cheap, and you only need a low-rank piece of it, which is the case where a full dense factorization is too slow or doesn't fit in memory.

Left: time to compute the leading singular values of a matrix-free 2D blur operator, dense versus this package. Right: the values it returns sit on top of the exact spectrum. Both plots come from the scripts in examples/.

What's here

• rsvd / rsvdvals: randomized SVD and singular values for general (possibly rectangular) operators.
• reigen_hermitian / reigvals_hermitian: randomized eigendecomposition and eigenvalues for Hermitian operators.
• trace: stochastic trace estimation. XTrace by default, or a streaming Hutchinson estimator when memory is tight.

All of it runs on CPU or GPU arrays, and the decomposition routines take a seed_Q keyword if you already have a basis to start from.

Installation

] add MatrixFreeRandomizedLinearAlgebra

Quick example

using MatrixFreeRandomizedLinearAlgebra, LinearAlgebra

A = randn(100, 50)            # any operator with size, *, and '
U, S, Vt = rsvd(A, 10)        # rank-10 randomized SVD
rel = opnorm(A - U * Diagonal(S) * Vt) / opnorm(A)
println("relative error: ", rel)

Learn more

• Documentation: what matrix-free means, how the algorithms work, and the API reference.
• examples/: runnable scripts that compare accuracy and runtime against dense references. They generate the plots above.