low_rank_toolbox.matrices.quasi_svd

QuasiSVD low-rank matrix class and functions.

Authors: Benjamin Carrel and Rik Vorhaar

University of Geneva, 2022-2025

Classes

QuasiSVD(U, S, V, **extra_data)

Quasi Singular Value Decomposition (Quasi-SVD) for low-rank matrices.
class low_rank_toolbox.matrices.quasi_svd.QuasiSVD(U, S, V, **extra_data)[source]

Bases: LowRankMatrix

Quasi Singular Value Decomposition (Quasi-SVD) for low-rank matrices.

A generalization of the SVD where the middle matrix S is not necessarily diagonal. This representation is useful for operations that preserve orthogonality of U and V but may produce non-diagonal or even singular middle matrices.

Mathematical Representation

X = U @ S @ V.T

where: - U ∈ ℝ^(m×r) should have orthonormal columns (U^T @ U = I_r) - V ∈ ℝ^(n×q) should have orthonormal columns (V^T @ V = I_q) - S ∈ ℝ^(r×q) is a general matrix (not necessarily diagonal, may be singular)

Key Differences from SVD

  • SVD: S is diagonal with non-negative singular values, guaranteed non-singular

  • QuasiSVD: S is a general matrix, may be non-diagonal and potentially singular

  • QuasiSVD can represent intermediate results of matrix operations efficiently

  • Converting QuasiSVD → SVD requires O(r³) operations (SVD of S)

Storage Efficiency

Full matrix: O(mn) storage QuasiSVD: O(mr + rq + nq) storage Memory savings when r, q << min(m, n)

Important Notes

  • U and V are ASSUMED to have orthonormal columns (not verified at initialization)

  • Use .is_orthogonal() to verify orthonormality if needed

  • After matrix operations, orthogonality of U and V is preserved when possible, otherwise returns LowRankMatrix

  • After addition/subtraction, S may become singular or ill-conditioned

  • Use .truncate() to convert to SVD and remove small/zero singular values

  • Use .to_svd() to convert to diagonal form without truncation

U

Left matrix (assumed to have orthonormal columns)

Type:

ndarray, shape (m, r)

S

Middle matrix (general matrix, may be singular)

Type:

ndarray, shape (r, q)

V

Right matrix (assumed to have orthonormal columns)

Type:

ndarray, shape (n, q)

Vh

Hermitian conjugate of V (V.T.conj())

Type:

ndarray, shape (q, n)

Vt

Transpose of V (without conjugate)

Type:

ndarray, shape (q, n)

Ut

Transpose of U (without conjugate)

Type:

ndarray, shape (r, m)

Uh

Hermitian conjugate of U (U.T.conj())

Type:

ndarray, shape (r, m)

K

Product U @ S, computed on demand and not cached

Type:

ndarray, shape (m, r)

L

Product V @ S.T, computed on demand and not cached

Type:

ndarray, shape (r, n)

Properties
----------
shape

Shape of the represented matrix (m, n)

Type:

tuple

rank

Current rank of the representation (not necessarily matrix rank)

Type:

int

sing_vals[source]

Singular values of S (computed on demand and cached)

Type:

ndarray

Return type:

ndarray

Methods Overview
----------------
Core Operations

  • __add__, __sub__ : Addition/subtraction maintaining low-rank structure

  • __mul__ : Scalar or Hadamard (element-wise) multiplication

  • dot : Matrix-vector or matrix-matrix multiplication

  • hadamard : Element-wise multiplication with another matrix

Conversion & Truncation

  • to_svd() : Convert to SVD with diagonal S

  • truncate() : Convert to SVD and remove small singular values

Properties & Checks

  • is_symmetric() : Check if matrix is symmetric

  • is_orthogonal() : Check if U and V are orthonormal

  • is_singular() : Check if S is singular

  • norm() : Compute matrix norm (Frobenius, 2-norm, nuclear)

