mathematicskit.linalg#

mathematicskit.linalg hub · History · Examples · API reference

LU decomposition with partial pivoting (scipy.linalg.lu); QR via Householder reflections (scipy.linalg.qr) and a hand-rolled Gram-Schmidt comparison; Cholesky decomposition (numpy.linalg.cholesky); eigenvalue computation (numpy.linalg.eigh/eig) alongside hand-rolled power/inverse iteration; SVD (numpy.linalg.svd); conjugate gradient and GMRES (scipy.sparse.linalg.cg/gmres); condition-number estimation (numpy.linalg.cond) and a least-squares stability comparison (normal equations vs. QR vs. numpy.linalg.lstsq).

mathematicskit.linalg: numerical linear algebra, built on numpy.linalg/scipy.linalg/scipy.sparse.linalg.

LU decomposition with partial pivoting (scipy.linalg.lu()); QR decomposition via Householder reflections (scipy.linalg.qr()), with a hand-rolled Gram-Schmidt variant kept only for the pedagogical Householder-vs-Gram-Schmidt stability comparison; Cholesky decomposition for SPD matrices (numpy.linalg.cholesky()); eigenvalue computation (numpy.linalg.eigh()/numpy.linalg.eig()); power iteration and inverse iteration for dominant/nearest eigenvalues (hand-rolled – the iteration itself is the pedagogical subject, with no library equivalent); SVD (numpy.linalg.svd()); conjugate gradient and GMRES (scipy.sparse.linalg.cg()/gmres()); condition number (numpy.linalg.cond()) and least-squares (numpy.linalg.lstsq()), with a comparison utility showing why normal equations are less stable than QR for ill-conditioned systems. Also: Cramer’s rule (numpy.linalg.slogdet()); characteristic polynomials (numpy.poly()) and the Cayley-Hamilton theorem; Jacobi, Gauss-Seidel, and SOR stationary iterations (hand-rolled – no library equivalent); the Moore-Penrose pseudoinverse (numpy.linalg.pinv()); Gershgorin discs; Lanczos extreme eigenpairs (scipy.sparse.linalg.eigsh()); and Hessenberg/Schur forms (scipy.linalg.hessenberg()/schur()).

mathematicskit’s systems/ classes/functions wrap these library calls in the package’s own dataclass results, docstrings, and stability commentary – they are thin, well-tested wrappers, not reimplementations of already- correct numerical primitives.

class mathematicskit.linalg.CholeskyResult(L)[source]#

Bases: object

Output of a Cholesky decomposition of an SPD matrix: A = L L^T.

Parameters:

L (ndarray)

L: ndarray#

Lower-triangular factor.

Type:

ndarray, shape (n, n)

class mathematicskit.linalg.ConjugateGradient(tol=1e-08, max_iter=1000)[source]#

Bases: IterativeLinearSolver

Conjugate gradient method for symmetric positive-definite systems.

Minimizes the quadratic \(\phi(x) = \tfrac12 x^T A x - b^T x\) (equivalent to solving A x = b for SPD A) by generating search directions \(A\)-conjugate to all previous ones, which guarantees exact convergence in at most n iterations in exact arithmetic. Calls scipy.sparse.linalg.cg() directly. See Golub & Van Loan, Matrix Computations, 4th ed., Ch. 11.3, Algorithm 11.3.1.

Parameters:

  • tol (float) – Relative residual-norm tolerance, ||r_k|| < tol * ||b|| (passed as rtol to scipy.sparse.linalg.cg()).

  • max_iter (int) – Maximum number of iterations (defaults to DEFAULT_MAX_ITER; n is normally sufficient).

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> result = ConjugateGradient().solve(A, b)
>>> np.allclose(A @ result.x, b, atol=1e-8)
True
solve(a, b, x0=None)[source]#

Solve A x = b iteratively.

Parameters:
Return type:

IterativeSolveResult

Returns:

IterativeSolveResult

class mathematicskit.linalg.EigenResult(eigenvalues, eigenvectors, iterations=0, converged=True, method='', extra=<factory>)[source]#

Bases: object

Output of an eigenvalue computation, whole-spectrum or single-pair.

The whole-spectrum solvers (eigen_symmetric(), eigen_general()) fill eigenvalues/eigenvectors with the full spectrum in one shot and report iterations=1, converged=True, since LAPACK does not expose its internal iteration count. The hand-rolled iterations (power_iteration(), inverse_iteration()) return a single eigenpair and a genuine iteration count and convergence flag.

Parameters:
converged: bool = True#

Whether the convergence tolerance was met (always True for the direct library solvers, which do not iterate visibly).

Type:

bool

eigenvalues: ndarray#

The whole spectrum, shape (n,), for the library-backed solvers; a single dominant (or shift-selected) eigenvalue, shape (1,), for power/inverse iteration. Real and ascending from eigh, possibly complex and unordered from eig.

Type:

ndarray

eigenvectors: ndarray#

