API Reference¶

Data Structures¶

class linearsystems.LinearSystemResult(x, residual, success=True, message='', n_iter=0)¶

Bases: object

Result of a linear system solve.

Parameters:

x (ndarray)
residual (float)
success (bool)
message (str)
n_iter (int)

x¶

Solution vector of shape (n,).

Type:: np.ndarray

residual¶

Euclidean norm of the residual ‖Ax − b‖.

Type:: float

success¶

Whether the solve succeeded.

Type:: bool

message¶

Descriptive status message.

Type:: str

n_iter¶

Number of iterations (iterative solvers only).

Type:: int

Dense Direct Solvers¶

linearsystems.lslrg(a, b)¶

Solve a general real linear system Ax = b.

Uses numpy.linalg.solve() (LU decomposition with partial pivoting).

Parameters:

a (np.ndarray) – Real coefficient matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,) or matrix of shape (n, k).

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.lscrg(a, b)¶

Solve a general complex linear system Ax = b.

Uses numpy.linalg.solve() on complex arrays.

Parameters:

a (np.ndarray) – Complex coefficient matrix of shape (n, n).
b (np.ndarray) – Complex right-hand side vector of shape (n,) or matrix of shape (n, k).

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.lslpb(a, b)¶

Solve a positive definite banded linear system Ax = b.

Accepts a as a full square matrix or as a banded storage array (lower + upper + 1, n) compatible with scipy.linalg.solve_banded(). When a full matrix is provided it is converted to banded form with bandwidth inferred from the non-zero pattern.

Parameters:

a (np.ndarray) – Coefficient matrix of shape (n, n) (full) or banded storage of shape (lower + upper + 1, n).
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If the system dimensions are inconsistent.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.lslsf(a, b)¶

Solve a symmetric positive definite linear system Ax = b.

Uses scipy.linalg.solve() with assume_a='pos' to exploit the Cholesky factorization.

Parameters:

a (np.ndarray) – Symmetric positive definite matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,) or matrix of shape (n, k).

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular or not positive definite.

linearsystems.lsltr(a, b, lower=True)¶

Solve a triangular linear system Ax = b.

Uses scipy.linalg.solve_triangular().

Parameters:

a (np.ndarray) – Triangular matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,) or matrix of shape (n, k).
lower (bool) – If True (default), treat a as lower triangular; otherwise upper triangular.

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.lsltc(a, b, lower=True)¶

Solve a complex triangular linear system Ax = b.

Uses scipy.linalg.solve_triangular() on complex arrays.

Parameters:

a (np.ndarray) – Complex triangular matrix of shape (n, n).
b (np.ndarray) – Complex right-hand side vector of shape (n,) or matrix of shape (n, k).
lower (bool) – If True (default), treat a as lower triangular; otherwise upper triangular.

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

Factorization-Based Solvers¶

linearsystems.lufact(a)¶

Compute the LU factorization of a square matrix.

Wraps scipy.linalg.lu_factor().

Parameters:

a (np.ndarray) – Real matrix of shape (n, n) to factor.

Returns:

(lu, piv) compatible with lusolve() /: scipy.linalg.lu_solve().

Return type:

tuple

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.lusolve(lu_piv, b)¶

Solve Ax = b using a pre-computed LU factorization.

Wraps scipy.linalg.lu_solve().

Parameters:

lu_piv (tuple) – (lu, piv) tuple returned by lufact().
b (np.ndarray) – Right-hand side vector of shape (n,) or matrix of shape (n, k).

Returns:

Solution vector x of shape (n,).

Return type:

np.ndarray

linearsystems.cholfact(a)¶

Compute the Cholesky factorization of a symmetric positive definite matrix.

Wraps scipy.linalg.cho_factor().

Parameters:

a (np.ndarray) – Symmetric positive definite matrix of shape (n, n).

Returns:

(c, lower) compatible with cholsolve() /: scipy.linalg.cho_solve().

Return type:

tuple

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is not positive definite.

linearsystems.cholsolve(c_low, b)¶

Solve Ax = b using a pre-computed Cholesky factorization.

Wraps scipy.linalg.cho_solve().

Parameters:

c_low (tuple) – (c, lower) tuple returned by cholfact().
b (np.ndarray) – Right-hand side vector of shape (n,) or matrix of shape (n, k).

Returns:

Solution vector x of shape (n,).

Return type:

np.ndarray

linearsystems.qrfact(a)¶