Class Methods

  • multi_add() : Efficient addition of multiple QuasiSVD matrices

  • multi_dot() : Efficient multiplication of multiple QuasiSVD matrices

  • generalized_nystroem() : Randomized low-rank approximation

Configuration
-------------
Default behavior controlled by
- DEFAULT_ATOL (machine precision)
Type:

Absolute tolerance for truncation

- DEFAULT_RTOL (None)
Type:

Relative tolerance for truncation

Examples

>>> import numpy as np
>>> from low_rank_toolbox.matrices import QuasiSVD
>>>
>>> # Create orthonormal matrices
>>> m, n, r = 100, 80, 10
>>> U, _ = np.linalg.qr(np.random.randn(m, r))
>>> V, _ = np.linalg.qr(np.random.randn(n, r))
>>> S = np.random.randn(r, r)
>>>
>>> # Create QuasiSVD
>>> X = QuasiSVD(U, S, V)
>>> print(X.shape)  # (100, 80)
>>> print(X.rank)   # 10
>>>
>>> # Operations preserve low-rank structure
>>> Y = X + X  # rank = 20 (sum of ranks)
>>> Z = X @ X.T  # matrix multiplication
>>>
>>> # Convert to SVD (diagonal S)
>>> X_svd = X.to_svd()
>>>
>>> # Truncate small singular values
>>> X_trunc = X.truncate(rtol=1e-10)
>>>
>>> # Addition of multiple matrices
>>> W = QuasiSVD.multi_add([X, -X])  # rank = 2*rank(X), call .truncate() to reduce

See also

SVD

Singular Value Decomposition with diagonal middle matrix

LowRankMatrix

Base class for low-rank matrix representations

References

Notes

  • This class is designed for numerical linear algebra applications

  • Operations may accumulate numerical errors; consider periodic truncation

  • Orthonormality of U and V is assumed but NOT enforced at initialization

  • S may become singular after operations like addition or subtraction

  • Use .is_orthogonal() to verify the orthonormality assumption

  • Use .is_singular() to check if S has become singular

  • For exact SVD with guaranteed diagonal non-singular S, use the SVD class

  • The representation is not unique: (U, S, V) and (UQ, Q^T S R, VR^T) represent the same matrix for any orthogonal Q, R

property H: QuasiSVD

Hermitian conjugate (conjugate transpose) of the QuasiSVD matrix.

For X = U @ S @ V.H, the Hermitian conjugate is:

X.H = (U @ S @ V.H).H = V @ S.H @ U.H = V @ S*.T @ U.H

Returns QuasiSVD(V, S.T.conj(), U) which stores [V, S*.T, U.H] and represents V @ S*.T @ U.H.

This is equivalent to X.T.conj() or X.conj().T.

Returns:

Hermitian conjugate in QuasiSVD format

Return type:

QuasiSVD

Notes

For real matrices, X.H is equivalent to X.T.

property K: ndarray

Compute and return the product U @ S on demand.

Returns:

Product U @ S, shape (m, r)

Return type:

ndarray

property L: ndarray

Compute and return the product V @ S.T on demand.

Returns:

Product V @ S.T, shape (n, r)

Return type:

ndarray

property S
property T: QuasiSVD

Transpose of the QuasiSVD matrix (without conjugation).

For X = U @ S @ V.H, the transpose is:

X.T = (U @ S @ V.H).T = V* @ S.T @ U.T

where V* denotes conjugate (not Hermitian).

Returns QuasiSVD(V.conj(), S.T, U.conj()) which stores [V*, S.T, U*.H] and represents V* @ S.T @ U.T.

Returns:

Transposed matrix in QuasiSVD format

Return type:

QuasiSVD

Notes

For real matrices, this is equivalent to QuasiSVD(V, S.T, U). For complex matrices, U and V must be conjugated.

property U
property Uh
property Ut
property V
property Vh
property Vt
__init__(U, S, V, **extra_data)[source]

Create a low-rank matrix stored by its SVD: Y = U @ S @ V.T

NOTE: U and V must be orthonormal NOTE: S is not necessarly diagonal and can be rectangular NOTE: The user must give V and not V.T or V.H