Shape (n, n) with eigenvectors as columns for the whole-spectrum solvers, so eigenvectors[:, k] pairs with eigenvalues[k]; a single unit eigenvector of shape (n,) for power/inverse iteration.

Type:

ndarray

extra: dict#

Free-form diagnostics slot, unused by the current solvers.

Type:

dict

iterations: int = 0#

Number of iterations performed (always 1 for the library-backed solvers, which do not expose LAPACK’s internal count).

Type:

int

method: str = ''#

One of "numpy_eigh", "numpy_eig", "power_iteration", or "inverse_iteration".

Type:

str

class mathematicskit.linalg.GMRES(tol=1e-08, max_iter=1000)[source]#

Bases: IterativeLinearSolver

GMRES: Generalized Minimal RESidual method for general (non-symmetric) systems.

Builds an orthonormal basis of the Krylov subspace via the Arnoldi process and picks the iterate minimizing the residual norm over it. Unlike conjugate gradient, GMRES applies to any nonsingular A (not just SPD). Calls scipy.sparse.linalg.gmres() directly with restart capped at n (so a single restart cycle already spans the full Krylov space, i.e. “full GMRES”, matching the textbook presentation). See Saad & Schultz (1986), GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comput. 7(3), and Golub & Van Loan, Matrix Computations, 4th ed., Ch. 9.3.3.

Parameters:

  • tol (float) – Relative residual-norm tolerance.

  • max_iter (int) – Maximum Krylov subspace dimension.

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 2.0, 5.0]])
>>> b = np.array([1.0, 2.0, 3.0])
>>> result = GMRES().solve(A, b)
>>> np.allclose(A @ result.x, b, atol=1e-8)
True
solve(a, b, x0=None)[source]#

Solve A x = b iteratively.

Parameters:
Return type:

IterativeSolveResult

Returns:

IterativeSolveResult

class mathematicskit.linalg.GaussSeidel(tol=1e-08, max_iter=1000)[source]#

Bases: SOR

Gauss-Seidel iteration (Gauss 1823, Seidel 1874): SOR with omega = 1.

Each unknown is solved from its own equation using the newest values of the unknowns already updated in this sweep. Converges for strictly diagonally dominant or symmetric positive-definite A.

Parameters:

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> result = GaussSeidel(tol=1e-12).solve(A, b)
>>> np.round(result.x, 8)
array([0.09090909, 0.63636364])
method = 'gauss_seidel'#
class mathematicskit.linalg.GershgorinResult(centers, radii)[source]#

Bases: object

The Gershgorin discs of a square matrix.

Disc \(k\) is centered at \(a_{kk}\) with radius \(R_k = \sum_{j \ne k} |a_{kj}|\); every eigenvalue lies in the union of the discs.

Parameters:
centers: ndarray#

Disc centers, the diagonal entries (possibly complex).

Type:

ndarray, shape (n,)

contains(z)[source]#

Whether the complex number z lies in the union of the discs (boundary included, up to rounding).

Return type:

bool

Parameters:

z (complex)

radii: ndarray#

Disc radii, the off-diagonal absolute row (or column) sums.

Type:

ndarray, shape (n,)

class mathematicskit.linalg.IterativeLinearSolver(tol=1e-08, max_iter=1000)[source]#

Bases: ABC

Common base for Krylov-subspace iterative linear-system solvers.

Parameters:

  • tol (float) – Relative residual-norm convergence tolerance.

  • max_iter (int) – Maximum number of iterations.

abstractmethod solve(a, b, x0=None)[source]#

Solve A x = b iteratively.

Parameters:
Return type:

IterativeSolveResult

Returns:

IterativeSolveResult

class mathematicskit.linalg.IterativeSolveResult(x, residual_history, iterations=0, converged=True, method='')[source]#

Bases: object

Output of an iterative linear-system solver (CG, GMRES).

Parameters:
converged: bool = True#

Whether the residual tolerance was met.

Type:

bool

iterations: int = 0#

Number of iterations performed.

Type:

int

method: str = ''#

"conjugate_gradient", "gmres", "jacobi", "gauss_seidel", or "sor".

Type:

str

residual_history: ndarray#

The relative residual norm ||b - A x_k|| / ||b|| at each iterate, starting from the initial guess. Normalizing by ||b|| is what makes the history directly comparable to method’s tol (itself a relative tolerance) and comparable across right-hand sides of different magnitudes.

Type:

ndarray, shape (iterations + 1,)

x: ndarray#

Approximate solution.

Type:

ndarray, shape (n,)

class mathematicskit.linalg.JacobiIteration(tol=1e-08, max_iter=1000)[source]#

Bases: SOR

Jacobi iteration (Jacobi 1845): every unknown updated simultaneously from the previous iterate.

\(x_{k+1} = D^{-1}(b - (L + U) x_k)\). Simpler and fully parallel, but for consistently ordered matrices it needs about twice as many sweeps as Gauss-Seidel.

