Skip to content

linalg

Linear algebra utilities (numpy / cupy compatible).

Functions#

truncated_qr #

truncated_qr(
    matrix: NDArray, cutoff: float, xp: ModuleType = np
) -> NDArray

QR with SVD-based rank truncation, returning only the isometry.

Decomposes matrix = Q @ R, then truncates via SVD on R, discarding singular values below cutoff * sigma_max.

When cutoff <= 0 this falls back to a plain QR (no truncation).

PARAMETER DESCRIPTION
matrix

Input matrix of shape (m, n).

TYPE: NDArray

cutoff

Relative singular-value threshold. Singular values satisfying s < cutoff * s_max are discarded.

TYPE: float

xp

Array module (numpy or cupy); defaults to numpy.

TYPE: ModuleType DEFAULT: numpy

RETURNS DESCRIPTION
NDArray

The truncated isometry Q of shape (m, rank).
Source code in src/src_method/utils/linalg.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def truncated_qr(matrix: NDArray, cutoff: float, xp: ModuleType = np) -> NDArray:
    """QR with SVD-based rank truncation, returning only the isometry.

    Decomposes ``matrix = Q @ R``, then truncates via SVD on R,
    discarding singular values below ``cutoff * sigma_max``.

    When ``cutoff <= 0`` this falls back to a plain QR (no truncation).

    Args:
        matrix: Input matrix of shape (m, n).
        cutoff: Relative singular-value threshold.  Singular values
            satisfying ``s < cutoff * s_max`` are discarded.
        xp: Array module (``numpy`` or ``cupy``); defaults to numpy.

    Returns:
        The truncated isometry Q of shape ``(m, rank)``.
    """
    m, n = matrix.shape
    transpose = m < n

    if cutoff <= 0:
        return xp.linalg.qr(matrix)[0]

    Q, R = xp.linalg.qr(matrix.T if transpose else matrix)
    R_np = R.get() if hasattr(R, "get") else np.asarray(R)
    U, S, _ = np.linalg.svd(R_np.T if transpose else R_np, full_matrices=False)
    rank = max(1, int((cutoff * S[0] <= S).sum()))

    if transpose:
        Q_trunc = U[:, :rank]
    else:
        U_trunc = xp.asarray(U[:, :rank]) if xp is not np else U[:, :rank]
        Q_trunc = Q @ U_trunc

    if xp is not np:
        return xp.asarray(Q_trunc)
    return Q_trunc