low_rank_toolbox.matrices.low_rank_matrix
Generic low-rank matrix base class and utilities.
- Authors: Benjamin Carrel and Rik Vorhaar
University of Geneva, 2022-2025
Classes
|
Meta class for dealing with low rank matrices in different formats. |
Exceptions
Warning for inefficient operations on low-rank matrices. |
|
Warning when low-rank representation uses more memory than dense storage. |
- exception low_rank_toolbox.matrices.low_rank_matrix.LowRankEfficiencyWarning[source]
Bases:
WarningWarning for inefficient operations on low-rank matrices.
- class low_rank_toolbox.matrices.low_rank_matrix.LowRankMatrix(*matrices, **extra_data)[source]
Bases:
LinearOperatorMeta class for dealing with low rank matrices in different formats.
Inherits from scipy.sparse.linalg.LinearOperator for seamless integration with iterative solvers and standard numerical linear algebra tools.
Do not use this class directly, but rather use its subclasses.
We always decompose a matrix as a product of smaller matrices. These smaller matrices are stored in
self._matrices.LinearOperator Integration
As a subclass of scipy.sparse.linalg.LinearOperator, LowRankMatrix provides efficient matrix-vector and matrix-matrix products without forming the full dense matrix. This enables direct use with scipy’s iterative solvers and other numerical algorithms.
Key LinearOperator features: - Matrix-vector multiplication: A @ v or A.matvec(v) - Adjoint multiplication: A.H @ v or A.rmatvec(v) - Matrix-matrix multiplication: A @ B or A.matmat(B) - Lazy composition: A + B, A @ B (when B is LinearOperator) - Shape and dtype attributes: A.shape, A.dtype
Examples
Using with scipy iterative solvers:
>>> from scipy.sparse.linalg import gmres, LinearOperator >>> from low_rank_toolbox import SVD >>> import numpy as np >>> >>> # Create a well-conditioned low-rank matrix >>> U, _ = np.linalg.qr(np.random.randn(1000, 10)) >>> s = np.logspace(0, -1, 10) # Better conditioned >>> V, _ = np.linalg.qr(np.random.randn(1000, 10)) >>> # Add diagonal regularization for better conditioning >>> A = SVD(U, s, V) >>> A_reg = A + 0.1 * LinearOperator((1000, 1000), matvec=lambda x: x) >>> >>> # Solve Ax = b using GMRES (never forms the full matrix) >>> b = np.random.randn(1000) >>> x, info = gmres(A_reg, b, rtol=1e-6, atol=1e-6, maxiter=100) >>> assert info == 0 # Verify convergence
Lazy composition with other operators:
>>> from scipy.sparse import diags >>> from scipy.sparse.linalg import aslinearoperator >>> >>> # Create a diagonal operator for better conditioning >>> D = diags([0.5 for i in range(1000)]) >>> D_op = aslinearoperator(D) >>> >>> # Lazy sum - doesn't form full matrix >>> B = A + D_op # Returns _SumLinearOperator >>> >>> # Use in iterative solver >>> x2, info2 = gmres(B, b, rtol=1e-6, atol=1e-6, maxiter=100) >>> assert info2 == 0 # Verify convergence
Matrix-vector products (efficient, no dense formation):
>>> v = np.random.randn(1000) >>> y = A @ v # Efficient: never forms full 1000x1000 matrix >>> >>> # Adjoint product >>> z = A.H @ v # Hermitian transpose product
Custom preconditioners:
>>> # Use the regularized matrix from above >>> def precondition(v): ... # Simple diagonal preconditioner ... return v >>> >>> # Create LinearOperator from function >>> M = LinearOperator((1000, 1000), matvec=precondition) >>> >>> # Use as preconditioner >>> x3, info3 = gmres(A_reg, b, M=M, rtol=1e-6, atol=1e-6, maxiter=100) >>> assert info3 == 0 # Verify convergence
Notes
All matrix-vector products are computed efficiently in O(rank) operations
Full matrix is never formed unless explicitly requested with .full()
Subclasses (SVD, QR, QuasiSVD) inherit all LinearOperator functionality
Compatible with all scipy.sparse.linalg iterative solvers (gmres, cg, bicgstab, etc.)
- property H
Hermitian transpose (conjugate transpose) of the matrix.
- property T: LowRankMatrix
Transpose of the matrix (reverse order and transpose each factor).
- __init__(*matrices, **extra_data)[source]
Initialize a low-rank matrix from a sequence of factor matrices.
- Parameters:
*matrices (ndarray) – Sequence of matrices whose product forms the low-rank matrix. At least two matrices must be provided, and their shapes must align for matrix multiplication (i.e., matrices[i].shape[1] == matrices[i+1].shape[0]).
**extra_data (dict) – Additional data to store with the matrix (e.g., poles, residues).
- Raises:
ValueError – If fewer than 2 matrices are provided or if matrix shapes do not align.
- allclose(other, rtol=1e-05, atol=1e-08)[source]
Check approximate equality with another matrix.
- Parameters:
other (LowRankMatrix or ndarray) – Matrix to compare with.
rtol (float, optional) – Relative tolerance, by default 1e-5.
atol (float, optional) – Absolute tolerance, by default 1e-8.
- Returns:
True if matrices are approximately equal within tolerances.
- Return type:
Notes
This forms the full matrix for comparison, which may be inefficient.
- approximation_error(reference, ord='fro')[source]
Compute approximation error ||X - reference|| efficiently.
- Parameters:
- Returns:
Approximation error in the specified norm.
- Return type:
Notes
This forms the full matrix, which may be inefficient for very large matrices. For Frobenius norm, consider using norm_squared() for efficiency.
- compress_()[source]
Compress the factorization in-place to maximize memory efficiency.
- Returns:
Self (modified in-place).
- Return type:
Notes
This modifies the matrix in-place. Use compress() for a non-destructive version.
- compression_ratio()[source]
Compute compression ratio relative to dense storage.
- Returns:
Ratio of low-rank storage to dense storage (< 1 means savings).
- Return type:
Notes
A ratio of 0.1 means the low-rank format uses 10% of the memory of the dense format, i.e., a 10x compression.
- cond_estimate(method='power_iteration', n_iter=10)[source]
Estimate condition number without computing full SVD.
- Parameters:
- Returns:
Estimated condition number (ratio of largest to smallest singular value).
- Return type:
Notes
The ‘power_iteration’ method estimates the largest singular value using power iteration. The smallest is estimated using the inverse (via solving linear systems). This is approximate but avoids full SVD.
The ‘norm_ratio’ method uses matrix norms as bounds.
- copy()[source]
Create a deep copy of the matrix.
- Returns:
A deep copy of this matrix with independent factor matrices.
- Return type:
- static create_data_alias(key)[source]
Create a property that provides access to an entry in extra_data.
- Parameters:
key (str) – Key for accessing the data in self._extra_data.
- Returns:
Property object for accessing the extra data.
- Return type:
Examples
Create a property to access extra data stored with the matrix:
>>> class RationalKrylov(LowRankMatrix): ... poles = LowRankMatrix.create_data_alias('poles') ... residues = LowRankMatrix.create_data_alias('residues')
Now my_matrix.poles accesses self._extra_data[‘poles’].
- static create_matrix_alias(index, transpose=False, conjugate=False)[source]
Create a property that provides access to a factor matrix with optional transformations.
- Parameters:
- Returns:
Property object for accessing the transformed matrix.
- Return type:
Examples
Create a property to access the first factor matrix:
>>> class MyLowRank(LowRankMatrix): ... U = LowRankMatrix.create_matrix_alias(0) ... V = LowRankMatrix.create_matrix_alias(1, transpose=True)
Now my_matrix.U accesses self._matrices[0] and my_matrix.V accesses self._matrices[1].T.
- property deepshape: tuple
Shape tuple including all intermediate dimensions of the factorization.
- Returns:
Tuple of dimensions showing the shape of each factor in the chain. For factors A (m×k), B (k×n), returns (m, k, n).
- Return type:
- diag()[source]
Extract diagonal without forming full matrix.
- Returns:
Diagonal elements of the matrix.
- Return type:
ndarray
Notes
This is computed efficiently by extracting only the diagonal elements during the matrix chain multiplication, avoiding forming the full matrix.
- dot(other, side='right', dense_output=False)[source]
Matrix and vector multiplication
- Parameters:
other (LowRankMatrix, ndarray, spmatrix) – Matrix or vector to multiply with.
side (str, optional) – Whether to multiply on the left or right, by default ‘right’.
dense_output (bool, optional) – Whether to return a dense matrix or a low-rank matrix, by default False.
- Return type:
- dot_sparse(sparse_other, side='usual', dense_output=False)[source]
Efficient multiplication with a sparse matrix.
- Parameters:
- Returns:
Result of the multiplication.
- Return type:
ndarray | LowRankMatrix
- expm_multiply(A, h, side='left', dense_output=False, **extra_args)[source]
Efficient action of sparse matrix exponential left: output = exp(h*A) @ self right: output = self @ exp(h*A)
- Parameters:
A (spmatrix) – Sparse matrix in the exponential.
h (float) – Time step.
side (str, optional) – Whether to multiply on the left or right, by default ‘left’.
dense_output (bool, optional) – Whether to return a dense matrix or a low-rank matrix, by default False.
extra_args (dict, optional) – Extra arguments to pass to scipy.sparse.linalg.expm_multiply.
- Returns:
Resulting matrix after applying the exponential action.
- Return type:
ndarray | LowRankMatrix
- classmethod from_dense(matrix)[source]
Alias for from_matrix(). Must be implemented by subclasses.
- Parameters:
matrix (ndarray)
- classmethod from_full(matrix)[source]
Alias for from_matrix(). Must be implemented by subclasses.
- Parameters:
matrix (ndarray)
- classmethod from_low_rank(low_rank_matrix)[source]
Convert a LowRankMatrix to this specific subclass format.
- Parameters:
low_rank_matrix (LowRankMatrix) – Existing low-rank matrix to convert.
- Returns:
New matrix of the target subclass type.
- Return type:
Notes
Subclasses should override this to perform format-specific conversions. For example, SVD.from_low_rank() would compute the SVD of the input. The base class implementation creates a generic LowRankMatrix.
Examples
>>> X_generic = LowRankMatrix(A, B) >>> X_svd = SVD.from_low_rank(X_generic) # Converts to SVD format
- classmethod from_matrix(matrix)[source]
Create a low-rank matrix from a full matrix.
This method must be implemented by subclasses (e.g., SVD, QR). The base LowRankMatrix class requires at least 2 matrices for factorization.
- Return type:
- Parameters:
matrix (ndarray)
- full()[source]
Form the full dense matrix by multiplying all factors in optimal order.
- Return type:
- gather(indices)[source]
Access specific matrix entries without forming the full matrix.
- Parameters:
indices (ndarray or list) – Index specification. For single element: [row_idx, col_idx]. For multiple elements (fancy indexing): [row_indices, col_indices].
- Returns:
The requested matrix element(s).
- Return type:
float | ndarray
Notes
This is faster and more memory-efficient than forming the full matrix. Very useful for matrix completion tasks or estimating reconstruction error.
- get_block(rows, cols)[source]
Extract block submatrix.
- Parameters:
- Returns:
The requested block submatrix.
- Return type:
ndarray
Notes
This method forms the full matrix, which may be inefficient. It’s provided for convenience and compatibility.
- hadamard(other)[source]
Element-wise (Hadamard) product with another matrix.
- Parameters:
other (LowRankMatrix or ndarray) – Matrix to multiply element-wise.
- Returns:
Dense result of element-wise multiplication.
- Return type:
ndarray
- property is_memory_efficient: bool
Check if the low-rank representation uses less memory than dense storage.
- Returns:
True if low-rank format uses less memory than dense format, False otherwise.
- Return type:
Notes
This compares the total number of elements stored in the factorization versus the full m×n matrix. For empty matrices (zero size), returns True.
Examples
>>> A = LowRankMatrix(np.random.randn(1000, 10), np.random.randn(10, 1000)) >>> A.is_memory_efficient # True: 20,000 elements vs 1,000,000 True >>> B = LowRankMatrix(np.random.randn(100, 90), np.random.randn(90, 100)) >>> B.is_memory_efficient # False: 18,000 elements vs 10,000 False
- is_symmetric()[source]
Check if the matrix is symmetric.
- Returns:
True if the matrix is square and symmetric, False otherwise.
- Return type:
Notes
This method forms the full matrix to check symmetry. For large matrices, this may be memory-intensive.
- is_well_conditioned(threshold=10000000000.0)[source]
Check if the factorization is numerically well-conditioned.
- Parameters:
threshold (float, optional) – Condition number threshold, by default 1e10.
- Returns:
True if estimated condition number is below threshold.
- Return type:
Notes
This uses an approximate condition number estimate and may not be accurate for all cases. Consider this a heuristic check.
- property length: int
Number of factor matrices in the factorization.
- Returns:
The number of matrices in the product chain.
- Return type:
- classmethod load(filename)[source]
Load low-rank matrix from disk.
- Parameters:
filename (str) – Path to the saved matrix file.
- Returns:
Loaded low-rank matrix.
- Return type:
- matvec(v)[source]
Matrix-vector product (public interface, backward compatibility).
- Parameters:
v (ndarray) – Vector to multiply with.
- Returns:
Result of matrix-vector multiplication.
- Return type:
ndarray
- multi_add(others)[source]
Add multiple low-rank matrices efficiently.
- Parameters:
others (Sequence[LowRankMatrix]) – Sequence of low-rank matrices to add.
- Returns:
Dense matrix representing the sum.
- Return type:
ndarray
Notes
While more efficient than repeated use of +, this still requires forming full dense matrices and may be memory-intensive.
- multi_dot(others)[source]
Matrix multiplication of a sequence of matrices.
- Parameters:
others (Sequence[LowRankMatrix | ndarray]) – Sequence of matrices to multiply.
- Returns:
Low rank matrix representing the product.
- Return type:
- norm_squared()[source]
Compute squared Frobenius norm efficiently: ||X||²_F = tr(X^H X).
- Returns:
Squared Frobenius norm of the matrix.
- Return type:
Notes
This is more efficient than computing the norm and squaring it, as it avoids the square root operation. For complex matrices, uses Hermitian transpose (X^H X) to ensure real result.
- power(n)[source]
Compute matrix power X^n efficiently using repeated squaring.
- Parameters:
n (int) – Power to raise the matrix to. Must be non-negative.
- Returns:
Matrix raised to the n-th power. Returns identity matrix (as ndarray) for n=0, self for n=1, and LowRankMatrix for n>1.
- Return type:
LowRankMatrix | ndarray
- Raises:
ValueError – If matrix is not square or n is negative.
Notes
This uses the binary exponentiation algorithm for efficiency. Time complexity is O(log n) matrix multiplications.
- property rank: int
Rank of the low-rank factorization.
- Returns:
The maximum possible rank given the factorization structure, computed as the minimum dimension across all factor matrices.
- Return type:
Notes
This is an upper bound on the true numerical rank. The actual rank may be lower if the factors are rank-deficient. For the true numerical rank, compute SVD of the full matrix.
- rmatvec(v)[source]
Adjoint matrix-vector product (public interface, backward compatibility).
- Parameters:
v (ndarray) – Vector to multiply with.
- Returns:
Result of adjoint matrix-vector multiplication.
- Return type:
ndarray
- save(filename)[source]
Save low-rank matrix to disk efficiently.
- Parameters:
filename (str) – Path to save the matrix. Extension ‘.npz’ will be added if not present.
Notes
The matrix is saved in compressed NumPy format with all factor matrices and extra data preserved.
- scale_(scalar)[source]
Scale the matrix in-place by a scalar.
- Parameters:
scalar (float) – Scalar to multiply the matrix by.
- Returns:
Self (modified in-place).
- Return type:
- property size: int
Total number of elements stored in the factorization.
- Returns:
Sum of elements across all factor matrices.
- Return type:
Notes
This returns the storage cost, not the matrix size (m×n). For matrix dimensions, use .shape. For compression ratio, use .compression_ratio().
- to_sparse(format='csr', threshold=1e-10)[source]
Convert to sparse format, zeroing small entries.
- Parameters:
- Returns:
Sparse matrix representation.
- Return type:
spmatrix
Notes
This forms the full dense matrix first, which may be inefficient. Only use this if you expect the resulting matrix to be sparse.
- trace()[source]
Compute trace efficiently using cyclic property: tr(ABC) = tr(CAB) = tr(BCA).
- Returns:
Trace of the matrix.
- Return type:
Notes
For a product of matrices A₁A₂…Aₙ, this method uses the cyclic property of the trace to minimize computational cost. The trace is computed as tr(A₂…AₙA₁) by cyclically permuting to minimize the intermediate matrix sizes.
- Parameters:
matrices (ndarray)