Parameters:

  • U (ndarray) – Left singular vectors, shape (m, r)

  • S (ndarray) – Non-singular matrix, shape (r, q)

  • V (ndarray) – Right singular vectors, shape (n, q)

Raises:

  • ValueError – If matrix dimensions do not match for multiplication.

  • TypeError – If S is provided as a 1D array or diagonal matrix (use SVD class instead).

Warning

MemoryEfficiencyWarning

If the low-rank representation uses more memory than dense storage.

cond_estimate()[source]

Condition number of the QuasiSVD matrix, estimated from the singular values.

Returns:

Condition number of the matrix.

Return type:

float

conj()[source]

Complex conjugate of the QuasiSVD matrix.

Returns QuasiSVD instead of generic LowRankMatrix.

Returns:

Complex conjugate in QuasiSVD format

Return type:

QuasiSVD

diag()[source]

Extract the diagonal of the matrix efficiently.

For QuasiSVD: diag[i] = U[i,:] @ S @ V[i,:].T This is O(mr²) instead of O(mn) for the full matrix.

Returns:

The diagonal elements

Return type:

ndarray

dot(other, side='right', dense_output=False)[source]

Matrix multiplication between SVD and other. The output is an SVD or a Matrix, depending on the type of other. If two QuasiSVD are multiplied, the new rank is the minimum of the two ranks.

If side is ‘right’ or ‘usual’, compute self @ other. If side is ‘left’ or ‘opposite’, compute other @ self.

Parameters:

  • other (QuasiSVD or ndarray) – Matrix to multiply

  • side (str, optional) – ‘left’ or ‘right’, by default ‘right’

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

Returns:

Result of the matrix multiplication

Return type:

QuasiSVD or ndarray

expm(inplace=False, **extra_data)[source]

Compute the matrix exponential exp(X) = e^X.

For QuasiSVD X = U @ S @ V.H, if X is Hermitian (V.H == U), then:

exp(X) = U @ exp(S) @ U.H

where exp(S) is the matrix exponential of S.

This is more efficient than general matrix exponentiation when the middle matrix S is small (r << n).

Parameters:

  • inplace (bool, optional) – If True, modify the current object in-place. Default is False.

  • **extra_data – Additional keyword arguments passed to __init__ for the new object.

Returns:

Matrix exponential in QuasiSVD format.

Return type:

QuasiSVD

Raises:

Notes

Computational cost: O(r³) for matrix exponential of r×r matrix S (vs O(n³) for general n×n matrix).

This method currently only supports Hermitian matrices where V.H == U. For general matrices, the eigenvalue decomposition would be needed.

Examples

>>> S = np.array([[1.0, 0.1], [0.1, 0.5]])
>>> U, _ = np.linalg.qr(np.random.randn(100, 2))
>>> X = QuasiSVD(U, S, U)  # Symmetric/Hermitian
>>> # Note: Current implementation has a bug - references self.s which doesn't exist
>>> # Use SVD.expm() for diagonal S matrices instead
>>> # X_exp = X.expm()  # Will raise AttributeError

See also

sqrtm

Matrix square root

SVD.expm

Optimized version for diagonal S

classmethod from_qr(qr)[source]

Convert QR format to QuasiSVD.

Parameters:

qr (QR) – QR representation to convert

Returns:

QuasiSVD representation

Return type:

QuasiSVD

Examples

>>> X_qr = QR(Q, R)
>>> X = QuasiSVD.from_qr(X_qr)
classmethod generate_random(shape, rank, seed=1234, is_symmetric=False, **extra_data)[source]

Generate a random QuasiSVD matrix with orthonormal U and V.

Parameters:

  • shape (tuple) – Shape of the matrix (m, n)

  • rank (int) – Rank of the matrix

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

  • is_symmetric (bool, optional) – If True, generate a symmetric matrix (only for square shapes), by default False

Returns:

Random QuasiSVD matrix

Return type:

QuasiSVD

Raises:

ValueError – If is_symmetric is True but shape is not square.

hadamard(other)[source]

Hadamard product between two QuasiSVD matrices

The new rank is the multiplication of the two ranks, at maximum. NOTE: if the expected rank is too large, dense matrices are used for the computation, but the output is still an SVD.

Parameters:

other (QuasiSVD, LowRankMatrix or ndarray) – Matrix to multiply

Returns:

Result of the Hadamard product.

Return type:

QuasiSVD or SVD or ndarray

is_orthogonal()[source]

Check if U and V have orthonormal columns.

Returns:

True if both U.H @ U ≈ I and V.H @ V ≈ I, False otherwise.

Return type:

bool

Notes

Result is cached after first computation. Uses np.allclose with default tolerances. Orthogonality is ASSUMED at initialization but not enforced. This method verifies the assumption.

is_singular()[source]

Check if middle matrix S is numerically singular.

Returns:

True if S is singular (condition number >= 1/machine_eps), False otherwise.

Return type:

bool

Notes

Result is cached after first computation. Uses condition number test: cond(S) >= 1/ε where ε is machine precision. S may become singular after operations like addition or subtraction.

is_symmetric()[source]

Check if the QuasiSVD matrix is symmetric.

Returns:

True if matrix is symmetric (square and U = V), False otherwise.

Return type:

bool

Notes

For QuasiSVD to be symmetric, U and V must be identical (within tolerance). This is a necessary condition when X = U @ S @ V.T. Uses np.allclose for comparison (tolerates small numerical errors).

lstsq(b, rtol=None, atol=np.float64(2.220446049250313e-14))[source]

Compute the least-squares solution to Xx ≈ b.

Minimizes ||Xx - b||₂ using the pseudoinverse: x = X⁺ @ b

Parameters:

  • b (ndarray) – Right-hand side vector or matrix. Shape (m,) or (m, k).

  • rtol (float, optional) – Relative tolerance for pseudoinverse computation. Default is None (uses machine precision * max(m,n)).

  • atol (float, optional) – Absolute tolerance for pseudoinverse. Default is DEFAULT_ATOL.

Returns:

Least-squares solution x. Shape (n,) or (n, k).

Return type:

ndarray

Examples

>>> X = QuasiSVD.generate_random((100, 80), 20)
>>> b = np.random.randn(100)
>>> x = X.lstsq(b)
>>> # x minimizes ||X @ x - b||
>>> residual = np.linalg.norm(X.dot(x) - b)

See also

pseudoinverse

Compute pseudoinverse

solve

Solve linear system

classmethod multi_add(matrices)[source]

Addition of multiple QuasiSVD matrices.

This is efficiently done by stacking the U, S, V matrices and then re-orthogonalizing. The resulting rank is at most the sum of all input ranks.

Parameters:

matrices (List[QuasiSVD]) – Matrices to add

Returns:

Sum of the matrices.

Return type:

QuasiSVD

Notes

  • Adding matrices increases rank: rank(X + Y) <= rank(X) + rank(Y)

  • X - X has rank 2*rank(X) but represents zero. Call .truncate() to reduce rank.

Examples

>>> Z = QuasiSVD.multi_add([X, -X])  # rank = 2*rank(X), but represents zero
>>> Z = Z.truncate(atol=1e-14)  # rank ~ 0
classmethod multi_dot(matrices)[source]

Matrix multiplication of several QuasiSVD matrices.

The rank of the output is the minimum of the ranks of the first and last inputs, at maximum. The matrices are multiplied all at once, so this method is more efficient than multiplying them one by one.

Parameters:

matrices (List[QuasiSVD]) – Matrices to multiply

Returns:

Product of the matrices.

Return type:

QuasiSVD

norm(ord='fro')[source]

Calculate matrix norm with caching and optimization for orthogonal case.

Parameters:

ord (str or int, default='fro') – Order of the norm. Supported values: - ‘fro’: Frobenius norm (default) - 2: Spectral norm (largest singular value) - ‘nuc’: Nuclear norm (sum of singular values) - Other values: computed via np.linalg.norm(self.full(), ord=ord)

