low_rank_toolbox.matrices.LowRankMatrix

class low_rank_toolbox.matrices.LowRankMatrix(*matrices, **extra_data)

Bases: LinearOperator

Meta 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.)

__init__(*matrices, **extra_data)

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.

Methods

__init__(*matrices, **extra_data)	Initialize a low-rank matrix from a sequence of factor matrices.
adjoint()	Hermitian adjoint.
allclose(other[, rtol, atol])	Check approximate equality with another matrix.
approximation_error(reference[, ord])	Compute approximation error ||X - reference|| efficiently.
compress()	Compress the factorization to maximize memory efficiency.
compress_()	Compress the factorization in-place to maximize memory efficiency.
compression_ratio()	Compute compression ratio relative to dense storage.
cond_estimate([method, n_iter])	Estimate condition number without computing full SVD.
conj()	Complex conjugate of the matrix.
copy()	Create a deep copy of the matrix.
create_data_alias(key)	Create a property that provides access to an entry in extra_data.
create_matrix_alias(index[, transpose, ...])	Create a property that provides access to a factor matrix with optional transformations.
diag()	Extract diagonal without forming full matrix.
dot(other[, side, dense_output])	Matrix and vector multiplication
dot_sparse(sparse_other[, side, dense_output])	Efficient multiplication with a sparse matrix.
expm_multiply(A, h[, side, dense_output])	Efficient action of sparse matrix exponential left: output = exp(h*A) @ self right: output = self @ exp(h*A)
flatten()	Flatten the matrix to a 1D array.
from_dense(matrix)	Alias for from_matrix().
from_full(matrix)	Alias for from_matrix().
from_low_rank(low_rank_matrix)	Convert a LowRankMatrix to this specific subclass format.
from_matrix(matrix)	Create a low-rank matrix from a full matrix.
full()	Form the full dense matrix by multiplying all factors in optimal order.
gather(indices)	Access specific matrix entries without forming the full matrix.
get_block(rows, cols)	Extract block submatrix.
hadamard(other)	Element-wise (Hadamard) product with another matrix.
is_symmetric()	Check if the matrix is symmetric.
is_well_conditioned([threshold])	Check if the factorization is numerically well-conditioned.
load(filename)	Load low-rank matrix from disk.
matmat(X)	Matrix-matrix multiplication.
matvec(v)	Matrix-vector product (public interface, backward compatibility).
memory_usage([unit])	Report actual memory used by the factorization.
multi_add(others)	Add multiple low-rank matrices efficiently.
multi_dot(others)	Matrix multiplication of a sequence of matrices.
norm([ord])	Default implementation, overload this for some subclasses
norm_squared()	Compute squared Frobenius norm efficiently: ||X||²_F = tr(X^H X).
power(n)	Compute matrix power X^n efficiently using repeated squaring.
rmatmat(X)	Adjoint matrix-matrix multiplication.
rmatvec(v)	Adjoint matrix-vector product (public interface, backward compatibility).
save(filename)	Save low-rank matrix to disk efficiently.
scale_(scalar)	Scale the matrix in-place by a scalar.
to_dense()	Convert to dense matrix (alias for full).
to_full()	Convert to dense matrix (alias for full).
to_sparse([format, threshold])	Convert to sparse format, zeroing small entries.
todense()	Convert to dense matrix (alias for full).
trace()	Compute trace efficiently using cyclic property: tr(ABC) = tr(CAB) = tr(BCA).
transpose()	Transpose the matrix (alias for .T property).

Attributes

H	Hermitian transpose (conjugate transpose) of the matrix.
T	Transpose of the matrix (reverse order and transpose each factor).
deepshape	Shape tuple including all intermediate dimensions of the factorization.
is_memory_efficient	Check if the low-rank representation uses less memory than dense storage.
length	Number of factor matrices in the factorization.
ndim	Number of dimensions (always 2 for matrices).
rank	Rank of the low-rank factorization.
size	Total number of elements stored in the factorization.

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)

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)

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:

bool

Notes

This forms the full matrix for comparison, which may be inefficient.

approximation_error(reference, ord='fro')

Compute approximation error ||X - reference|| efficiently.

Parameters:

• reference (ndarray) – Reference matrix to compare against.

• ord (str or int, optional) – Norm type, by default ‘fro’ (Frobenius).

Returns:

Approximation error in the specified norm.

Return type:

float

Notes

This forms the full matrix, which may be inefficient for very large matrices. For Frobenius norm, consider using norm_squared() for efficiency.

compress()

Compress the factorization to maximize memory efficiency.

Return type:

LowRankMatrix

compress_()

Compress the factorization in-place to maximize memory efficiency.

Returns:

Self (modified in-place).

Return type:

LowRankMatrix

Notes

This modifies the matrix in-place. Use compress() for a non-destructive version.

compression_ratio()

Compute compression ratio relative to dense storage.

Returns:

Ratio of low-rank storage to dense storage (< 1 means savings).

Return type:

float

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)

Estimate condition number without computing full SVD.

Parameters:

• method (str, optional) – Method to use: ‘power_iteration’ (default) or ‘norm_ratio’.

• n_iter (int, optional) – Number of iterations for power iteration method, by default 10.

Returns:

Estimated condition number (ratio of largest to smallest singular value).

Return type:

float

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.

conj()

Complex conjugate of the matrix.

Return type:

LowRankMatrix

copy()

Create a deep copy of the matrix.

Returns:

A deep copy of this matrix with independent factor matrices.

Return type:

LowRankMatrix

static create_data_alias(key)

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:

property

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)

Create a property that provides access to a factor matrix with optional transformations.

Parameters:

• index (int) – Index of the factor matrix in self._matrices.

• transpose (bool, optional) – Whether to transpose the matrix when accessed, by default False.

• conjugate (bool, optional) – Whether to conjugate the matrix when accessed, by default False.

Returns:

Property object for accessing the transformed matrix.

Return type:

property

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:

tuple

diag()

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)

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:

ndarray | LowRankMatrix

dot_sparse(sparse_other, side='usual', dense_output=False)

Efficient multiplication with a sparse matrix.

Parameters:

• sparse_other (spmatrix) – Sparse matrix to multiply with.

• side (str, optional) – ‘right’ or ‘usual’: output = self @ sparse_other ‘left’ or ‘opposite’: output = sparse_other @ self

• dense_output (bool, optional) – Whether to return a dense matrix, by default False.

Returns:

Result of the multiplication.

Return type:

ndarray | LowRankMatrix

expm_multiply(A, h, side='left', dense_output=False, **extra_args)

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

flatten()

Flatten the matrix to a 1D array.

Return type:

ndarray

classmethod from_dense(matrix)

Alias for from_matrix(). Must be implemented by subclasses.

Parameters:

matrix (ndarray)

classmethod from_full(matrix)

Alias for from_matrix(). Must be implemented by subclasses.

Parameters:

matrix (ndarray)

classmethod from_low_rank(low_rank_matrix)

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:

LowRankMatrix

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)

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:

LowRankMatrix

Parameters:

matrix (ndarray)

full()

Form the full dense matrix by multiplying all factors in optimal order.

Return type:

ndarray

gather(indices)

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)

Extract block submatrix.

Parameters:

• rows (slice) – Row slice specification.

• cols (slice) – Column slice specification.

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)

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:

bool

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()

Check if the matrix is symmetric.

Returns:

True if the matrix is square and symmetric, False otherwise.

Return type:

bool

Notes

This method forms the full matrix to check symmetry. For large matrices, this may be memory-intensive.

is_well_conditioned(threshold=10000000000.0)

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:

bool

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:

int

classmethod load(filename)

Load low-rank matrix from disk.

Parameters:

filename (str) – Path to the saved matrix file.

Returns:

Loaded low-rank matrix.

Return type:

LowRankMatrix

matvec(v)

Matrix-vector product (public interface, backward compatibility).

Parameters:

v (ndarray) – Vector to multiply with.

Returns:

Result of matrix-vector multiplication.

Return type:

ndarray

memory_usage(unit='MB')

Report actual memory used by the factorization.

Parameters:

unit (str, optional) – Unit for memory size: ‘B’, ‘KB’, ‘MB’, ‘GB’, by default ‘MB’.

Returns:

Memory usage in the specified unit.

Return type:

float

multi_add(others)

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)

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:

LowRankMatrix

property ndim: int

Number of dimensions (always 2 for matrices).

norm(ord='fro')

Default implementation, overload this for some subclasses

Return type:

float

Parameters:

ord (str | int)

norm_squared()

Compute squared Frobenius norm efficiently: ||X||²_F = tr(X^H X).

Returns:

Squared Frobenius norm of the matrix.

Return type:

float

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)

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:

int

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)

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)

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)

Scale the matrix in-place by a scalar.

Parameters:

scalar (float) – Scalar to multiply the matrix by.

Returns:

Self (modified in-place).

Return type:

LowRankMatrix

property size: int

Total number of elements stored in the factorization.

Returns:

Sum of elements across all factor matrices.

Return type:

int

Notes

This returns the storage cost, not the matrix size (m×n). For matrix dimensions, use .shape. For compression ratio, use .compression_ratio().

to_dense()

Convert to dense matrix (alias for full).

Return type:

ndarray

to_full()

Convert to dense matrix (alias for full).

Return type:

ndarray

to_sparse(format='csr', threshold=1e-10)

Convert to sparse format, zeroing small entries.

Parameters:

• format (str, optional) – Sparse matrix format: ‘csr’, ‘csc’, ‘coo’, by default ‘csr’.

• threshold (float, optional) – Entries with absolute value below threshold are set to zero, by default 1e-10.

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.

todense()

Convert to dense matrix (alias for full).

Return type:

ndarray

trace()

Compute trace efficiently using cyclic property: tr(ABC) = tr(CAB) = tr(BCA).

Returns:

Trace of the matrix.

Return type:

float

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.

transpose()

Transpose the matrix (alias for .T property).

Return type:

LowRankMatrix

Parameters:

matrices (ndarray)