low_rank_toolbox.randomized

Randomized algorithms for low-rank matrix approximation.

This submodule implements state-of-the-art randomized algorithms for computing low-rank approximations of large matrices, based on the seminal work by Halko, Martinsson, and Tropp (2010).

Functions

rangefinder : Approximate the range of a dense matrix with a target rank adaptive_rangefinder : Approximate the range of a dense matrix with a target tolerance randomized_svd : Implementation of the randomized SVD algorithm with a prescribed rank adaptive_randomized_svd : Adaptive randomized SVD with tolerance-based rank selection generalized_nystrom : Generalized Nyström method for low-rank approximation

References

Author: Benjamin Carrel, University of Geneva

low_rank_toolbox.randomized.adaptive_randomized_svd(X, tol=1e-06, failure_prob=1e-06, max_rank=None, seed=1234, **extra_data)[source]

Adaptive randomized SVD algorithm.

Reference:

“Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions”, Halko, Martinsson and Tropp 2010.

The adaptive randomized SVD computes an approximate SVD of a matrix X with the following steps: 1. Estimate the range with Q = adaptive_rangefinder(X, tol, failure_prob, seed) 2. Form the smaller matrix C = Q^H X 3. Compute the SVD of C: C = U_C S V^H 4. Form the approximate SVD of X: X_approx = Q U_C S V^H

Parameters:
  • X (SVD | LowRankMatrix | ndarray | LinearOperator) – Matrix or operator to be decomposed

  • tol (float, optional) – Target approximation error (default: 1e-6)

  • failure_prob (float, optional) – Failure probability for the approximation error, must be in (0, 1) (default: 1e-6)

  • max_rank (int, optional) – Maximum rank for the approximation (default: None, no limit)

  • seed (int, optional) – Random seed for reproducibility (default: 1234)

  • **extra_data (dict) – Additional arguments passed to SVD constructor

Returns:

Near-optimal approximation of X with ||X - X_approx||_F <= tol with probability >= 1 - failure_prob

Return type:

SVD

low_rank_toolbox.randomized.adaptive_rangefinder(A, tol=1e-06, failure_prob=1e-06, seed=1234)[source]

The adaptive (randomized) rangefinder method. The tolerance is the error made by the approximation space ||A - QQ^H A||_F <= tol The failure probability is the probability that the error is larger than tol

Reference:

“Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions”, Halko, Martinsson and Tropp 2010.

NOTE: For efficiency, the method performs computations in blocks of size r (defined by the failure probability). Blocking allows to use BLAS3 operations.

Parameters:
  • A (LinearOperator) – The matrix to sketch.

  • tol (float, optional) – The tolerance for the approximation (default: 1e-6).

  • failure_prob (float, optional) – The failure probability, must be in (0, 1) (default: 1e-6).

  • seed (int, optional) – The seed for the random number generator (default: 1234).

Returns:

Q – The sketched matrix with orthonormal columns.

Return type:

ndarray

low_rank_toolbox.randomized.generalized_nystrom(X, r, oversampling_params=(10, 15), epsilon=None, seed=1234, **extra_data)[source]

Generalized Nyström method

Reference:

“Fast and stable randomized low-rank matrix approximation” Nakatsukasa, 2019

Approximation with formula

X ~= X J (K^T X J)^{dagger} K^T X

Parameters:
  • X (LinearOperator) – Matrix to approximate (real-valued only, complex matrices not supported)

  • r (int) – Rank of approximation, must be positive

  • oversampling_params (tuple, optional) – Oversampling parameters (p1, p2) for the two sketch matrices (default: (10, 15))

  • epsilon (float, optional) – When given, perform stable GN with epsilon-truncation for SVD (default: None)

  • seed (int, optional) – Random seed for reproducibility (default: 1234)

  • **extra_data (dict) – Additional arguments passed to QuasiSVD constructor

Returns:

Near optimal best rank r approximation of X in QuasiSVD format

Return type:

QuasiSVD

Notes

This method returns a QuasiSVD (not SVD) because the middle matrix S is typically inverted, making it non-diagonal. Convert to SVD if needed:

result = generalized_nystrom(X, r).to_svd()

low_rank_toolbox.randomized.randomized_svd(X, r, p=10, q=0, truncate=True, seed=1234, **extra_data)[source]

Randomized SVD algorithm.

Reference:

“Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions”, Halko, Martinsson and Tropp 2010.

The randomized SVD computes an approximate SVD of a matrix X with the following steps: 1. Estimate the range with Q = rangefinder(X, r, p, q, seed) 2. Form the smaller matrix C = Q^H X 3. Compute the SVD of C: C = U_C S V^H 4. Form the approximate SVD of X: X_approx = Q U_C S V^H

Parameters:
  • X (SVD | LowRankMatrix | ndarray | LinearOperator) – Matrix or operator to be decomposed

  • r (int) – Rank of the decomposition, must be positive

  • p (int, optional) – Oversampling parameter (default: 10)

  • q (int, optional) – Number of power iterations (default: 0)

  • truncate (bool, optional) – If True, truncate the SVD to rank r (default: True)

  • seed (int, optional) – Random seed for reproducibility (default: 1234)

  • **extra_data (dict) – Additional arguments passed to SVD constructor. If ‘Omega’ is provided, use it as the sketching matrix.

Returns:

Near-optimal best rank r approximation of X

Return type:

SVD

low_rank_toolbox.randomized.rangefinder(A, r, p=5, q=0, seed=1234, **extra_args)[source]

The (randomized) rangefinder method, also called HMT method.

Reference:

“Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions”, Halko, Martinsson and Tropp 2010.

The mathematical formulation is:

Q, _ = qr((A A^H)^q A Omega, mode=’economic’)

where Omega is a Gaussian random matrix of shape (n, r + p).

When A is a low-rank matrix, the rangefinder provides a good approximation of the range of A with high probability.

Parameters:
  • A (LinearOperator) – The matrix to sketch.

  • r (int) – The target rank.

  • p (int, optional) – The number of over-sampling (default: 5).

  • q (int, optional) – Number of power iterations (default: 0).

  • seed (int, optional) – The seed for the random number generator (default: 1234).

  • **extra_args (dict) – Additional arguments. If ‘Omega’ is provided, use it as the sketching matrix.

Returns:

Q – Estimation of the range of A.

Return type:

ndarray

Modules

generalized_nystrom(X, r[, ...])

Generalized Nyström method

randomized_svd(X, r[, p, q, truncate, seed])

Randomized SVD algorithm.

rangefinder(A, r[, p, q, seed])

The (randomized) rangefinder method, also called HMT method.