Compute the QR factorization of a matrix.

Wraps scipy.linalg.qr() with mode='economic'.

Parameters:

a (np.ndarray) – Matrix of shape (m, n) to factor.

Returns:

(Q, R) where Q is (m, k) orthogonal and: R is (k, n) upper triangular, k = min(m, n).

Return type:

tuple

linearsystems.qrsolve(qr_result, b)¶

Solve (or least-squares fit) Ax = b using a pre-computed QR factorization.

Parameters:

qr_result (tuple) – (Q, R) tuple returned by qrfact().
b (np.ndarray) – Right-hand side vector of shape (m,).

Returns:

Solution (or least-squares solution) vector of shape (n,).

Return type:

np.ndarray

Iterative Solvers¶

linearsystems.lsitcg(a, b, tol=1e-10, maxiter=None)¶

Solve a symmetric positive definite system using Conjugate Gradient.

Wraps scipy.sparse.linalg.cg().

Parameters:

a (np.ndarray) – Symmetric positive definite matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.
maxiter (int | None) – Maximum number of iterations. Defaults to None (scipy default).

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.

linearsystems.lsibcg(a, b, tol=1e-10, maxiter=None)¶

Solve a linear system using BiConjugate Gradient.

Wraps scipy.sparse.linalg.bicg().

Parameters:

a (np.ndarray) – Square matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.
maxiter (int | None) – Maximum number of iterations. Defaults to None (scipy default).

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.

linearsystems.lsigmr(a, b, tol=1e-10, maxiter=None)¶

Solve a linear system using GMRES (Generalized Minimal Residual).

Wraps scipy.sparse.linalg.gmres().

Parameters:

a (np.ndarray) – Square matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.
maxiter (int | None) – Maximum number of restarts. Defaults to None (scipy default).

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.

linearsystems.lsiccg(a, b, tol=1e-10, maxiter=None)¶

Solve a symmetric positive definite system using CG with IC preconditioner.

Computes an incomplete Cholesky factorisation via scipy.sparse.linalg.spilu and uses it as a left preconditioner for scipy.sparse.linalg.cg().

Parameters:

a (np.ndarray) – Symmetric positive definite matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.
maxiter (int | None) – Maximum number of iterations. Defaults to None (scipy default).

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square.

linearsystems.lsisor(a, b, omega=1.5, tol=1e-10, maxiter=1000)¶

Solve a linear system using Successive Over-Relaxation (SOR).

Parameters:

a (np.ndarray) – Square matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
omega (float) – Relaxation factor (0 < omega < 2). Defaults to 1.5.
tol (float) – Convergence tolerance on the residual norm. Defaults to 1e-10.
maxiter (int) – Maximum number of iterations. Defaults to 1000.

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Raises:

ValueError – If a is not square or omega is out of range.

Least Squares¶

linearsystems.lsqrr(a, b)¶

Solve the least-squares problem min‖Ax − b‖₂.

Wraps numpy.linalg.lstsq().

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
b (np.ndarray) – Right-hand side vector of shape (m,) or matrix of shape (m, k).

Returns:

(x, residuals, rank, sv) matching the numpy.linalg.lstsq(): return signature. x is the least-squares solution of shape (n,), residuals are the squared residual norms, rank is the effective rank, and sv are the singular values.

Return type:

tuple

linearsystems.lsqnr(a, b)¶

Solve the non-negative least-squares problem min‖Ax − b‖₂ subject to x ≥ 0.

Wraps scipy.optimize.nnls().

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
b (np.ndarray) – Right-hand side vector of shape (m,).

Returns:

(x, residual_norm) where x is the non-negative solution: of shape (n,) and residual_norm is ‖Ax − b‖₂.

Return type:

tuple

linearsystems.lsqcr(a, b, constraints=None)¶

Solve a (possibly constrained) least-squares problem.

Without constraints delegates to scipy.linalg.lstsq(). With equality constraints C @ x = d the problem is solved via the augmented normal equations.

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
b (np.ndarray) – Right-hand side vector of shape (m,).
constraints (dict | None) – Optional equality constraints. Expected keys are "C" (np.ndarray, shape (p, n)) and "d" (np.ndarray, shape (p,)). Defaults to None.

Returns:

(x, residuals, rank, sv) matching the: scipy.linalg.lstsq() return signature.

Return type:

tuple

linearsystems.lsbrr(a, b, lb=None, ub=None)¶

