Linear Algebra¶

Linear algebra primitives with multi-backend support.

This module provides backend-agnostic linear algebra operations for NumPy, PyTorch, and JAX, with robust fallbacks.

Decompositions¶

QR decomposition (with column pivoting option)
Cholesky decomposition (with safe fallback)
SVD (singular value decomposition)
Eigendecomposition for symmetric matrices

Solvers¶

solve_triangular: Efficient triangular system solver
solve_cholesky: Solve via Cholesky factorization
solve_qr: Solve via QR decomposition
lstsq: Least squares solution

Matrix Operations¶

safe_inverse: Matrix inverse with regularization
woodbury_inverse: Efficient update of inverse
quadratic_form: x’Ax computation
log_determinant: Log of matrix determinant

Examples

>>> from aurora.core.linalg import safe_cholesky, solve_cholesky
>>>
>>> # Robust Cholesky with fallback
>>> L = safe_cholesky(A)
>>> x = solve_cholesky(L, b)
>>>
>>> # QR for least squares
>>> Q, R = qr_decomposition(X)
>>> beta = solve_qr(Q, R, y)

aurora.core.linalg.qr_decomposition(A, *, mode='reduced', pivoting=False)[source]¶

QR decomposition with optional column pivoting.

Decomposes A = QR where Q is orthogonal and R is upper triangular. With pivoting, AP = QR where P is a permutation matrix.

Parameters:

A (ndarray of shape (m, n)) – Matrix to decompose.
mode ({'reduced', 'complete', 'r', 'raw'}, default='reduced') –
- ‘reduced’: Return Q (m, k) and R (k, n), where k = min(m, n)
- ’complete’: Return Q (m, m) and R (m, n)
- ’r’: Return only R
- ’raw’: Return (H, tau) for advanced use
pivoting (bool, default=False) – If True, use column pivoting and return permutation indices.

Returns:

Q (ndarray) – Orthogonal matrix (unless mode=’r’ or ‘raw’).
R (ndarray) – Upper triangular matrix.
P (ndarray (only if pivoting=True)) – Permutation indices such that A[:, P] = Q @ R.

Return type:

tuple[ndarray, ndarray] | tuple[ndarray, ndarray, ndarray]

Examples

>>> A = np.array([[1, 2], [3, 4], [5, 6]])
>>> Q, R = qr_decomposition(A)
>>> np.allclose(A, Q @ R)
True

>>> # With pivoting for rank-deficient matrices
>>> Q, R, P = qr_decomposition(A, pivoting=True)
>>> np.allclose(A[:, P], Q @ R)
True

aurora.core.linalg.safe_cholesky(A, *, lower=True, max_tries=5, jitter=1e-6)[source]¶

Cholesky decomposition with robustness to near-singular matrices.

Attempts standard Cholesky, and if it fails, adds progressively larger diagonal jitter until successful.

Parameters:

A (ndarray of shape (n, n)) – Symmetric positive semi-definite matrix.
lower (bool, default=True) – If True, return lower triangular L. If False, return upper.
max_tries (int, default=5) – Maximum number of attempts with increasing jitter.
jitter (float, default=1e-6) – Initial diagonal jitter to add on failure.

Returns:

L – Cholesky factor such that A ≈ L @ L.T.

Return type:

ndarray of shape (n, n)

Raises:

np.linalg.LinAlgError – If decomposition fails after max_tries attempts.

Examples

>>> A = np.array([[4, 2], [2, 1.001]])  # Nearly singular
>>> L = safe_cholesky(A)
>>> np.allclose(A, L @ L.T, atol=1e-5)
True

aurora.core.linalg.svd(A, *, full_matrices=False, compute_uv=True)[source]¶

Singular Value Decomposition.

Decomposes A = U @ S @ V.T where U and V are orthogonal and S is diagonal.

Parameters:

A (ndarray of shape (m, n)) – Matrix to decompose.
full_matrices (bool, default=False) – If True, return full U (m, m) and Vh (n, n). If False, return reduced U (m, k) and Vh (k, n).
compute_uv (bool, default=True) – If True, return U, S, Vh. If False, return only S.

Returns:

U (ndarray) – Left singular vectors (if compute_uv=True).
S (ndarray) – Singular values (sorted in descending order).
Vh (ndarray) – Right singular vectors transposed (if compute_uv=True).

Return type:

tuple[ndarray, ndarray, ndarray] | ndarray

Examples

>>> A = np.random.randn(5, 3)
>>> U, S, Vh = svd(A)
>>> np.allclose(A, U @ np.diag(S) @ Vh)
True

aurora.core.linalg.eigh(A, *, subset_by_index=None)[source]¶

Eigendecomposition for symmetric/Hermitian matrices.

Parameters:

A (ndarray of shape (n, n)) – Symmetric matrix.
subset_by_index (tuple of (lo, hi), optional) – Only compute eigenvalues/vectors in index range [lo, hi].

Returns:

eigenvalues (ndarray of shape (n,) or (hi-lo+1,)) – Eigenvalues in ascending order.
eigenvectors (ndarray of shape (n, n) or (n, hi-lo+1)) – Eigenvectors as columns.

Return type:

tuple[ndarray, ndarray]

Examples

>>> A = np.array([[2, 1], [1, 2]])
>>> vals, vecs = eigh(A)
>>> np.allclose(A @ vecs, vecs @ np.diag(vals))
True

aurora.core.linalg.solve_triangular(A, b, *, lower=True, trans=False)[source]¶

Solve a triangular system Ax = b.

Parameters:

A (ndarray of shape (n, n)) – Triangular matrix.
b (ndarray of shape (n,) or (n, k)) – Right-hand side.
lower (bool, default=True) – True if A is lower triangular, False if upper.
trans (bool, default=False) – If True, solve A.T @ x = b instead.

Returns:

x – Solution to the system.

Return type:

ndarray

aurora.core.linalg.solve_cholesky(L, b, *, lower=True)[source]¶

Solve a system using pre-computed Cholesky factor.

Solves A x = b where A = L @ L.T (lower=True) or A = U.T @ U (lower=False).

Parameters:

L (ndarray of shape (n, n)) – Cholesky factor.
b (ndarray of shape (n,) or (n, k)) – Right-hand side.
lower (bool, default=True) – True if L is lower triangular.

Returns:

x – Solution to the system.

Return type:

ndarray

Examples

>>> A = np.array([[4, 2], [2, 5]])
>>> L = safe_cholesky(A)
>>> b = np.array([1, 2])
>>> x = solve_cholesky(L, b)
>>> np.allclose(A @ x, b)
True

aurora.core.linalg.solve_qr(Q, R, b)[source]¶

Solve a least squares problem using QR factors.

Solves min ||Ax - b||² where A = QR.

Parameters:

Q (ndarray of shape (m, k)) – Orthogonal factor.
R (ndarray of shape (k, n)) – Upper triangular factor.
b (ndarray of shape (m,)) – Right-hand side.

Returns:

x – Least squares solution.

Return type:

ndarray of shape (n,)

aurora.core.linalg.lstsq(A, b, *, rcond=None)[source]¶

Least squares solution to Ax = b.

Parameters:

A (ndarray of shape (m, n)) – Design matrix.
b (ndarray of shape (m,) or (m, k)) – Right-hand side.
rcond (float, optional) – Cutoff ratio for small singular values.

Returns:

x (ndarray) – Least squares solution.
residuals (ndarray) – Sum of squared residuals (if m > n).
rank (int) – Effective rank of A.
s (ndarray) – Singular values of A.

Return type:

tuple[ndarray, ndarray, int, ndarray]

aurora.core.linalg.safe_inverse(A, *, rcond=1e-15)[source]¶

Compute matrix inverse with regularization for near-singular matrices.

Uses SVD-based pseudoinverse with condition number thresholding.

Parameters:

A (ndarray of shape (n, n)) – Matrix to invert.
rcond (float, default=1e-15) – Cutoff for small singular values. Singular values smaller than rcond * largest_singular_value are set to zero.

Returns:

A_inv – (Pseudo)inverse of A.

Return type:

ndarray of shape (n, n)

Examples

>>> A = np.array([[1, 2], [2, 4.001]])  # Nearly singular
>>> A_inv = safe_inverse(A)

aurora.core.linalg.woodbury_inverse(A_inv, U, C, V)[source]¶

Woodbury matrix identity for efficient inverse updates.

Computes (A + UCV)^{-1} given A^{-1} efficiently.

The Woodbury identity states: (A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1} + VA^{-1}U)^{-1}VA^{-1}

Parameters:

A_inv (ndarray of shape (n, n)) – Inverse of A.
U (ndarray of shape (n, k)) – Update matrix.
C (ndarray of shape (k, k)) – Core matrix.
V (ndarray of shape (k, n)) – Update matrix.

Returns:

result – (A + UCV)^{-1}

Return type:

ndarray of shape (n, n)

Examples

>>> A = np.eye(3)
>>> U = np.array([[1], [0], [0]])
>>> C = np.array([[1]])
>>> V = np.array([[1, 0, 0]])
>>> new_inv = woodbury_inverse(np.eye(3), U, C, V)

aurora.core.linalg.quadratic_form(x, A, y=None)[source]¶

Compute quadratic form x’Ay or x’Ax.

Parameters:

x (ndarray of shape (n,)) – Left vector.
A (ndarray of shape (n, n)) – Matrix.
y (ndarray of shape (n,), optional) – Right vector. If None, computes x’Ax.

Returns:

result – The quadratic form value.

Return type:

float

aurora.core.linalg.log_determinant(A, *, method='auto')[source]¶

Compute log of matrix determinant in a numerically stable way.

Parameters:

A (ndarray of shape (n, n)) – Positive definite matrix.
method ({'cholesky', 'svd', 'auto'}, default='auto') – Method to use. ‘auto’ tries Cholesky first, falls back to SVD.

Returns:

logdet – Log of determinant of A.

Return type:

float

Examples

>>> A = np.array([[2, 1], [1, 2]])
>>> logdet = log_determinant(A)
>>> np.isclose(logdet, np.log(np.linalg.det(A)))
True

aurora.core.linalg.matrix_rank(A, *, tol=None)[source]¶

Compute matrix rank.

Parameters:

A (ndarray) – Input matrix.
tol (float, optional) – Threshold below which singular values are considered zero.

Returns:

rank – Numerical rank of the matrix.

Return type:

int

aurora.core.linalg.condition_number(A, *, p=2)[source]¶

Compute matrix condition number.

Parameters:

A (ndarray) – Input matrix.
p ({1, 2, inf, 'fro'}, default=2) – Order of the norm.

Returns:

cond – Condition number of A.

Return type:

float

aurora.core.linalg.weighted_condition_number(X, weights=None)[source]¶

Compute condition number of weighted design matrix X^T W X.

Parameters:

X (ndarray, shape (n, p)) – Design matrix.
weights (ndarray, shape (n,), optional) – Observation weights. If None, computes condition of X^T X.

Returns:

kappa – Condition number κ = σ_max / σ_min of X^T W X.

Return type:

float