truncated_svd#

t3toolbox.backend.linalg.truncated_svd(A, min_rank=None, max_rank=None, rtol=None, atol=None)#
def truncated_svd(
        A: NDArray, # shape=(...,N,M)
        min_rank: int = None,  # 1 <= min_rank <= max_rank <= minimum(N, M)
        max_rank: int = None,  # 1 <= min_rank <= max_rank <= minimum(N, M)
        rtol: float = None,  # removes singular values satisfying sigma < maximum(atol, rtol*sigma1)
        atol: float = None,  # removes singular values satisfying sigma < maximum(atol, rtol*sigma1)
) -> typ.Tuple[
    NDArray, # U, shape=(...,N,r)
    NDArray, # ss, shape=(...,r)
    NDArray, # Vt, shape=(...,r,M)
]:

Compute (truncated) singular value decomposition of matrix A.

A = U @ diag(ss) @ Vt Equality may be approximate if truncation is used.

Parameters:

  • A (NDArray) – Matrix. shape=(…, N, M)

  • min_rank (int) – Minimum rank for truncation. Should have 1 <= min_rank <= max_rank <= minimum(N, M).

  • min_rank – Maximum rank for truncation. Should have 1 <= min_rank <= max_rank <= minimum(N, M).

  • rtol (float) – Relative tolerance for truncation. Remove singular values satisfying sigma < maximum(atol, rtol*sigma1). Cannot be used for stacked A (len(A.shape) > 2).

  • atol (float) – Absolute tolerance for truncation. Remove singular values satisfying sigma < maximum(atol, rtol*sigma1). Cannot be used for stacked A (len(A.shape) > 2).

  • max_rank (int)

Returns:

  • U (NDArray) – Left singular vectors. shape=(…, N, r). U.T @ U = identity matrix

  • ss (NDArray) – Singular values. Non-negative. shape=(…, r).

  • Vt (NDArray) – Right singular vectors. shape=(…, r, M) Vt @ Vt.T = identity matrix

Return type:

t3toolbox.backend.common.typ.Tuple[NDArray, NDArray, NDArray]

Examples

Default (no truncation): a full thin SVD, applied over the leading stack axes. The factors reconstruct A and are orthonormal:

>>> import numpy as np
>>> import t3toolbox.backend.linalg as linalg
>>> np.random.seed(0)
>>> A = np.random.randn(2, 3, 4, 55, 70)              # stack_shape=(2,3,4), matrices 55x70
>>> U, ss, Vt = linalg.truncated_svd(A)
>>> print(U.shape, ss.shape, Vt.shape)               # r = min(55, 70) = 55
(2, 3, 4, 55, 55) (2, 3, 4, 55) (2, 3, 4, 55, 70)
>>> A2 = np.einsum('...ix,...x,...xj->...ij', U, ss, Vt)
>>> print(np.allclose(A, A2))                         # U @ diag(ss) @ Vt == A
True
>>> print(np.allclose(np.einsum('...ix,...iy->...xy', U, U), np.eye(U.shape[-1])))   # U^T U = I
True
>>> print(np.allclose(np.einsum('...xj,...yj->...xy', Vt, Vt), np.eye(Vt.shape[-2])))  # Vt Vt^T = I
True
>>> print(bool(np.all(ss >= 0.0)))                    # singular values non-negative
True

max_rank caps the kept rank r (changes the output shapes), and works on a stack:

>>> np.random.seed(0)
>>> A = np.random.randn(2, 3, 4, 55, 70)
>>> U, ss, Vt = linalg.truncated_svd(A, max_rank=5)
>>> print(U.shape, ss.shape, Vt.shape)               # r capped at 5
(2, 3, 4, 55, 5) (2, 3, 4, 5) (2, 3, 4, 5, 70)

rtol drops singular values below rtol * sigma1 – a deliberate approximation, so we show the kept rank and assert the accuracy bound rather than equality. Feed a graded spectrum (a Hilbert-like matrix) so the tolerance actually truncates; rtol/atol require an unstacked A:

>>> A = np.array([[1.0 / (ii + jj) for jj in range(1, 70)] for ii in range(1, 55)])  # graded spectrum
>>> U, ss, Vt = linalg.truncated_svd(A, rtol=1e-2)
>>> print(ss.shape[-1])                               # rtol=1e-2 keeps 3 singular values
3
>>> A2 = np.einsum('ix,x,xj->ij', U, ss, Vt)
>>> rel_err = np.linalg.norm(A - A2, 2) / np.linalg.norm(A, 2)
>>> print(bool(rel_err < 1e-2))                       # relative 2-norm error below rtol
True
>>> ss_full = np.linalg.svd(A, compute_uv=False)      # accuracy bound (generalized Oseledets):
>>> dropped = np.linalg.norm(ss_full[ss.shape[-1]:])  #   ||A - A2||_F <= sqrt(dropped energy)
>>> print(bool(np.linalg.norm(A - A2) <= dropped + 1e-12))
True

min_rank is a floor that overrides the tolerance – here it forces rank 10 even though rtol alone would keep only 3, driving the error far below the tolerance:

>>> U, ss, Vt = linalg.truncated_svd(A, rtol=1e-2, min_rank=10)
>>> print(ss.shape[-1])                               # floored at min_rank=10
10

Gotcha: rtol/atol on a stacked A raises – the kept rank could differ per stack element, giving ragged shapes. Unstack first, then truncate each matrix:

>>> np.random.seed(0)
>>> A = np.random.randn(2, 3, 4, 55, 70)
>>> linalg.truncated_svd(A, rtol=1e-2)
Traceback (most recent call last):
    ...
ValueError