Solve a bounded least-squares problem using SLSQP.

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
b (np.ndarray) – Right-hand side vector of shape (m,).
lb (np.ndarray | None) – Lower bounds of shape (n,). Defaults to -inf for each component.
ub (np.ndarray | None) – Upper bounds of shape (n,). Defaults to +inf for each component.

Returns:

(x, residual_norm) where x is the bounded: least-squares solution of shape (n,) and residual_norm is ‖Ax − b‖₂.

Return type:

tuple

linearsystems.ridge_regression(a, b, alpha=1.0)¶

Solve a Tikhonov-regularised (ridge) least-squares problem.

Minimises ‖Ax − b‖₂² + alpha * ‖x‖₂².

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
b (np.ndarray) – Right-hand side vector of shape (m,).
alpha (float) – Regularisation parameter (must be positive). Defaults to 1.0.

Returns:

Ridge solution vector of shape (n,).

Return type:

np.ndarray

Raises:

ValueError – If alpha is not positive.

Matrix Inverse and Utility¶

linearsystems.matrix_inverse(a)¶

Compute the inverse of a square matrix.

Wraps numpy.linalg.inv().

Parameters:

a (np.ndarray) – Square real matrix of shape (n, n).

Returns:

Inverse matrix of shape (n, n).

Return type:

np.ndarray

Raises:

ValueError – If a is not square.
numpy.linalg.LinAlgError – If a is singular.

linearsystems.pseudo_inverse(a, rcond=None)¶

Compute the Moore-Penrose pseudoinverse of a matrix.

Wraps numpy.linalg.pinv().

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
rcond (float | None) – Cutoff for small singular values. Values smaller than rcond * largest_sv are set to zero. Defaults to None (numpy default).

Returns:

Pseudoinverse of shape (n, m).

Return type:

np.ndarray

linearsystems.matrix_power(a, n)¶

Raise a square matrix to an integer power.

Wraps numpy.linalg.matrix_power().

Parameters:

a (np.ndarray) – Square real matrix of shape (k, k).
n (int) – Integer exponent (may be negative).

Returns:

Matrix a**n of shape (k, k).

Return type:

np.ndarray

Raises:

ValueError – If a is not square.

linearsystems.solve_upper_triangular(a, b)¶

Solve an upper triangular system by back substitution.

Parameters:

a (np.ndarray) – Upper triangular matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution vector of shape (n,).

Return type:

np.ndarray

Raises:

ValueError – If a is not square.

linearsystems.solve_lower_triangular(a, b)¶

Solve a lower triangular system by forward substitution.

Parameters:

a (np.ndarray) – Lower triangular matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution vector of shape (n,).

Return type:

np.ndarray

Raises:

ValueError – If a is not square.

Norms and Conditioning¶

linearsystems.vector_norm(x, ord=2)¶

Compute the norm of a vector.

Wraps numpy.linalg.norm().

Parameters:

x (np.ndarray) – Input vector of shape (n,).
ord (int | float) – Order of the norm (e.g., 1, 2, np.inf). Defaults to 2.

Returns:

Norm value.

Return type:

float

linearsystems.matrix_norm(a, ord=None)¶

Compute the norm of a matrix.

Wraps numpy.linalg.norm().

Parameters:

a (np.ndarray) – Matrix of shape (m, n).
ord (int | float | str | None) – Order of the norm. Supported values: None (Frobenius), 'fro', 1, -1, 2, -2, np.inf, -np.inf. Defaults to None.

Returns:

Norm value.

Return type:

float

linearsystems.condition_number(a, ord=None)¶

Compute the condition number of a matrix.

Wraps numpy.linalg.cond().

Parameters:

a (np.ndarray) – Square matrix of shape (n, n).
ord (int | float | str | None) – Norm order used to compute the condition number. Defaults to None (2-norm).

Returns:

Condition number.

Return type:

float

linearsystems.rcond_estimate(a)¶

Estimate the reciprocal of the condition number of a matrix.

Uses the ratio of the smallest to the largest singular value.

Parameters:

a (np.ndarray) – Matrix of shape (m, n).

Returns:

Reciprocal condition number estimate in [0, 1]. A value: close to 0 indicates near-singularity.

Return type:

float

Banded Matrix Utilities¶

linearsystems.to_banded(a, lower, upper)¶

Convert a full matrix to banded (LAPACK) storage.

The output has shape (lower + upper + 1, n) where ab[u + i - j, j] stores a[i, j].

Parameters:

