deepbox/linalg

Decompositions

LU, Cholesky, QR, SVD, eigendecompositions, Schur, Hessenberg and polar factorizations of 2-D matrices.

Factorizations

Overview

Every function takes 2-D tensors, reads them as float64 and returns float64 tensors on the default device, whatever the input dtype. NaN or infinite input throws DataValidationError, and string or complex input throws DTypeError.
svd(a, fullMatrices = true) returns [U, s, Vt] with the singular values as a vector. qr returns [Q, R] in reduced mode by default. lu returns [P, L, U] with P A = L U. cholesky returns the lower factor, or the upper one with { upper: true }, and throws for matrices that are not symmetric positive definite. hessenberg returns [H, Q] and polar returns [U, P].
eig supports real eigenvalues only. A non-symmetric matrix with complex eigenvalues throws InvalidParameterError. It returns [values, vectors] with unit-norm eigenvectors as columns. Symmetric input takes a faster path and gives ascending eigenvalues, and eigh and eigvalsh are the dedicated functions when you know the matrix is symmetric.

Types

CholeskyOptions

type

Options for cholesky.

- upper: if true, return the upper triangular factor U with A = Uᵀ * U instead of the lower triangular factor L with A = L * Lᵀ (default: false).

readonly upper?: boolean

EigOptions

type

Options for eig and eigvals.

- maxIter: maximum number of QR sweeps spent on a single eigenvalue before a ConvergenceError is thrown (default: 300). A trailing 2×2 block is solved in closed form and needs no sweeps. - tol: relative threshold below which a sub-diagonal entry is treated as zero (default: machine epsilon, about 2.2e-16). Larger values deflate earlier at the cost of accuracy.

readonly maxIter?: number
readonly tol?: number

Functions

cholesky

function cholesky(a: Tensor, options?: CholeskyOptions): Tensor;

Cholesky decomposition.

Details

Factorizes symmetric positive definite matrix A into L * L^T where L is lower triangular.

Algorithm: Cholesky-Banachiewicz algorithm Time Complexity: O(N³/3) - approximately half the cost of LU decomposition Space Complexity: O(N²) for output matrix

Mathematical Background: For a symmetric positive definite matrix A, the Cholesky decomposition is: A = L * L^T where L is lower triangular with positive diagonal elements.

Numerical Stability: - Only works for positive definite matrices (all eigenvalues > 0) - More stable than LU for this class of matrices - No pivoting required due to positive definiteness

Parameters:

Parameters

a Tensor: Symmetric positive definite matrix of shape (N, N)
optionsoptional CholeskyOptions: Optional settings (see CholeskyOptions)

Returns

Tensor

Throws

ShapeErrorIf input is not 2D square matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf matrix is not symmetric
DataValidationErrorIf matrix is not positive definite
DataValidationErrorIf input contains non-finite values (NaN, Infinity)

ExampleTypeScript

import { cholesky } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; // Positive definite matrixconst A = tensor([[4, 12, -16], [12, 37, -43], [-16, -43, 98]]);const L = cholesky(A); // Verify: A ≈ L @ L^T

eig

function eig(a: Tensor, options?: EigOptions): [ Tensor, Tensor ];

Compute eigenvalues and eigenvectors of a square matrix.

Details

Solves A * v = λ * v where λ are eigenvalues and v are eigenvectors.

Algorithm: - Symmetric matrices: Householder tridiagonalization + implicit QL (same as eigh). Eigenvalues are returned in ascending order. - General matrices: balancing, Hessenberg reduction, Francis double-shift QR to the real Schur form, then back-substitution for the eigenvectors. Eigenvalues are returned in the order they appear on the diagonal of the Schur form.

Limitations: - Only real eigenvalues are supported. Non-symmetric matrices whose spectrum includes complex eigenvalues cause an InvalidParameterError to be thrown. - A conjugate pair whose imaginary part is at rounding level (at most 10 * N * eps * max|A| after balancing, as for the repeated zero eigenvalue of a low-rank matrix) is reported as a repeated real eigenvalue. - For symmetric/Hermitian matrices, use eigh() for better performance - Eigenvectors of a defective matrix (fewer independent eigenvectors than eigenvalues) are nearly parallel, as with numpy.linalg.eig.

Parameters:

Parameters

a Tensor: Square matrix of shape (N, N)
optionsoptional EigOptions: Optional configuration overrides (see EigOptions)