Parameters:

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> result = JacobiIteration(tol=1e-12).solve(A, b)
>>> np.allclose(A @ result.x, b)
True
method = 'jacobi'#
class mathematicskit.linalg.LUResult(L, U, P, num_row_swaps=0)[source]#

Bases: object

Output of an LU decomposition with partial pivoting: P A = L U.

Parameters:
L: ndarray#

Unit lower-triangular factor.

Type:

ndarray, shape (n, n)

P: ndarray#

Row-permutation matrix, P @ A == L @ U.

Type:

ndarray, shape (n, n)

U: ndarray#

Upper-triangular factor.

Type:

ndarray, shape (n, n)

num_row_swaps: int = 0#

Number of row swaps performed (parity gives det(P)).

Type:

int

class mathematicskit.linalg.LeastSquaresResult(coefficients, residual_norm=0.0, method='', condition_number=None)[source]#

Bases: object

Output of a least-squares solve (normal equations or QR-based).

Parameters:
coefficients: ndarray#

Least-squares solution.

Type:

ndarray, shape (n,)

condition_number: float | None = None#

2-norm condition number of A ("qr") or of the normal-equations matrix A^T A ("normal_equations", which squares it).

Type:

float, optional

method: str = ''#

"normal_equations" or "qr".

Type:

str

residual_norm: float = 0.0#

||b - A x||_2 at the solution.

Type:

float

class mathematicskit.linalg.QRResult(Q, R, method='')[source]#

Bases: object

Output of a QR decomposition: A = Q R.

Parameters:
Q: ndarray#

Orthonormal-columns factor.

Type:

ndarray, shape (m, n)

R: ndarray#

Upper-triangular factor.

Type:

ndarray, shape (n, n)

method: str = ''#

"householder" or "gram_schmidt".

Type:

str

class mathematicskit.linalg.SOR(omega=1.0, tol=1e-08, max_iter=1000)[source]#

Bases: IterativeLinearSolver

Successive over-relaxation (Frankel 1950, Young 1950).

One sweep updates the unknowns in order, blending the Gauss-Seidel value with the old one:

\[x_i \leftarrow (1 - \omega)\, x_i + \frac{\omega}{a_{ii}} \Big(b_i - \sum_{j < i} a_{ij} x_j^{\text{new}} - \sum_{j > i} a_{ij} x_j^{\text{old}}\Big).\]

omega = 1 is Gauss-Seidel; 1 < omega < 2 over-relaxes. For symmetric positive-definite A it converges for every 0 < omega < 2 (Ostrowski-Reich).

Parameters:

  • omega (float) – Relaxation factor, 0 < omega < 2.

  • tol (float) – Relative residual tolerance, ||b - A x_k|| < tol * ||b||.

  • max_iter (int) – Maximum number of sweeps.

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> result = SOR(omega=1.1, tol=1e-12).solve(A, b)
>>> np.allclose(A @ result.x, b), result.method
(True, 'sor')
method = 'sor'#
solve(a, b, x0=None)[source]#

Solve A x = b iteratively.

Parameters:
Return type:

IterativeSolveResult

Returns:

IterativeSolveResult

class mathematicskit.linalg.SVDResult(U, S, Vt)[source]#

Bases: object

Output of a singular value decomposition: A = U diag(S) V^T.

Parameters:
S: ndarray#

Singular values, descending.

Type:

ndarray, shape (k,)

U: ndarray#

Left singular vectors.

Type:

ndarray, shape (m, k)

Vt: ndarray#

Right singular vectors, transposed.

Type:

ndarray, shape (k, n)

class mathematicskit.linalg.SchurResult(T, Z)[source]#

Bases: object

Output of a Schur decomposition: A = Z T Z^H.

Parameters:
T: ndarray#

Schur form – upper triangular (complex output) or quasi-upper-triangular with 2x2 blocks for complex-conjugate eigenvalue pairs (real output). Its diagonal (or block) entries are the eigenvalues of A.

Type:

ndarray, shape (n, n)

Z: ndarray#

Unitary (orthogonal, for real output) Schur vectors.

Type:

ndarray, shape (n, n)

property eigenvalues: ndarray#

Eigenvalues read off T (2x2 diagonal blocks resolved with numpy.linalg.eigvals()).

Type:

ndarray, shape (n,)

mathematicskit.linalg.characteristic_polynomial(a)[source]#

Coefficients of \(p_A(\lambda) = \det(\lambda I - A)\), highest degree first, via numpy.poly().

The polynomial is monic, the coefficient of \(\lambda^{n-1}\) is \(-\operatorname{tr} A\), and the constant term is \((-1)^n \det A\).

Parameters:

a (ndarray)

Return type:

ndarray

Returns:

ndarray, shape (n + 1,) – Real if a is real (imaginary round-off from complex-conjugate eigenvalue pairs is discarded).

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> np.round(characteristic_polynomial(A), 8)  # l^2 - 5 l - 2
array([ 1., -5., -2.])
mathematicskit.linalg.cholesky_decompose(a)[source]#