Returns:

The computed norm value.

Return type:

float

Notes

Result is cached after first computation for each norm type. For orthogonal U and V, common norms are computed efficiently from singular values of S (O(r³)) instead of forming full matrix (O(mnr)). - Frobenius: ||X||_F = ||S||_F - Spectral (2-norm): ||X||_2 = σ_max(S) - Nuclear: ||X||_* = Σ σ_i(S)

norm_squared()[source]

Compute the squared Frobenius norm efficiently.

For orthogonal U and V: ||X||²_F = ||S||²_F This is O(r²) instead of O(mnr) for generic computation.

Returns:

The squared Frobenius norm

Return type:

float

numerical_health_check(verbose=True)[source]

Check the numerical health of the QuasiSVD representation.

This method checks: - Orthogonality of U and V - Condition number of S - Presence of near-zero singular values - Memory efficiency

Parameters:

verbose (bool, optional) – If True, print a summary. Default is True.

Returns:

Dictionary with health metrics: - ‘orthogonal_U’: bool, whether U is orthogonal - ‘orthogonal_V’: bool, whether V is orthogonal - ‘orthogonality_error_U’: float, ||U.H @ U - I||_F - ‘orthogonality_error_V’: float, ||V.H @ V - I||_F - ‘condition_number_S’: float, condition number of S - ‘is_singular’: bool, whether S is numerically singular - ‘min_singular_value’: float, smallest singular value of S - ‘max_singular_value’: float, largest singular value of S - ‘singular_value_ratio’: float, max/min singular values - ‘compression_ratio’: float, storage efficiency - ‘memory_efficient’: bool, whether low-rank is beneficial - ‘recommendations’: list of str, suggested actions

Return type:

dict

Examples

>>> X = QuasiSVD.generate_random((100, 80), 10)
>>> health = X.numerical_health_check()
>>> if not health['orthogonal_U']:
...     X = X.reorthogonalize()
property numerical_rank: int

Numerical rank of the QuasiSVD matrix.

The numerical rank is defined as the number of non-negligible singular values up to the machine precision.

Returns:

Numerical rank of the matrix.

Return type:

int

plot_singular_value_decay(semilogy=True, show=True, **kwargs)[source]

Plot the singular value decay of the middle matrix S.

This helps visualize the numerical rank and identify a good truncation threshold.

Parameters:

  • semilogy (bool, optional) – If True, use logarithmic scale for y-axis. Default is True.

  • show (bool, optional) – If True, call plt.show(). Default is True.

  • **kwargs – Additional keyword arguments passed to plt.plot()

Returns:

Matplotlib figure and axes objects

Return type:

fig, ax

Examples

>>> X = QuasiSVD.generate_random((100, 80), 20)
>>> fig, ax = X.plot_singular_value_decay()
>>> ax.axhline(1e-10, color='r', linestyle='--', label='Threshold')
>>> ax.legend()
project_onto_column_space(other, dense_output=False)[source]

Projection of other onto the column space of self.

The rank of the output is the rank of self, at maximum. Output is typically a QR object unless dense_output=True, in which case it is a dense ndarray.

The formula is given by:

P_U Y = UUh Y

where X = U S Vh is the SVD of matrix self and Y is the input (other) matrix to project.

Parameters:

  • other (ndarray or LowRankMatrix) – Matrix to project

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

Returns:

Projection of other onto the column space of self. Returns QR if dense_output=False, ndarray otherwise.

Return type:

QR or ndarray

project_onto_interpolated_tangent_space(mode='online', dense_output=False, **kwargs)[source]

Oblique projection onto the tangent space at self using interpolation (DEIM-like methods).

The formula is given by:

P_X Y = M_u Y_u - M_u Y_uv M_v^* + Y_v M_v^*

where M_u = U (P_U^* U)^{-1}, M_v = V (P_V^* V)^{-1}, Y_u = Y[p_u, :], Y_v = Y[:, p_v] and Y_uv = Y[p_u, p_v] with p_u and p_v being the interpolation indices for U and V, respectively.