Returns

[ Tensor, Tensor ]

Throws

ShapeErrorIf input is not square matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
InvalidParameterErrorIf matrix has complex eigenvalues, or an option is invalid
ConvergenceErrorIf the QR iteration does not converge within maxIter sweeps

ExampleTypeScript

import { eig } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [2, 1]]);const [eigenvalues, eigenvectors] = eig(A); // Verify: A @ eigenvectors[:,i] ≈ eigenvalues[i] * eigenvectors[:,i]

eigh

function eigh(a: Tensor): [ Tensor, Tensor ];

Compute eigenvalues and eigenvectors of a symmetric/Hermitian matrix.

More efficient than eig() for symmetric matrices.

Parameters:

Parameters

a Tensor: Symmetric matrix of shape (N, N) Returns: [eigenvalues, eigenvectors] - Eigenvalues are real and in ascending order - Eigenvectors are orthonormal columns

Returns

[ Tensor, Tensor ]

Throws

ShapeErrorIf input is not square matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input is not symmetric
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
ConvergenceErrorIf the QL iteration does not converge

ExampleTypeScript

import { eigh } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [2, 1]]);  // Symmetricconst [eigenvalues, eigenvectors] = eigh(A);

eigvals

function eigvals(a: Tensor, options?: EigOptions): Tensor;

Compute eigenvalues only (faster than eig).

Skips the eigenvector computation. Eigenvalues of symmetric input are returned in ascending order; for general input they follow the order of the real Schur form.

Parameters:

Parameters

a Tensor: Square matrix of shape (N, N)
optionsoptional EigOptions: Optional configuration overrides (see EigOptions) Returns: eigenvalues - Array of real eigenvalues Limitations: - Matrices with complex eigenvalues are not supported and will throw InvalidParameterError

Returns

Tensor

Throws

ShapeErrorIf input is not square matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
InvalidParameterErrorIf matrix has complex eigenvalues, or an option is invalid
ConvergenceErrorIf the QR iteration does not converge within maxIter sweeps

ExampleTypeScript

import { eigvals } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [2, 1]]);const eigenvalues = eigvals(A);console.log(eigenvalues);  // [-1, 3]

eigvalsh

function eigvalsh(a: Tensor): Tensor;

Compute eigenvalues only of symmetric matrix (faster than eigh).

Parameters:

Parameters

a Tensor: Symmetric square matrix of shape (N, N) Returns: eigenvalues - Array of real eigenvalues in ascending order

Returns

Tensor

Throws

ShapeErrorIf input is not square matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input is not symmetric
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
ConvergenceErrorIf the QL iteration does not converge

ExampleTypeScript

import { eigvalsh } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [2, 1]]);const eigenvalues = eigvalsh(A);console.log(eigenvalues);  // [-1, 3]

hessenberg

function hessenberg(a: Tensor): [ Tensor, Tensor ];

Upper Hessenberg decomposition of a square matrix.

Decomposes A into A = Q * H * Qᵀ where: - Q is an orthogonal matrix - H is an upper Hessenberg matrix (zero below the first sub-diagonal)

Uses Householder reflections. Unlike scipy.linalg.hessenberg, both factors are always returned.

Time Complexity: O(N³)

Parameters

a Tensor: Input square matrix of shape (N, N)

Returns

[ Tensor, Tensor ][H, Q] where H is upper Hessenberg and Q is orthogonal

Throws

ShapeErrorIf input is not a square 2D matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)

ExampleTypeScript

import { hessenberg } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2, 3], [4, 5, 6], [7, 8, 9]]);const [H, Q] = hessenberg(A);// Q * H * Q^T ≈ A// H is upper Hessenberg (zero below first sub-diagonal)

lu

function lu(a: Tensor): [ Tensor, Tensor, Tensor ];

LU decomposition with partial pivoting.

Details

Factorizes matrix A into P * A = L * U where: - P is a permutation matrix - L is lower triangular with unit diagonal - U is upper triangular

Algorithm: Gaussian elimination with partial pivoting

Convention: P is the row permutation applied to A (P @ A = L @ U). This is the transpose of the p returned by scipy.linalg.lu, which satisfies A = p @ l @ u.

Parameters:

Parameters