Factor a symmetric positive-definite A = L L^T via numpy.linalg.cholesky().

Roughly half the cost of a general LU factorization (exploits symmetry); the underlying LAPACK routine (?potrf) raises when a diagonal pivot would require a square root of a negative number, which doubles as a certificate that a is not positive-definite. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 6.6, Theorem 6.25.

Parameters:

a (ndarray) – Symmetric positive-definite matrix.

Return type:

CholeskyResult

Returns:

CholeskyResult

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 2.0], [2.0, 5.0]])
>>> result = cholesky_decompose(A)
>>> np.allclose(result.L @ result.L.T, A)
True
mathematicskit.linalg.cholesky_solve(result, b)[source]#

Solve A x = b given a precomputed CholeskyResult.

Forward substitution with L, then back substitution with L^T, via scipy.linalg.solve_triangular().

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 2.0], [2.0, 5.0]])
>>> b = np.array([2.0, 1.0])
>>> x = cholesky_solve(cholesky_decompose(A), b)
>>> np.allclose(A @ x, b)
True
mathematicskit.linalg.condition_number_2norm(a)[source]#

2-norm condition number \(\kappa_2(A) = \sigma_{\max}/\sigma_{\min}\), via numpy.linalg.cond().

A large \(\kappa_2\) means small relative perturbations in A or b can produce large relative changes in the solution of A x = b – see least_squares_normal_equations() for why this matters concretely for least squares. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 12.

Parameters:

a (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> round(condition_number_2norm(np.diag([10.0, 1.0])), 6)
10.0
mathematicskit.linalg.cramer_solve(a, b)[source]#

Solve A x = b by Cramer’s rule, \(x_i = \det(A_i) / \det(A)\).

\(A_i\) is A with its i-th column replaced by b.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

Raises:

numpy.linalg.LinAlgError – If det(A) is (numerically) zero.

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 1.0], [1.0, 3.0]])
>>> b = np.array([3.0, 5.0])
>>> np.round(cramer_solve(A, b), 8)
array([0.8, 1.4])
mathematicskit.linalg.eigen_general(a)[source]#

Eigenvalues/eigenvectors of a general (possibly non-symmetric) matrix, via numpy.linalg.eig().

Eigenvalues/eigenvectors may be complex even for a real input matrix (e.g. a rotation matrix). See Golub & Van Loan, Matrix Computations, 4th ed., Ch. 7.

Parameters:

a (ndarray)

Return type:

EigenResult

Returns:

EigenResult – method="numpy_eig".

Examples

>>> import numpy as np
>>> A = np.array([[0.0, -1.0], [1.0, 0.0]])  # 90-degree rotation
>>> result = eigen_general(A)
>>> np.round(np.sort(np.abs(result.eigenvalues.imag)), 8)
array([1., 1.])
mathematicskit.linalg.eigen_symmetric(a)[source]#

Eigenvalues/eigenvectors of a symmetric (Hermitian) matrix, via numpy.linalg.eigh().

eigh exploits symmetry to guarantee real eigenvalues and an orthonormal eigenvector basis, returned in ascending order, at roughly half the cost of the general eigenvalue problem. See Golub & Van Loan, Matrix Computations, 4th ed., Ch. 8.

Parameters:

a (ndarray) – Symmetric matrix.

Return type:

EigenResult

Returns:

EigenResult – eigenvalues shape (n,) ascending, eigenvectors shape (n, n) with columns the corresponding eigenvectors, method="numpy_eigh".

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 1.0], [1.0, 2.0]])
>>> result = eigen_symmetric(A)
>>> np.round(result.eigenvalues, 8)
array([1., 3.])
mathematicskit.linalg.frobenius_norm(a)[source]#

Frobenius norm \(\|A\|_F = \sqrt{\sum_{ij} A_{ij}^2}\), via numpy.linalg.norm().

Parameters:

a (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> round(frobenius_norm(np.array([[3.0, 0.0], [0.0, 4.0]])), 6)
5.0
mathematicskit.linalg.gershgorin_discs(a, by='rows')[source]#

The Gershgorin discs of a square matrix.

Disc \(k\) has center \(a_{kk}\) and radius \(R_k = \sum_{j \ne k} |a_{kj}|\) (by="rows") or \(\sum_{j \ne k} |a_{jk}|\) (by="columns"). A consequence: a strictly diagonally dominant matrix (every \(|a_{kk}| > R_k\)) has no eigenvalue at 0 and is therefore nonsingular.

Parameters:

  • a (ndarray) – Real or complex square matrix.

  • by (str)

Return type:

GershgorinResult

Returns:

GershgorinResult

Examples

>>> import numpy as np
>>> A = np.array([[10.0, 1.0, 0.0], [0.2, 8.0, 0.2], [1.0, 1.0, 2.0]])
>>> discs = gershgorin_discs(A)
>>> discs.centers, discs.radii
(array([10.,  8.,  2.]), array([1. , 0.4, 2. ]))
>>> all(discs.contains(lam) for lam in np.linalg.eigvals(A))
True
mathematicskit.linalg.gram_schmidt_qr(a, modified=True)[source]#

QR decomposition by (modified or classical) Gram-Schmidt.

Orthogonalizes A’s columns one at a time against the previously-computed Q columns. Modified Gram-Schmidt (modified=True, the default) subtracts each projection immediately as it’s computed rather than all at once from the original vector (classical Gram-Schmidt, modified=False) – algebraically equivalent in exact arithmetic, but modified Gram-Schmidt loses far less orthogonality in floating point, especially for near-collinear columns. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 8, and orthogonality_error() for a direct numerical comparison.

Parameters:

  • a (ndarray) – Full-column-rank matrix.

  • modified (bool) – Whether to use modified (default) or classical Gram-Schmidt.

Return type:

QRResult

Returns:

QRResult – With method="gram_schmidt".

Examples

>>> import numpy as np
>>> A = np.array([[1.0, -1.0], [1.0, 1.0], [0.0, 1.0]])
>>> result = gram_schmidt_qr(A)
>>> np.allclose(result.Q @ result.R, A)
True
mathematicskit.linalg.hessenberg_reduce(a)[source]#

Upper-Hessenberg form A = Q H Q^T via scipy.linalg.hessenberg().

H is zero below its first subdiagonal. Reducing to it first (by Householder reflections) cuts each QR step from \(O(n^3)\) to \(O(n^2)\), and the QR iteration preserves the Hessenberg shape.

Parameters:

a (ndarray)

Return type:

tuple[ndarray, ndarray]

Returns:

  • H (ndarray, shape (n, n))

  • Q (ndarray, shape (n, n)) – Orthogonal (unitary) similarity.

Examples

>>> import numpy as np
>>> A = np.arange(16.0).reshape(4, 4) + np.eye(4)
>>> H, Q = hessenberg_reduce(A)
>>> bool(np.allclose(np.tril(H, -2), 0.0)), bool(np.allclose(Q @ H @ Q.T, A))
(True, True)
mathematicskit.linalg.householder_qr(a)[source]#

QR decomposition by Householder reflections, via scipy.linalg.qr().

Each column below the diagonal is zeroed by an orthogonal reflector \(H_k = I - 2 v_k v_k^T / (v_k^T v_k)\) chosen to map that column’s trailing part onto a multiple of \(e_1\); the product of all the \(H_k\) is Q. mode="economic" requests the thin factorization (Q shape (m, n) rather than (m, m)) to match gram_schmidt_qr()’s shape. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 10.

Parameters:

a (ndarray) – Full-column-rank matrix.

Return type:

QRResult

Returns:

QRResult – With method="householder", Q shape (m, n), R shape (n, n).

Examples

>>> import numpy as np
>>> A = np.array([[1.0, -1.0], [1.0, 1.0], [0.0, 1.0]])
>>> result = householder_qr(A)
>>> np.allclose(result.Q @ result.R, A)
True
>>> np.allclose(result.Q.T @ result.Q, np.eye(2), atol=1e-10)
True
mathematicskit.linalg.inverse_iteration(a, mu, tol=1e-12, max_iter=1000, v0=None)[source]#

Inverse iteration for the eigenvalue nearest a shift mu.

Power iteration applied to \((A - \mu I)^{-1}\): since that matrix’s dominant eigenvalue is \(1/(\lambda - \mu)\) for the \(\lambda\) closest to mu, iterating it converges to the corresponding eigenvector, with the Rayleigh quotient of the original A giving that eigenvalue directly. The linear solve at each step uses numpy.linalg.solve() rather than an explicit inverse; only the iteration itself is hand-rolled. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 9.3, and Golub & Van Loan, Matrix Computations, 4th ed., Ch. 8.2.2.

Parameters:
Return type:

EigenResult

Returns:

EigenResult

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 0.0], [0.0, 5.0]])
>>> result = inverse_iteration(A, mu=4.5)
>>> round(float(result.eigenvalues[0]), 6)
5.0
mathematicskit.linalg.is_symmetric(a, atol=1e-10)[source]#

Check whether a matrix is (numerically) symmetric.

Parameters:
Return type:

bool

Returns:

bool

Examples

>>> import numpy as np
>>> is_symmetric(np.array([[1.0, 2.0], [2.0, 1.0]]))
True
mathematicskit.linalg.is_symmetric_positive_definite(a)[source]#

Check symmetric positive-definiteness by attempting a Cholesky factorization.

Parameters:

a (ndarray)

Return type:

bool

Returns:

bool

Examples

>>> import numpy as np
>>> is_symmetric_positive_definite(np.array([[4.0, 2.0], [2.0, 5.0]]))
True
>>> is_symmetric_positive_definite(np.array([[1.0, 2.0], [2.0, 1.0]]))
False
mathematicskit.linalg.jacobi_spectral_radius(a)[source]#

