GitHub

Sketching

Sketching methods project a large matrix onto a low-dimensional subspace before any expensive dense factorization is performed. The result is an approximate low-rank decomposition that can be substantially cheaper than the corresponding full decomposition, depending on the quality of the sketch and the spectrum of the original matrix.

This page first describes how to compute a sketch as a standalone operation, then how to plug a sketch into a partial decomposition such as a truncated SVD.

Standalone Sketches

The basic sketching primitives are left_sketch and right_sketch (with their bang counterparts). They return an isometric range/co-range factor together with a corresponding core matrix, independent of any subsequent decomposition.

left_sketch computes an isometric matrix Q whose column span approximates the range of A, together with the core factor B = Q' * A:

using MatrixAlgebraKit
using MatrixAlgebraKit: diagview
using LinearAlgebra: norm
using Random: MersenneTwister

# A rank-3 matrix
A = randn(MersenneTwister(0), 8, 3) * randn(MersenneTwister(1), 3, 6);

Q, B = left_sketch(A; howmany = 3, rng = MersenneTwister(42));
isisometric(Q) && A ≈ Q * B

right_sketch is the dual operation, returning a right-isometric matrix Pᴴ (orthonormal rows) and core factor B = A * Pᴴ':

B, Pᴴ = right_sketch(A; howmany = 3, rng = MersenneTwister(42));
isisometric(Pᴴ') && A ≈ B * Pᴴ

These functions follow the same conventions as other decompositions: left_sketch! / right_sketch! may destroy the input, and the bang forms optionally accept pre-allocated output tuples.

MatrixAlgebraKit.left_sketch — Function
left_sketch(A; howmany, kwargs...) -> Q, B
left_sketch(A, alg::AbstractAlgorithm) -> Q, B
left_sketch!(A, [QB]; howmany, kwargs...) -> Q, B
left_sketch!(A, [QB], alg::AbstractAlgorithm) -> Q, B

Compute an isometric matrix Q (orthonormal columns) of size m×k, whose column span approximates the range of A of size m×n. Also create the core factor B = Q' * A. Here k = howmany is the sketch dimension.

The keyword arguments construct a GaussianSketching strategy unless an explicit alg::SketchingStrategy is supplied. howmany is required.

Note

The bang method left_sketch! optionally accepts the output matrices Q, B and possibly destroys the input matrix A. Always use the return value of the function as it may not always be possible to use the provided Q, B as output.

See also right_sketch(!) and SketchedAlgorithm.

source
MatrixAlgebraKit.left_sketch! — Function
left_sketch(A; howmany, kwargs...) -> Q, B
left_sketch(A, alg::AbstractAlgorithm) -> Q, B
left_sketch!(A, [QB]; howmany, kwargs...) -> Q, B
left_sketch!(A, [QB], alg::AbstractAlgorithm) -> Q, B

Compute an isometric matrix Q (orthonormal columns) of size m×k, whose column span approximates the range of A of size m×n. Also create the core factor B = Q' * A. Here k = howmany is the sketch dimension.

The keyword arguments construct a GaussianSketching strategy unless an explicit alg::SketchingStrategy is supplied. howmany is required.

Note

The bang method left_sketch! optionally accepts the output matrices Q, B and possibly destroys the input matrix A. Always use the return value of the function as it may not always be possible to use the provided Q, B as output.

See also right_sketch(!) and SketchedAlgorithm.

source
MatrixAlgebraKit.right_sketch — Function
right_sketch(A; howmany, kwargs...) -> B, Pᴴ
right_sketch(A, alg::AbstractAlgorithm) -> B, Pᴴ
right_sketch!(A, [BPᴴ]; howmany, kwargs...) -> B, Pᴴ
right_sketch!(A, [BPᴴ], alg::AbstractAlgorithm) -> B, Pᴴ

Compute a right-isometric matrix Pᴴ (orthonormal rows) of size k×n, whose row span approximates the range of A of size m×n. Also create the core factor B = A * Pᴴ' Here k = howmany is the sketch dimension.

The keyword arguments construct a GaussianSketching strategy unless an explicit alg::SketchingStrategy is supplied. howmany is required.

Note

The bang method right_sketch! optionally accepts the output matrices BPᴴ and possibly destroys the input matrix A. Always use the return value of the function as it may not always be possible to use the provided BPᴴ as output.

See also left_sketch(!) and SketchedAlgorithm.

source
MatrixAlgebraKit.right_sketch! — Function
right_sketch(A; howmany, kwargs...) -> B, Pᴴ
right_sketch(A, alg::AbstractAlgorithm) -> B, Pᴴ
right_sketch!(A, [BPᴴ]; howmany, kwargs...) -> B, Pᴴ
right_sketch!(A, [BPᴴ], alg::AbstractAlgorithm) -> B, Pᴴ

Compute a right-isometric matrix Pᴴ (orthonormal rows) of size k×n, whose row span approximates the range of A of size m×n. Also create the core factor B = A * Pᴴ' Here k = howmany is the sketch dimension.

The keyword arguments construct a GaussianSketching strategy unless an explicit alg::SketchingStrategy is supplied. howmany is required.

Note

The bang method right_sketch! optionally accepts the output matrices BPᴴ and possibly destroys the input matrix A. Always use the return value of the function as it may not always be possible to use the provided BPᴴ as output.

See also left_sketch(!) and SketchedAlgorithm.

source

Available Sketching Strategies

The keyword arguments accepted by left_sketch / right_sketch are forwarded to the default sketching strategy for the input type (currently GaussianSketching for all AbstractMatrix). For full control, construct a strategy directly and pass it as the second positional argument:

Q, B = left_sketch(A, GaussianSketching(3; numiter = 4, rng = MersenneTwister(42)));
A ≈ Q * B
MatrixAlgebraKit.GaussianSketching — Type
GaussianSketching(howmany; numiter, rng)

Sketching strategy using a Gaussian random matrix with optional power iterations to improve accuracy on slowly-decaying spectra.

Fields

  • howmany::Int: number of singular values to compute.
  • numiter::Int: number of power iterations (numiter ≥ 1; the first counts as the initial sketch).
  • rng::AbstractRNG: random number generator used to draw the Gaussian sketch matrix.
source

The numiter keyword controls the number of power iterations. The first iteration is the initial sketch; additional iterations apply A * A' to improve accuracy on slowly-decaying spectra at the cost of two extra matrix multiplications per iteration. The default numiter = 2 is a conservative choice; values of 4–5 often improve accuracy significantly when the singular values decay slowly.

Additional strategies

GaussianSketching is currently the only built-in SketchingStrategy. Additional strategies (for example, structured or subsampled sketches) may be added in the future; the interface is deliberately written against the abstract SketchingStrategy supertype so that new strategies plug in without changes to the downstream decomposition code.

Sketched Partial Decompositions

A sketch can be combined with a small dense decomposition of the core factor to obtain an approximate partial decomposition of the original matrix. At present, this is supported for the truncated SVD via svd_trunc / svd_trunc! / svd_trunc_no_error.

Not yet supported

Sketched variants of eigh_trunc and eig_trunc are natural extensions of the same machinery but are not implemented yet. The SketchedAlgorithm wrapper and the sketch = keyword are currently only accepted by the truncated SVD functions.

There are two equivalent ways to request a sketched truncated SVD, paralleling the two-form syntax used for Truncations.

1. Using the sketch keyword with a NamedTuple

The simplest form is to pass a NamedTuple of sketch parameters together with the desired truncation:

U, S, Vᴴ, ϵ = svd_trunc(A;
    sketch = (; howmany = 3, rng = MersenneTwister(42)),
    trunc = truncrank(3),
);
size(diagview(S), 1) == 3 && A ≈ U * S * Vᴴ

The NamedTuple keywords are forwarded to the default sketching strategy for the input type, exactly as for left_sketch above.

2. Using an explicit SketchedAlgorithm

For full control, construct a SketchedAlgorithm value directly and pass it as the alg argument:

alg = SketchedAlgorithm(;
    sketch = GaussianSketching(3; rng = MersenneTwister(42)),
    trunc = truncrank(3),
);
U, S, Vᴴ, ϵ = svd_trunc(A, alg);
A ≈ U * S * Vᴴ

When an alg::SketchedAlgorithm is supplied, the sketch and trunc keywords cannot also be specified at the call site; doing so raises ArgumentError. All configuration must instead live inside the algorithm constructor.

The SketchedAlgorithm Wrapper

MatrixAlgebraKit.SketchedAlgorithm — Type
SketchedAlgorithm(;
    alg::AbstractAlgorithm, sketch::SketchingStrategy,
    trunc::TruncationStrategy, driver::Driver = DefaultDriver()
)

Generic wrapper type for self-truncating algorithms that produce an approximate low-rank factorization by first applying a sketching operation specified by sketch, then computing a small dense decomposition of the projected matrix using alg. The driver selects the backend implementing the sketched factorization (e.g. Native() for the generic sketch-then-decompose pipeline, CUSOLVER() for the fused gesvdr kernel).

source

SketchedAlgorithm differs from TruncatedAlgorithm in that it is self-truncating: the sketch step itself produces a small dense problem of size sketch.howmany, and any further trunc is applied to the result of the inner factorization rather than to a full dense decomposition.

The driver field selects the backend implementing the sketched pipeline:

  • Native() (the default for CPU array types) runs the generic sketch-then-decompose pipeline using the standard left_sketch! / right_sketch! building blocks followed by the inner alg.
  • CUSOLVER() (the default for CUDA array types, with the appropriate extension loaded) dispatches to cuSOLVER's fused gesvdr kernel, which performs the sketch and the small SVD in a single device call.

To force a particular driver, set it explicitly:

using MatrixAlgebraKit: CUSOLVER  # driver types are not exported by default

alg = SketchedAlgorithm(;
    sketch = GaussianSketching(k; numiter = 4),
    trunc = truncrank(k),
    driver = CUSOLVER(),
)
U, S, Vᴴ, ϵ = svd_trunc(A_cuda, alg)