a Tensor: Input matrix of shape (M, N) Returns: [P, L, U] - P: Permutation matrix of shape (M, M) - L: Lower triangular matrix of shape (M, K) where K = min(M, N) - U: Upper triangular matrix of shape (K, N) Requirements: - Input must be 2D matrix Properties: - P @ A = L @ U - L has unit diagonal (L[i,i] = 1) - Partial pivoting ensures numerical stability

Returns

[ Tensor, Tensor, Tensor ]

Throws

ShapeErrorIf input is not 2D
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity) Rank-deficient matrices (with zero pivots) are accepted. The factorization still produces valid P, L, U factors.

ExampleTypeScript

import { lu } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[2, 1, 1], [4, 3, 3], [8, 7, 9]]);const [P, L, U] = lu(A); // Verify: P @ A ≈ L @ U

polar

function polar(a: Tensor, side?: "right" | "left"): [ Tensor, Tensor ];

Polar decomposition of a matrix.

Details

Decomposes A into A = U * P where: - U has orthonormal columns (a unitary matrix when A is square, M >= N) or orthonormal rows (when M < N) - P is a symmetric positive semi-definite matrix (right polar factor)

Or A = S * U where: - S is a symmetric positive semi-definite matrix (left polar factor) - U is as above

Uses the SVD-based algorithm: given the reduced SVD A = W Σ Vᵀ, U = W Vᵀ (orthogonal factor) P = V Σ Vᵀ (right polar factor) S = W Σ Wᵀ (left polar factor)

For a rank-deficient A the orthogonal factor U is not unique; the one returned follows from the singular vectors that belong to the zero singular values.

Time Complexity: O(M * N * min(M,N)) dominated by SVD

Parameters

a Tensor: Input matrix of shape (M, N)
sideoptional "right" | "left": 'right' (default) returns [U, P] where A = U*P, 'left' returns [S, U] where A = S*U

Returns

[ Tensor, Tensor ][U, P] of shapes (M, N) and (N, N), or [S, U] of shapes (M, M) and (M, N)

Throws

ShapeErrorIf input is not 2D
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
InvalidParameterErrorIf side is not 'right' or 'left'
ConvergenceErrorIf the underlying SVD does not converge

ExampleTypeScript

import { polar } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [3, 4]]);const [U, P] = polar(A);// U is unitary: U * U^T ≈ I// P is positive semi-definite// U * P ≈ A

qr

function qr(a: Tensor, mode?: "reduced" | "complete"): [ Tensor, Tensor ];

QR decomposition.

Details

Factorizes matrix A into Q * R where Q is orthogonal and R is upper triangular.

Algorithm: Householder reflections, with the sign convention of LAPACK (and therefore NumPy): the reflection for column k maps it to -sign(x₀) * ‖x‖ * e₁, and a column that is already zero below the diagonal is left alone. As a result the diagonal of R is not necessarily positive.

Parameters:

Parameters

a Tensor: Input matrix of shape (M, N)
modeoptional "reduced" | "complete": Decomposition mode: - 'reduced' (default): Q has shape (M, K), R has shape (K, N) where K = min(M, N) - 'complete': Q has shape (M, M), R has shape (M, N) Returns: [Q, R] - Q: Orthogonal matrix (Q^T * Q = I) - R: Upper triangular matrix Requirements: - Input must be 2D matrix Properties: - A = Q @ R - Q is orthogonal: Q^T @ Q = I - R is upper triangular - In 'reduced' mode only the first K columns of Q are formed, so the cost is O(M * N * K) rather than O(M²) memory for tall inputs.

Returns

[ Tensor, Tensor ]

Throws

ShapeErrorIf input is not 2D
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
InvalidParameterErrorIf mode is not 'reduced' or 'complete'

ExampleTypeScript

import { qr } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]);const [Q, R] = qr(A); // Verify: A ≈ Q @ R// Verify: Q is orthogonal (Q^T @ Q ≈ I)

schur

function schur(a: Tensor): [ Tensor, Tensor ];

Real Schur decomposition of a square matrix.

Details

Decomposes A into A = Q * T * Qᵀ where: - Q is an orthogonal matrix (Schur vectors) - T is a quasi-upper triangular matrix (Schur form) - Real eigenvalues appear on the diagonal - Complex conjugate pairs appear in 2×2 blocks on the diagonal; such a block has equal diagonal entries and off-diagonal entries of opposite sign