Spectral radius \(\rho(G_J)\) of the Jacobi iteration matrix \(G_J = I - D^{-1} A\).

Computed with numpy.linalg.eigvals(). The Jacobi iteration converges from every starting vector iff this is below 1.

Parameters:

a (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> n = 5
>>> A = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
>>> bool(np.isclose(jacobi_spectral_radius(A), np.cos(np.pi / (n + 1))))
True
mathematicskit.linalg.lanczos_eigsh(a, k=6, which='LA', tol=0.0, sigma=None)[source]#

The k extreme eigenpairs of a symmetric matrix, via scipy.sparse.linalg.eigsh() (Lanczos/ARPACK).

Lanczos builds an orthonormal basis \(Q_m\) of the Krylov subspace \(\operatorname{span}\{v, Av, \dots, A^{m-1}v\}\) with a three-term recurrence, in which \(T_m = Q_m^T A Q_m\) is tridiagonal; the eigenvalues of \(T_m\) (Ritz values) converge first to the extreme eigenvalues of A.

Parameters:

  • a (ndarray, sparse matrix, or LinearOperator, shape (n, n)) – Symmetric (Hermitian) matrix.

  • k (int) – Number of eigenpairs, k < n.

  • which (str) – Largest/smallest algebraic, largest/smallest magnitude, or both ends of the spectrum.

  • tol (float) – Relative accuracy for the Ritz values; 0 means machine precision.

  • sigma (float | None) – Shift-invert mode: run Lanczos on \((A - \sigma I)^{-1}\) (factored once by sparse LU), which returns the k eigenvalues nearest sigma when which="LM". The standard route to interior or tightly clustered eigenvalues, which plain Lanczos resolves only slowly.

Return type:

EigenResult

Returns:

EigenResult – eigenvalues shape (k,) ascending, eigenvectors shape (n, k), method="lanczos_arpack".

Examples

>>> import numpy as np
>>> A = np.diag(np.arange(1.0, 21.0))
>>> np.round(lanczos_eigsh(A, k=3).eigenvalues, 8)
array([18., 19., 20.])
>>> np.round(lanczos_eigsh(A, k=2, sigma=7.2, which="LM").eigenvalues, 8)
array([7., 8.])
mathematicskit.linalg.least_squares_lstsq(a, b)[source]#

Least squares via numpy.linalg.lstsq() (SVD-based, via LAPACK ?gelsd).

The most robust of the three routes here – it handles rank-deficient a gracefully (returning the minimum-norm solution) where both least_squares_normal_equations() and least_squares_qr() would fail or become ill-conditioned – because it never needs a to have full column rank in the first place. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 11.

Parameters:
Return type:

LeastSquaresResult

Returns:

LeastSquaresResult – method="lstsq".

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
>>> b = np.array([2.0, 3.0, 5.0])
>>> result = least_squares_lstsq(A, b)
>>> np.round(result.coefficients, 4)
array([0.3333, 1.5   ])
mathematicskit.linalg.least_squares_normal_equations(a, b)[source]#

Least squares via the normal equations \(A^T A x = A^T b\).

A^T A and A^T b are formed explicitly and solved with numpy.linalg.solve(). Simple, but \(\kappa_2(A^T A) = \kappa_2(A)^2\) – squaring the condition number of A itself – so this can lose roughly twice as many digits of accuracy as the QR-based approach (least_squares_qr()) for an ill-conditioned A. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 11.

Parameters:
Return type:

LeastSquaresResult

Returns:

LeastSquaresResult

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
>>> b = np.array([2.0, 3.0, 5.0])
>>> result = least_squares_normal_equations(A, b)
>>> np.round(result.coefficients, 4)
array([0.3333, 1.5   ])
mathematicskit.linalg.least_squares_qr(a, b)[source]#

Least squares via A = Q R (scipy.linalg.qr()): solve R x = Q^T b.

Avoids ever forming \(A^T A\), so the effective conditioning is \(\kappa_2(A)\) rather than its square – the standard remedy for least_squares_normal_equations()’s instability. (For very ill-conditioned or rank-deficient problems, numpy.linalg.lstsq() – SVD-based – is even more robust; the QR route is used here specifically to make the comparison with the normal equations concrete.) See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 11.

Parameters:
Return type:

LeastSquaresResult

Returns:

LeastSquaresResult

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
>>> b = np.array([2.0, 3.0, 5.0])
>>> result = least_squares_qr(A, b)
>>> np.round(result.coefficients, 4)
array([0.3333, 1.5   ])
mathematicskit.linalg.lu_decompose(a)[source]#

Factor P A = L U by Gaussian elimination with partial pivoting.

Thin wrapper around scipy.linalg.lu(), which brings the largest-magnitude entry in each column to the pivot position by a row swap before eliminating below it – the standard remedy for the numerical instability of naive (unpivoted) Gaussian elimination. scipy.linalg.lu() returns a permutation matrix P0 with A = P0 L U; this function returns its transpose so that, matching the textbook convention, result.P @ a == result.L @ result.U. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 6.5 (“Matrix Factorizations”), and Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 21.

Parameters:

a (ndarray) – Square, nonsingular matrix.

Return type:

LUResult

Returns:

LUResult

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]])
>>> result = lu_decompose(A)
>>> np.allclose(result.P @ A, result.L @ result.U)
True
>>> # Singularity is judged from U's pivots, not from det(A), which
>>> # underflows to 0.0 for large well-conditioned matrices.
>>> np.allclose(lu_decompose(0.01 * np.eye(200)).U, 0.01 * np.eye(200))
True
mathematicskit.linalg.lu_det(a)[source]#