Here, self = U S Vh is the SVD of matrix self and Y is the input matrix to project.

Two usage modes:

  1. Offline mode (provide pre-computed interpolatory matrices): Provide Y_u, Y_v, Y_uv, M_u, M_v in kwargs.

  2. Online mode (compute interpolation on-the-fly): Provide Y in kwargs. Optionally provide cssp_method_u and/or cssp_method_v. The indices and interpolatory matrices are computed using the specified CSSP methods. If not provided, QDEIM is used by default.

Parameters:

  • mode (str, optional) – Either ‘online’ or ‘offline’. Default is ‘online’.

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

  • **kwargs (dict) –

    For offline mode: Y_u, Y_v, Y_uv, M_u, M_v For online mode: Y (required), cssp_method_u (optional, default QDEIM),

    cssp_method_v (optional, default QDEIM), cssp_kwargs_u (optional), cssp_kwargs_v (optional)

Returns:

Oblique projection of Y onto the tangent space at self using interpolation.

Return type:

QuasiSVD

Examples

Offline mode:

>>> # Pre-compute interpolatory matrices
>>> p_u, M_u = DEIM(X.U, return_projector=True)
>>> p_v, M_v = DEIM(X.V, return_projector=True)
>>> Y_full = Y.full()
>>> Y_u = Y_full[p_u, :]
>>> Y_v = Y_full[:, p_v]
>>> Y_uv = Y_full[np.ix_(p_u, p_v)]
>>> result = X.project_onto_interpolated_tangent_space(
...     mode='offline',
...     Y_u=Y_u, Y_v=Y_v, Y_uv=Y_uv, M_u=M_u, M_v=M_v
... )

Online mode (with default QDEIM):

>>> result = X.project_onto_interpolated_tangent_space(
...     mode='online',
...     Y=Y
... )

Online mode (with custom CSSP methods):

>>> from low_rank_toolbox.cssp import DEIM, QDEIM
>>> result = X.project_onto_interpolated_tangent_space(
...     mode='online',
...     Y=Y,
...     cssp_method_u=DEIM,
...     cssp_method_v=QDEIM
... )
project_onto_row_space(other, dense_output=False)[source]

Projection of other onto the row space of self.

The rank of the output is the rank of self, at maximum. Output is typically a QR object unless dense_output=True, in which case it is a dense ndarray.

The formula is given by:

P_V Y = Y VVh

where X = U S Vh is the SVD of matrix self and Y is the input (other) matrix to project.

Parameters:

  • other (ndarray or LowRankMatrix) – Matrix to project

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

Returns:

Projection of other onto the row space of self. Returns QR if dense_output=False, ndarray otherwise.

Return type:

QR or ndarray

project_onto_tangent_space(other, dense_output=False)[source]

Projection of other onto the tangent space at self.

The rank of the output is two times the rank of self, at maximum.

The formula is given by:

P_X Y = UUh Y - UUh Y VVh + Y VVh

where X = U S Vh is the SVD of matrix self and Y is the input (other) matrix to project.

Parameters:

  • other (ndarray or LowRankMatrix) – Matrix to project

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

Returns:

Projection of other onto the tangent space at self.

Return type:

QuasiSVD

pseudoinverse(rtol=None, atol=np.float64(2.220446049250313e-14))[source]

Compute the Moore-Penrose pseudoinverse X⁺ efficiently.

For QuasiSVD X = U @ S @ V.T, the pseudoinverse is:

X⁺ = V @ S⁺ @ U.T

where S⁺ is the pseudoinverse of S computed via SVD.

Small singular values (< max(rtol*σ_max, atol)) are treated as zero.

Parameters:

  • rtol (float, optional) – Relative tolerance for determining zero singular values. Default is None (uses machine precision * max(m,n)).

  • atol (float, optional) – Absolute tolerance. Default is DEFAULT_ATOL (machine precision).

Returns:

Pseudoinverse in QuasiSVD format

Return type:

QuasiSVD

Examples

>>> X = QuasiSVD.generate_random((100, 80), 10)
>>> X_pinv = X.pseudoinverse()
>>> # Check: X @ X_pinv @ X ≈ X
>>> reconstruction = X.dot(X_pinv.dot(X.full(), dense_output=True), dense_output=True)
>>> np.allclose(X.full(), reconstruction)
True

See also

solve

Solve linear system Xx = b

lstsq

Least squares solution

rank_one_update(u, v, alpha=1.0)[source]

Efficient rank-1 update: X_new = X + alpha * u @ v.T

This adds a rank-1 matrix to the current QuasiSVD without forming the full matrix. The result has rank at most rank(X) + 1.

Parameters:

  • u (ndarray, shape (m,)) – Left vector for rank-1 update

  • v (ndarray, shape (n,)) – Right vector for rank-1 update

  • alpha (float, optional) – Scaling factor, default is 1.0

Returns:

Updated matrix with rank increased by at most 1

Return type:

QuasiSVD

Examples

>>> X = QuasiSVD.generate_random((100, 80), 5)
>>> u = np.random.randn(100)
>>> v = np.random.randn(80)
>>> X_new = X.rank_one_update(u, v, alpha=0.5)
>>> # X_new represents X + 0.5 * u @ v.T with rank at most 6
reorthogonalize(method='qr', inplace=False)[source]

Re-orthogonalize U and V if numerical drift has occurred.

After many operations, U and V may lose orthogonality due to numerical errors. This method restores orthogonality.

Parameters:

  • method (str, optional) – Method for re-orthogonalization: - ‘qr’: QR decomposition (default, stable and efficient) - ‘svd’: Full SVD (more expensive but more accurate)

  • inplace (bool, optional) – If True, modify the current object in-place. Otherwise, return a new object. Default is False.

Returns:

Re-orthogonalized QuasiSVD with same matrix values

Return type:

QuasiSVD

Examples

>>> X = QuasiSVD.generate_random((100, 80), 10)
>>> # ... many operations ...
>>> if not X.is_orthogonal():
...     X = X.reorthogonalize()
sing_vals()[source]

Singular values of the QuasiSVD matrix. Computed from S, then cached.

Returns:

Singular values of the middle matrix S in descending order. Shape is (min(r, q),) where S has shape (r, q).

Return type:

ndarray

Notes

Result is cached after first computation for efficiency. Computes via np.linalg.svdvals(S), which is O(r²q) for S with shape (r, q).

solve(b, method='lstsq')[source]

Solve the linear system Xx = b.

For square full-rank matrices, solves Xx = b. For rectangular or rank-deficient matrices, computes the least-squares solution.

Parameters:

  • b (ndarray) – Right-hand side vector or matrix. Shape (m,) or (m, k).

  • method (str, optional) –

    Solution method: - ‘lstsq’: Use pseudoinverse (default, works for any shape) - ‘direct’: Use factored form assuming orthogonality (faster but requires

    orthogonal U, V and invertible S)

Returns:

Solution x. Shape (n,) or (n, k).

Return type:

ndarray

Raises:

ValueError – If matrix is not square and method=’direct’. If dimensions are incompatible.

Examples

>>> X = QuasiSVD.generate_random((100, 100), 100)  # Full rank for unique solution
>>> b = np.random.randn(100)
>>> x = X.solve(b)
>>> # Check: X @ x ≈ b
>>> np.allclose(X.dot(x), b)
True

>>> # For rank-deficient systems, use lstsq instead:
>>> X_deficient = QuasiSVD.generate_random((100, 100), 20)
>>> x_ls = X_deficient.lstsq(b)

See also

pseudoinverse

Compute pseudoinverse

lstsq

Least squares solution

Notes

The default method is ‘lstsq’ which uses the pseudoinverse and is robust to non-orthogonal factors and rank deficiency. Use method=’direct’ only when you are certain that U and V are orthogonal and S is invertible.