a (np.ndarray) – Dense matrix of shape (n, n).
lower (int) – Number of sub-diagonals to store.
upper (int) – Number of super-diagonals to store.

Returns:

Banded storage array of shape: (lower + upper + 1, n).

Return type:

np.ndarray

Raises:

ValueError – If a is not square.

linearsystems.from_banded(ab, lower, upper, n)¶

Convert banded (LAPACK) storage back to a full matrix.

Parameters:

ab (np.ndarray) – Banded storage array of shape (lower + upper + 1, n).
lower (int) – Number of sub-diagonals stored.
upper (int) – Number of super-diagonals stored.
n (int) – Dimension of the output square matrix.

Returns:

Dense matrix of shape (n, n).

Return type:

np.ndarray

linearsystems.lslbd(ab, b, lower, upper)¶

Solve a banded linear system from banded storage.

Wraps scipy.linalg.solve_banded().

Parameters:

ab (np.ndarray) – Banded storage array of shape (lower + upper + 1, n) as used by LAPACK.
b (np.ndarray) – Right-hand side vector of shape (n,).
lower (int) – Number of sub-diagonals.
upper (int) – Number of super-diagonals.

Returns:

Solution with x, residual, success,: and message fields.

Return type:

LinearSystemResult

Sparse System Utilities¶

linearsystems.sparse_solve(a_sparse, b)¶

Solve a sparse linear system Ax = b.

Wraps scipy.sparse.linalg.spsolve().

Parameters:

a_sparse (scipy.sparse.spmatrix) – Sparse square matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution vector of shape (n,).

Return type:

np.ndarray

linearsystems.sparse_cg(a_sparse, b, tol=1e-10)¶

Solve a sparse SPD system using Conjugate Gradient.

Wraps scipy.sparse.linalg.cg().

Parameters:

a_sparse (scipy.sparse.spmatrix) – Sparse symmetric positive definite matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

linearsystems.sparse_gmres(a_sparse, b, tol=1e-10)¶

Solve a sparse linear system using GMRES.

Wraps scipy.sparse.linalg.gmres().

Parameters:

a_sparse (scipy.sparse.spmatrix) – Sparse square matrix of shape (n, n).
b (np.ndarray) – Right-hand side vector of shape (n,).
tol (float) – Convergence tolerance. Defaults to 1e-10.

Returns:

Solution with x, residual, success,: message, and n_iter fields.

Return type:

LinearSystemResult

Tridiagonal Solvers¶

linearsystems.tridiag_solve(lower, main, upper, b)¶

Solve a tridiagonal linear system using the Thomas algorithm.

Parameters:

lower (np.ndarray) – Sub-diagonal of length n-1.
main (np.ndarray) – Main diagonal of length n.
upper (np.ndarray) – Super-diagonal of length n-1.
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution vector of shape (n,).

Return type:

np.ndarray

Raises:

ValueError – If the diagonal lengths are inconsistent.

linearsystems.cyclic_tridiag_solve(lower, main, upper, b)¶

Solve a cyclic tridiagonal system using the Sherman-Morrison formula.

The cyclic system has extra corner entries: A[0, n-1] = lower[-1] and A[n-1, 0] = upper[-1].

Parameters:

lower (np.ndarray) – Sub-diagonal entries including the wrap-around corner; length n.
main (np.ndarray) – Main diagonal of length n.
upper (np.ndarray) – Super-diagonal entries including the wrap-around corner; length n.
b (np.ndarray) – Right-hand side vector of shape (n,).

Returns:

Solution vector of shape (n,).

Return type:

np.ndarray

Raises:

ValueError – If the diagonal lengths are inconsistent.

linearsystems.block_tridiag_solve(blocks, b)¶

Solve a block tridiagonal system via block-LU elimination (numpy).

Each element of blocks is a triple (L_i, D_i, U_i) where D_i is the diagonal block and L_i, U_i are the off-diagonal blocks below and above the diagonal respectively. The first block has no L_0 (pass None) and the last block has no U_last (pass None).

Parameters:

blocks (list[tuple[np.ndarray, np.ndarray, np.ndarray]]) – Block triplets (L_i, D_i, U_i) of length N (number of blocks). L_0 and U_{N-1} should be None.
b (np.ndarray) – Right-hand side vector of shape (total_rows,) where total_rows is the sum of all diagonal block row counts.

Returns:

Solution vector of shape (total_rows,).

Return type:

np.ndarray