Determinant via LU decomposition: \(\det A = (-1)^s \prod_i U_{ii}\).

Thin wrapper around scipy.linalg.det(), which itself factors a via LAPACK’s ?getrf and takes the signed product of U’s diagonal – \(O(n^3)\), versus the \(O(n!)\) cofactor expansion. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 6.5.

Parameters:

a (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> round(lu_det(A), 10)
-2.0
mathematicskit.linalg.lu_solve(result, b)[source]#

Solve A x = b given a precomputed LUResult for A.

Solves L y = P b then U x = y by triangular substitution (scipy.linalg.solve_triangular(), \(O(n^2)\)) rather than refactoring at \(O(n^3)\) – the whole point of reusing a factorization across several right-hand sides.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 1.0], [1.0, 3.0]])
>>> b = np.array([3.0, 5.0])
>>> x = lu_solve(lu_decompose(A), b)
>>> np.allclose(A @ x, b)
True
mathematicskit.linalg.lu_solve_system(a, b)[source]#

Convenience: factor a and solve in one call.

Uses scipy.linalg.lu_factor()/scipy.linalg.lu_solve() directly (LAPACK’s ?getrf/?getrs pair) rather than routing through lu_decompose()’s explicit P/L/U, which is marginally faster since it skips reconstructing the dense factors.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> import numpy as np
>>> A = np.array([[3.0, 2.0, -1.0], [2.0, -2.0, 4.0], [-1.0, 0.5, -1.0]])
>>> x_true = np.array([1.0, -2.0, -2.0])
>>> b = A @ x_true
>>> np.allclose(lu_solve_system(A, b), x_true, atol=1e-8)
True
mathematicskit.linalg.matrix_polynomial(coeffs, a)[source]#

Evaluate \(p(A) = c_0 A^d + c_1 A^{d-1} + \dots + c_d I\) by Horner’s rule.

With coeffs = characteristic_polynomial(a) the result is the zero matrix, up to rounding – the Cayley-Hamilton theorem.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n, n)

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> np.allclose(matrix_polynomial(characteristic_polynomial(A), A), 0.0)
True
mathematicskit.linalg.optimal_sor_omega(a)[source]#

Young’s optimal SOR relaxation factor \(\omega^* = 2 / (1 + \sqrt{1 - \rho(G_J)^2})\).

Exact for consistently ordered matrices whose Jacobi iteration matrix has real eigenvalues with \(\rho(G_J) < 1\) (e.g. the discrete Poisson matrices); at \(\omega^*\) the SOR spectral radius drops to \(\omega^* - 1\). See D. M. Young, Trans. AMS 76 (1954), 92-111.

Parameters:

a (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> A = 2 * np.eye(3) - np.eye(3, k=1) - np.eye(3, k=-1)  # rho_J = cos(pi/4)
>>> round(optimal_sor_omega(A), 6)  # 2 / (1 + sqrt(1 - cos(pi/4)^2))
1.171573
mathematicskit.linalg.orthogonality_error(q)[source]#

Measure how far Q’s columns are from exactly orthonormal.

\(\max_{ij}\left|(Q^TQ - I)_{ij}\right|\): the largest entry of Q’s Gram matrix minus the identity. (This is the entrywise max norm, not the induced \(\infty\)-norm, which would be the largest row sum.) Exactly zero for a perfectly orthonormal Q, and a direct numerical demonstration of which QR method loses orthogonality on ill-conditioned input: on a degree-12 Vandermonde matrix, classical Gram-Schmidt reaches order 1, modified Gram-Schmidt stays near 1e-10, and Householder near machine epsilon. See Trefethen & Bau, Numerical Linear Algebra, 1997, Lecture 9 (Figure 9.2’s classical-vs-modified comparison).

Parameters:

q (ndarray)

Return type:

float

Returns:

float

Examples

>>> import numpy as np
>>> orthogonality_error(np.eye(3)) == 0.0
True
mathematicskit.linalg.penrose_residuals(a, x)[source]#

Frobenius-norm residuals of Penrose’s four equations for a candidate X = A^+.

\[AXA = A, \qquad XAX = X, \qquad (AX)^H = AX, \qquad (XA)^H = XA.\]

Penrose (1955) proved exactly one X satisfies all four.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (4,) – ||AXA - A||_F, ||XAX - X||_F, ||(AX)^H - AX||_F, ||(XA)^H - XA||_F.

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 2.0], [2.0, 4.0], [0.0, 1.0]])
>>> bool(np.all(penrose_residuals(A, pseudoinverse(A)) < 1e-12))
True
mathematicskit.linalg.power_iteration(a, tol=1e-12, max_iter=1000, v0=None)[source]#