sqrtm(inplace=False, **extra_data)[source]

Compute the matrix square root X^{1/2} such that X^{1/2} @ X^{1/2} = X.

For QuasiSVD X = U @ S @ V.T, the square root is:

X^{1/2} = U @ S^{1/2} @ V.T

where S^{1/2} is the matrix square root of S.

Parameters:

inplace (bool, optional) – If True, modify the current object in-place. Default is False.

Returns:

Matrix square root in QuasiSVD format

Return type:

QuasiSVD

property svd_type: str

Classify the type of SVD representation based on S dimensions.

Returns:

One of: - ‘full’: S has shape (m, n) - full SVD with all singular vectors - ‘reduced’: S has shape (min(m,n), min(m,n)) - reduced/economic SVD - ‘truncated’: S has shape (r, r) where r < min(m,n) - truncated SVD - ‘unconventional’: S has shape (r, k) where r ≠ k - general quasi-SVD

Return type:

str

Examples

>>> X = QuasiSVD.generate_random((100, 80), 50)
>>> print(X.svd_type)  # 'reduced' (50 = min(100,80))
>>> X_trunc = X.truncate(r=10)
>>> print(X_trunc.svd_type)  # 'truncated' (10 < min(100,80))
to_qr()[source]

Convert QuasiSVD to QR format.

Returns X = Q @ R where Q is orthogonal and R is upper triangular. This is useful for efficient column-space operations.

Returns:

QR representation of the matrix

Return type:

QR

Examples

>>> X = QuasiSVD.generate_random((100, 80), 10)
>>> X_qr = X.to_qr()
to_svd()[source]

Convert QuasiSVD to SVD by computing the SVD of the middle matrix S.

This operation computes U_s, s, V_s = svd(S) and returns:

SVD(self.U @ U_s, s, self.V @ V_s)

If you want to truncate small singular values, use the .truncate() method instead.

Returns:

An SVD object with diagonal S matrix

Return type:

SVD

Notes

This method imports SVD only at runtime to avoid circular imports. The computational cost is O(r³) where r is the rank of S.

trace()[source]

Compute the trace of the matrix efficiently.

For square QuasiSVD: trace(X) = trace(U @ S @ V.T) = trace(V.T @ U @ S) This is O(r³) instead of O(mn) for the full matrix.

For non-square matrices, trace is only defined on the square part.

Returns:

The trace of the matrix

Return type:

float

Raises:

ValueError – If matrix is not square

transpose()[source]

Transpose the matrix (returns QuasiSVD).

Return type:

QuasiSVD

truncate(r=None, rtol=None, atol=np.float64(2.220446049250313e-14), inplace=False)[source]

Truncate the QuasiSVD by converting to SVD and removing small singular values.

The QuasiSVD is first converted to SVD (diagonal S), then singular values are truncated based on rank or tolerance criteria.

Parameters:

  • r (int, optional) – Target rank. If specified, keep only the r largest singular values.

  • rtol (float, optional) – Relative tolerance. Singular values < rtol * max(singular_values) are removed.

  • atol (float, optional) – Absolute tolerance. Singular values < atol are removed. Default is DEFAULT_ATOL (machine precision).

  • inplace (bool, optional) – If True, this parameter is ignored (QuasiSVD cannot be truncated in-place, must convert to SVD). Kept for API compatibility. Default is False.

Returns:

Truncated SVD object

Return type:

SVD

Notes

Priority of truncation criteria (from highest to lowest): 1. r (explicit rank) 2. rtol (relative tolerance) 3. atol (absolute tolerance)

Examples

>>> X = QuasiSVD(U, S, V)
>>> X_trunc = X.truncate(r=10)  # Keep top 10 singular values
>>> X_trunc = X.truncate(rtol=1e-6)  # Keep s_i > 1e-6 * s_max
>>> X_trunc = X.truncate(atol=1e-10)  # Keep s_i > 1e-10
Parameters:

  • U (ndarray)

  • S (ndarray)

  • V (ndarray)