Reduces A to Hessenberg form with Householder reflections, then runs the QR algorithm with implicit double shifts (Francis iteration) and deflation. The result has the same structure as scipy.linalg.schur(a, output="real"); the Schur vectors are unique only up to sign.

Time Complexity: O(N³)

Parameters

a Tensor: Input square matrix of shape (N, N)

Returns

[ Tensor, Tensor ][T, Q] where T is the Schur form and Q contains the Schur vectors

Throws

ShapeErrorIf input is not a square 2D matrix
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
ConvergenceErrorIf the QR iteration does not converge

ExampleTypeScript

import { schur } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [3, 4]]);const [T, Q] = schur(A);// Q * T * Q^T ≈ A// T is quasi-upper triangular

svd

function svd(a: Tensor, fullMatrices?: boolean): [ Tensor, Tensor, Tensor ];

Singular Value Decomposition.

Details

Factorizes matrix A into three matrices: A = U * Σ * V^T

Algorithm: Golub–Reinsch. Householder reflections reduce A to bidiagonal form, then implicit-shift QR iteration diagonalizes the bidiagonal matrix. Wide matrices (M < N) are handled by decomposing Aᵀ and swapping the factors.

Numerical Stability: Works on A directly and never forms AᵀA, so small singular values keep their relative accuracy far better than a normal-equations approach. Matrices with extreme entry magnitudes are prescaled by a power of two (an exact operation). When fullMatrices is true, the extra columns of U (or rows of Vᵀ) are an orthonormal completion and are not unique.

Parameters:

Parameters

a Tensor: Input matrix of shape (M, N)
fullMatricesoptional boolean: If true (default), U has shape (M, M) and Vt has shape (N, N). If false, U has shape (M, K) and Vt has shape (K, N) where K = min(M, N) Returns: [U, s, Vt] - U: Left singular vectors of shape (M, M) or (M, K) - s: Singular values of shape (K,) in descending order - Vt: Right singular vectors transposed, of shape (N, N) or (K, N) Requirements: - Input must be 2D matrix - Works with any M x N matrix, including rank-deficient and empty ones Properties: - A = U @ diag(s) @ Vt (with the first K columns of U and rows of Vt) - U and Vt have orthonormal columns and rows respectively - Singular values are non-negative and sorted in descending order - Singular vectors are unique only up to a common sign flip per pair

Returns

[ Tensor, Tensor, Tensor ]

Throws

ShapeErrorIf input is not 2D
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
ConvergenceErrorIf the QR iteration does not converge

ExampleTypeScript

import { svd } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const A = tensor([[1, 2], [3, 4], [5, 6]]);const [U, s, Vt] = svd(A); console.log(U.shape);   // [3, 3] if fullMatrices=true, [3, 2] if falseconsole.log(s.shape);   // [2]console.log(Vt.shape);  // [2, 2] // Reconstruction: A ≈ U[:, :2] @ diag(s) @ Vt

svdvals

function svdvals(a: Tensor): Tensor;

Singular values of a only (no U/V).

Skips the accumulation of the singular vectors, so it is considerably cheaper than svd. Entries of extreme magnitude are prescaled by a power of two, as in svd.

Parameters

a Tensor: Input matrix of shape (M, N)

Returns

TensorSingular values of shape (min(M, N),) in descending order

Throws

ShapeErrorIf input is not 2D
DTypeErrorIf input has string or complex dtype
DataValidationErrorIf input contains non-finite values (NaN, Infinity)
ConvergenceErrorIf the QR iteration does not converge

ExampleTypeScript

import { svdvals } from 'deepbox/linalg';import { tensor } from 'deepbox/ndarray'; const s = svdvals(tensor([[3, 0], [0, -2]]));// [3, 2]

linalg-decompositions.tsTypeScript

import { cholesky, eig, hessenberg, lu, polar, qr, schur, svd } from "deepbox/linalg";import { tensor } from "deepbox/ndarray"; const a = tensor([  [4, 1],  [1, 3],]); const [q, r] = qr(a);const [u, s, vt] = svd(a);const [values, vectors] = eig(a);const [p, l, upper] = lu(a); console.log(q.shape, r.shape, u.shape, vt.shape);console.log(s.toString(), values.toString(), vectors.shape);console.log(p.shape, l.shape, upper.shape);console.log(cholesky(a).toString());console.log(schur(a)[0].shape, hessenberg(a)[0].shape, polar(a)[0].shape);

Found a mistake? Open an issue.