Power iteration for the dominant (largest-magnitude) eigenpair.

\(v_{k+1} = A v_k / \|A v_k\|\) converges to the eigenvector of the eigenvalue of largest magnitude (assuming it is unique and the initial vector has a nonzero component along it), with the Rayleigh quotient \(v_k^T A v_k\) converging to the eigenvalue itself. No library equivalent for the iteration itself – watching this converge is the point. See Burden & Faires, Numerical Analysis, 10th ed., Ch. 9.3.

Parameters:

  • a (ndarray)

  • tol (float) – Convergence tolerance on the Rayleigh-quotient estimate.

  • max_iter (int)

  • v0 (ndarray | None) – Initial vector; a fixed pseudo-random default is used if omitted.

Return type:

EigenResult

Returns:

EigenResult – eigenvalues a length-1 array, eigenvectors shape (n,).

Examples

>>> import numpy as np
>>> A = np.array([[2.0, 0.0], [0.0, 5.0]])
>>> result = power_iteration(A)
>>> round(float(result.eigenvalues[0]), 6)
5.0
mathematicskit.linalg.pseudoinverse(a, rcond=None)[source]#

Moore-Penrose pseudoinverse \(A^+ = V \Sigma^+ U^T\), via numpy.linalg.pinv().

For any right-hand side b, \(x = A^+ b\) is the minimum-norm least-squares solution of A x = b – it exists and is unique even when A is rectangular or rank-deficient.

Parameters:

  • a (ndarray)

  • rcond (float | None) – Singular values below rcond * max(S) are treated as zero; defaults to numpy’s max(m, n) * eps.

Return type:

ndarray

Returns:

ndarray, shape (n, m)

Examples

>>> import numpy as np
>>> A = np.array([[1.0, 1.0], [1.0, 1.0]])  # rank 1
>>> pseudoinverse(A)
array([[0.25, 0.25],
       [0.25, 0.25]])
mathematicskit.linalg.random_spd_matrix(n, seed=0, condition_scale=1.0)[source]#

Generate a random symmetric positive-definite matrix.

Constructs \(A = M^T M + \epsilon I\) for a random \(M\), which is SPD for any nonsingular \(M\) (\(x^T A x = \|Mx\|^2 + \epsilon\|x\|^2 > 0\) for \(x \neq 0\)); condition_scale scales \(M\)’s entries to make the resulting matrix better- or worse-conditioned for stability demonstrations.

Parameters:

  • n (int) – Matrix size.

  • seed (int) – Random seed, for reproducibility.

  • condition_scale (float) – Multiplies the generating matrix’s entries.

Return type:

ndarray

Returns:

ndarray, shape (n, n)

Examples

>>> import numpy as np
>>> A = random_spd_matrix(4, seed=1)
>>> np.allclose(A, A.T)
True
>>> bool(np.all(np.linalg.eigvalsh(A) > 0))
True
mathematicskit.linalg.schur_decompose(a, output='real')[source]#

Schur decomposition A = Z T Z^H via scipy.linalg.schur().

Parameters:

  • a (ndarray)

  • output (str) – "real" keeps a real T with 2x2 diagonal blocks for complex-conjugate eigenvalue pairs; "complex" returns a truly upper-triangular T.

Return type:

SchurResult

Returns:

SchurResult

Examples

>>> import numpy as np
>>> A = np.array([[4.0, 1.0], [2.0, 3.0]])  # eigenvalues 2 and 5
>>> result = schur_decompose(A)
>>> np.round(np.sort(np.diag(result.T)), 8)
array([2., 5.])
>>> bool(np.allclose(result.Z @ result.T @ result.Z.T, A))
True
mathematicskit.linalg.svd_decompose(a)[source]#

Compute the (thin) SVD A = U diag(S) V^T via numpy.linalg.svd().

full_matrices=False requests the thin/economy factorization (U shape (m, k), Vt shape (k, n) with k = min(m, n)), which is what every mathematicskit.linalg caller needs. Singular values are returned in descending order (numpy’s convention).

Parameters:

a (ndarray)

Return type:

SVDResult

Returns:

SVDResult

Examples

>>> import numpy as np
>>> A = np.array([[3.0, 0.0], [0.0, -2.0]])
>>> result = svd_decompose(A)
>>> np.round(result.S, 6)
array([3., 2.])
>>> np.allclose(result.U @ np.diag(result.S) @ result.Vt, A)
True
On this page
Show Source