Advanced Linear Algebra

SVD decomposes any m×n matrix A = UΣV^T, where U is m×m orthogonal (left singular vectors), Σ is m×n diagonal with non-negative singular values σ₁≥σ₂≥...≥σ_min(m,n)≥0, and V is n×n orthogonal (right singular vectors). SVD reveals: rank (number of non-zero singular values), spectral norm (σ₁), Frobenius norm (√Σσᵢ²), condition number (σ₁/σ_min), and provides the best low-rank approximation by truncating to the top k singular values. Every matrix has an SVD, unlike eigendecomposition (requires square matrix). Geometrically, SVD shows A maps the unit sphere to an ellipsoid with semi-axes σᵢ along directions uᵢ.

import numpy as np
A = np.random.randn(5, 3)
U, s, Vh = np.linalg.svd(A, full_matrices=False)
print(U.shape, s.shape, Vh.shape)  # (5,3) (3,) (3,3)
print(np.allclose(A, U @ np.diag(s) @ Vh))  # True

The Eckart-Young-Mirsky theorem states that the best rank-k approximation of matrix A in both the Frobenius and spectral norm is A_k = U_k Σ_k V_k^T, where U_k, Σ_k, V_k are the truncated SVD keeping only the top k singular values/vectors. Formally: ||A - B||_F ≥ ||A - A_k||_F for any rank-k matrix B. The approximation error in Frobenius norm is √(σ_{k+1}² + ... + σ_r²). This theorem justifies PCA (keep top k principal components), image compression (store only k singular triplets instead of m×n pixels), and low-rank matrix factorization in recommendation systems—all retain the most information per rank unit.

import numpy as np
def low_rank_approx(A, k):
    U, s, Vh = np.linalg.svd(A, full_matrices=False)
    return U[:,:k] @ np.diag(s[:k]) @ Vh[:k,:]
A = np.random.randn(100,50)
A5 = low_rank_approx(A, 5)
print(np.linalg.matrix_rank(A5))  # 5

Eigendecomposition of a square matrix A produces A = PΛP^{-1}, where Λ is diagonal (eigenvalues) and P's columns are eigenvectors. A has a complete eigendecomposition (is diagonalizable) if and only if it has n linearly independent eigenvectors—equivalently, every eigenvalue's geometric multiplicity equals its algebraic multiplicity. Sufficient conditions: A is symmetric (eigenvectors orthogonal, real eigenvalues by spectral theorem); A has n distinct eigenvalues. Matrices that are not diagonalizable require the Jordan normal form. In ML: covariance matrices are symmetric positive semi-definite, always diagonalizable; transition matrices in Markov chains must be analyzed with eigendecomposition.

import numpy as np
A = np.array([[4,1],[2,3]], dtype=float)
vals, vecs = np.linalg.eig(A)
print(vals)  # eigenvalues
print(np.allclose(A @ vecs, vecs @ np.diag(vals)))  # True

A symmetric matrix A is positive definite (PD) if x^TAx > 0 for all non-zero vectors x. Equivalent characterizations: (1) All eigenvalues are strictly positive. (2) All leading principal minors (determinants of top-left k×k submatrices) are positive (Sylvester's criterion). (3) A has a Cholesky decomposition A = LL^T with positive diagonal entries in L. A is positive semi-definite (PSD) if x^TAx ≥ 0 (eigenvalues ≥ 0, Cholesky with non-negative diagonal). PD matrices arise as covariance matrices, Gram matrices of linearly independent vectors, and Hessians at strict local minima. Covariance matrices are always PSD, strictly PD when features are linearly independent.

import numpy as np
def is_positive_definite(A):
    try:
        np.linalg.cholesky(A)
        return True
    except np.linalg.LinAlgError:
        return False
A = np.array([[4,2],[2,3.]])
print(is_positive_definite(A))  # True

QR decomposes an m×n matrix A (m≥n) as A=QR, where Q is m×m orthogonal and R is m×n upper triangular. Applications: solving least squares (QR is numerically stable), eigenvalue algorithms (QR iteration), and orthonormal basis computation. Two main algorithms: (1) Gram-Schmidt orthogonalization: process columns of A one by one, orthogonalizing against previous columns. Classical GS is unstable; Modified GS is numerically stable. O(mn²). (2) Householder reflections: apply n reflections to zero out below-diagonal entries. More numerically stable than GS, preferred in practice (used in LAPACK). O(2mn²-2n³/3). Both produce the same mathematical result but differ in stability.

import numpy as np
A = np.random.randn(10, 4)
Q, R = np.linalg.qr(A)
print(Q.shape, R.shape)  # (10,4) (4,4) reduced QR
print(np.allclose(A, Q@R))  # True
print(np.allclose(Q.T@Q, np.eye(4)))  # True: orthonormal cols

LU factorizes a square matrix as A = LU (L lower triangular with 1s on diagonal, U upper triangular). Without pivoting, if any diagonal pivot is zero or small, division by a small number causes numerical instability. Partial pivoting permutes rows so the largest magnitude element in each column becomes the pivot: PA = LU, where P is a permutation matrix. This bounds all entries of L to be ≤ 1 in magnitude, controlling error growth. LU with pivoting solves Ax=b in O(n³) and is numerically stable for most practical matrices. Full pivoting (permuting both rows and columns) is theoretically more stable but rarely necessary. Used in LAPACK's dgetrf and scipy.linalg.lu.

import numpy as np
from scipy.linalg import lu
A = np.array([[2.,1.],[1.,3.]])
P, L, U = lu(A)
print(np.allclose(P@A, L@U))  # True

Cholesky decomposes a symmetric positive definite (SPD) matrix A = LL^T, where L is lower triangular with positive diagonal. Since A is symmetric, only the lower triangular part needs to be stored (n(n+1)/2 entries vs n²). The algorithm runs in n³/3 FLOPs—roughly half of LU's 2n³/3 FLOPs—because symmetry is exploited. Cholesky is always numerically stable for SPD matrices (no pivoting needed) because positive definiteness bounds diagonal entries. Applications: solving SPD linear systems (Gaussian processes, Newton step with PD Hessian), multivariate Gaussian sampling (x=μ+Lz, z~N(0,I)), and log-determinant computation (log det A = 2 sum log(diag L)).

import numpy as np
A = np.array([[4.,2.],[2.,3.]])
L = np.linalg.cholesky(A)
print(L)  # lower triangular
print(np.allclose(L @ L.T, A))  # True
log_det = 2*np.sum(np.log(np.diag(L)))  # log det A

The Moore-Penrose pseudoinverse A⁺ extends the concept of matrix inverse to non-square and rank-deficient matrices. Via SVD A=UΣV^T: A⁺ = V Σ⁺ U^T, where Σ⁺ inverts non-zero diagonal entries of Σ and sets zero entries to zero. Properties: AA⁺A=A, A⁺AA⁺=A⁺, (AA⁺)^T=AA⁺, (A⁺A)^T=A⁺A. For full-column-rank A: A⁺=(A^TA)^{-1}A^T (left inverse, least-squares solution). For full-row-rank A: A⁺=A^T(AA^T)^{-1} (right inverse, minimum-norm solution). The least-squares solution to Ax=b with minimum norm is x=A⁺b. np.linalg.pinv uses SVD with a tolerance to threshold small singular values.

import numpy as np
A = np.array([[1,2],[3,4],[5,6.]])
Ap = np.linalg.pinv(A)
print(Ap.shape)  # (2,3)
# Least-norm least-squares solution:
b = np.array([1.,2.,3.])
x = Ap @ b
print(x)  # minimum norm solution to Ax=b

The spectral theorem states that every real symmetric matrix A is orthogonally diagonalizable: A = QΛQ^T, where Q is orthogonal (columns are eigenvectors) and Λ is diagonal (eigenvalues are all real). This means: (1) Symmetric matrices always have real eigenvalues. (2) Eigenvectors corresponding to distinct eigenvalues are orthogonal. (3) Eigenvectors form a complete orthonormal basis of R^n. The theorem extends to complex Hermitian matrices (A=A^H) with unitary U and real eigenvalues. In ML: the Hessian of a loss at any stationary point is symmetric, and its eigenvalues determine the curvature. Covariance matrices are symmetric PSD and diagonalized by their eigenvectors (principal components).

import numpy as np
A = np.array([[3,1,0],[1,4,1],[0,1,5.]])
vals, Q = np.linalg.eigh(A)  # eigh for symmetric
print(np.allclose(A, Q @ np.diag(vals) @ Q.T))  # True
print(np.allclose(Q.T @ Q, np.eye(3)))  # True: orthogonal

The trace of a square matrix tr(A) = Σ A_{ii} (sum of diagonal entries). Key properties: (1) tr(A) = Σλᵢ (sum of eigenvalues). (2) tr(AB) = tr(BA) (cyclic property; extends to tr(ABC) = tr(BCA) = tr(CAB)). (3) tr(A^T) = tr(A). Additional: tr(A+B) = tr(A)+tr(B) (linearity). The trace is used in: computing expected squared norms (E[||Ax||²] = tr(A^TA) for unit-variance x), the Hutchinson estimator for stochastic trace computation, Gaussian process marginal likelihood (tr(K^{-1}y y^T)), and as a matrix norm proxy. The Frobenius norm ||A||_F = √tr(A^TA).

import numpy as np
A = np.random.randn(4,4)
B = np.random.randn(4,4)
print(np.allclose(np.trace(A@B), np.trace(B@A)))  # True
vals = np.linalg.eigvals(A)
print(np.allclose(np.trace(A), np.sum(vals.real)))  # True for real A

The determinant det(A) of an n×n matrix measures the signed volume scaling factor of the linear transformation A: the image of the unit hypercube has volume |det(A)|. det(A) = Π λᵢ (product of eigenvalues). det(AB) = det(A)·det(B). det(A^T) = det(A). If det(A)=0, A is singular (not invertible), and one eigenvalue is zero. log det(A) is used in multivariate Gaussian log-likelihoods: log p(x|μ,Σ) = -½(log det(Σ) + (x-μ)^TΣ^{-1}(x-μ) + n log 2π). Efficient computation via Cholesky: log det(A) = 2·Σlog(diag(L)) for SPD A=LL^T, avoiding O(n³) direct determinant computation.

import numpy as np
A = np.array([[4.,2.],[2.,3.]])
L = np.linalg.cholesky(A)
log_det = 2*np.sum(np.log(np.diag(L)))
print(log_det, np.log(np.linalg.det(A)))  # same value

The Kronecker product A⊗B (where A is m×n and B is p×q) produces an mp×nq block matrix where each aᵢⱼ·B forms a p×q block. Key properties: (A⊗B)(C⊗D) = (AC)⊗(BD) (mixed product); (A⊗B)^T = A^T⊗B^T; (A⊗B)^{-1} = A^{-1}⊗B^{-1} (when invertible). In KFAC (Kronecker-Factored Approximate Curvature), the Fisher Information Matrix block for a layer is approximated as A⊗G (activation statistics ⊗ gradient statistics), enabling efficient O(n^{1.5}) inversion instead of O(n³). Also used in tensorized neural networks (replacing weight matrices with Kronecker products of smaller matrices) to reduce parameters.

import numpy as np
A = np.array([[1,2],[3,4]])
B = np.array([[0,5],[6,7]])
K = np.kron(A, B)
print(K.shape)  # (4,4)
# Property: (A⊗B)(C⊗D) = (AC)⊗(BD)
C = np.eye(2); D = np.eye(2)
print(np.allclose(np.kron(A@C, B@D), np.kron(A,B)@np.kron(C,D)))  # True

The rank of a matrix A (m×n) is the dimension of its column space (= row space dimension). The null space is {x: Ax=0}; its dimension is the nullity. Rank-nullity theorem: rank(A) + nullity(A) = n (number of columns). Implications: if rank(A)=n, A has trivial null space (full column rank), and Ax=b has at most one solution. If rank(A)=m, A has full row rank, and Ax=b always has a solution. Rank can be computed via SVD (count non-zero singular values, with tolerance ε·max singular value for numerical rank). A matrix is rank-deficient if rank < min(m,n); columns (or rows) are linearly dependent.

import numpy as np
A = np.array([[1,2,3],[4,5,6],[7,8,9.]])
print(np.linalg.matrix_rank(A))  # 2 (rank-deficient)
_, s, _ = np.linalg.svd(A)
print(s)  # [1.68e+01, 1.07e+00, 4.21e-16] -> 2 non-zero

Matrix calculus extends scalar calculus to matrix-valued functions. For scalar f(X) where X is m×n: ∂f/∂X is an m×n matrix with (i,j) entry ∂f/∂X_{ij}. Key identities: ∂tr(AX)/∂X = A^T. ∂tr(X^TAX)/∂X = (A+A^T)X = 2AX (if A symmetric). ∂log det(X)/∂X = X^{-T} = (X^{-1})^T. ∂||X||_F²/∂X = 2X. These are essential for deriving backpropagation equations, computing gradients of regularization terms, and deriving closed-form solutions (e.g., linear regression: setting ∂||Xw-y||²/∂w = 0 gives w=(X^TX)^{-1}X^Ty).

import numpy as np
# Verify: d/dX tr(AX) = A^T
A = np.random.randn(3,3)
X = np.random.randn(3,3)
eps = 1e-7
numerical_grad = np.zeros_like(X)
for i in range(3):
    for j in range(3):
        X[i,j]+=eps; fp=np.trace(A@X)
        X[i,j]-=2*eps; fm=np.trace(A@X)
        X[i,j]+=eps
        numerical_grad[i,j]=(fp-fm)/(2*eps)
print(np.allclose(numerical_grad, A.T, atol=1e-5))  # True

The nuclear norm ||A||_* = Σσᵢ (sum of all singular values) is the convex envelope of matrix rank—it is the tightest convex lower bound on the rank function over matrices with spectral norm ≤ 1. Minimizing the nuclear norm promotes low-rank solutions, serving as a convex relaxation of rank minimization (analogous to L1 norm as a convex relaxation of L0 for vectors). Applications: matrix completion (Netflix prize: fill in missing ratings by minimizing nuclear norm), robust PCA (decompose matrix as low-rank + sparse), and collaborative filtering. The subdifferential of ||A||_* at A = UΣV^T is UV^T + W where ||W||≤1 and U^TW=WV=0.

import numpy as np
A = np.random.randn(5, 4)
_, s, _ = np.linalg.svd(A)
nuclear_norm = np.sum(s)
print(nuclear_norm)

The Gram matrix G of a set of vectors {x₁,...,xₙ} is G_{ij} = xᵢ·xⱼ, i.e., G = X^TX (or XX^T depending on orientation). Properties: (1) G is always symmetric: G^T = G. (2) G is always positive semi-definite: for any vector v, v^TGv = ||Xv||² ≥ 0. (3) rank(G) = rank(X) (same rank). (4) If X has linearly independent rows, G is positive definite. Gram matrices are central to kernel methods: the kernel matrix K_{ij} = k(xᵢ,xⱼ) must be a valid Gram matrix (PSD) for Mercer's theorem to apply. Gram matrices enable computing inner products in feature space without explicitly computing the features (kernel trick).

import numpy as np
X = np.random.randn(50, 10)  # 50 samples, 10 features
G = X @ X.T  # Gram matrix (50×50)
eigenvalues = np.linalg.eigvalsh(G)
print(np.all(eigenvalues >= -1e-10))  # True: PSD
print(G.shape)  # (50, 50)

The Woodbury identity states: (A + UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}, where A is n×n, U is n×k, C is k×k, V is k×n, with k << n. This converts an n×n matrix inversion into a k×k inversion, saving O(n³-k³) FLOPs. Special case (Sherman-Morrison): (A + uv^T)^{-1} = A^{-1} - A^{-1}uv^TA^{-1}/(1+v^TA^{-1}u) for rank-1 update. ML applications: Gaussian process posterior updates (adding a new data point), ridge regression with kernel trick (n×n vs d×d inversion), online Bayesian updating, and low-rank corrections to large systems.

import numpy as np
n,k=100,5
A=np.eye(n); U=np.random.randn(n,k); V=U.T
# (I + U V)^-1 via Woodbury:
C=np.eye(k)
Ainv=A  # A=I so A^-1=I
result = Ainv - Ainv@U@np.linalg.inv(C+V@Ainv@U)@V@Ainv
print(np.allclose(np.linalg.inv(A+U@V), result))  # True

The spectral radius ρ(A) = max|λᵢ| (maximum absolute eigenvalue). For a square matrix A, the sequence ||A^k||^{1/k} → ρ(A) as k → ∞ (Gelfand formula). Power iteration estimates the dominant eigenvector: repeatedly multiply a random vector by A and normalize. It converges at rate (|λ₂|/|λ₁|)^k—fast when the spectral gap (|λ₁|-|λ₂|) is large. Shift-and-invert accelerates this for other eigenvalues. The spectral radius determines whether a matrix iteration x_{t+1} = Ax_t converges (ρ<1), diverges (ρ>1), or oscillates (ρ=1). In recurrent networks, the spectral radius of the weight matrix is crucial for stability.

import numpy as np
def power_iteration(A, n_iter=1000):
    x = np.random.randn(A.shape[0])
    for _ in range(n_iter):
        x = A @ x
        x /= np.linalg.norm(x)
    eigenval = x @ A @ x
    return eigenval, x
A = np.random.randn(5,5); A = A.T@A
val, vec = power_iteration(A)
print(val, np.linalg.eigvalsh(A)[-1])  # should match

The matrix exponential exp(A) = Σ A^k/k! (k from 0 to ∞). For a diagonalizable A = PΛP^{-1}: exp(A) = P·diag(exp(λᵢ))·P^{-1}. The system of ODEs ẋ = Ax with x(0) = x₀ has exact solution x(t) = exp(At)·x₀. This is used in continuous-time dynamical systems, Markov chain analysis (transition rates), and neural ODEs (where the dynamics are parameterized by a neural network). Computing exp(At) efficiently for large sparse A uses Krylov subspace methods (Expokit). For stability: x(t) remains bounded iff all eigenvalues of A have non-positive real parts—the spectral abscissa determines stability.

import numpy as np
from scipy.linalg import expm
A = np.array([[-1,2],[0,-3.]])
t = 0.5
print(expm(A*t))  # solution operator for x' = Ax

The Schur decomposition states every square matrix A = QTQ^H, where Q is unitary and T is upper triangular (Schur form). Diagonal entries of T are the eigenvalues of A. Unlike eigendecomposition, Schur decomposition always exists (even for defective matrices) and Q is always unitary (orthogonal for real normal matrices). For normal matrices (AA^H = A^HA, including symmetric, unitary, skew-symmetric), the Schur form is diagonal—equivalent to eigendecomposition. The real Schur decomposition of a real matrix uses orthogonal Q and quasi-upper-triangular T (with 2×2 blocks for complex conjugate eigenvalue pairs). LAPACK's schur routines are the basis for most numerical eigenvalue computations.

from scipy.linalg import schur
import numpy as np
A = np.array([[2,1,0],[0,3,1],[0,0,4.]])
T, Z = schur(A)
print(T)  # upper triangular
print(np.allclose(Z@T@Z.T, A))  # True

Given centered data X (n×d), PCA finds directions of maximum variance. Covariance matrix C = X^TX/(n-1). Eigendecomposition C = QΛQ^T gives principal components (columns of Q, ordered by eigenvalue). SVD of X: X = UΣV^T. The principal components are V (right singular vectors of X). The eigenvalues of C equal Σ²/(n-1). Using SVD directly on X avoids forming X^TX (which loses precision for nearly-collinear features) and is more numerically stable. The proportion of variance explained by k components = Σᵢ₌₁ᵏ σᵢ²/Σ σᵢ². Sklearn's PCA uses truncated SVD (randomized SVD for large datasets).

import numpy as np
X = np.random.randn(100, 10)
X -= X.mean(0)  # center
U, s, Vh = np.linalg.svd(X, full_matrices=False)
# Principal components are rows of Vh (or cols of Vh.T)
# Projected data (scores):
scores = X @ Vh.T  # (100, 10)
var_explained = s**2 / (len(X)-1)
print(var_explained / var_explained.sum())

For an m×n matrix A with rank r, the four fundamental subspaces (Gilbert Strang) are: (1) Column space C(A): span of columns, dimension r. (2) Null space N(A): {x: Ax=0}, dimension n-r (by rank-nullity). (3) Row space C(A^T): span of rows = column space of A^T, dimension r. (4) Left null space N(A^T): {y: A^Ty=0}, dimension m-r. Orthogonality: N(A) ⊥ C(A^T) in R^n, and N(A^T) ⊥ C(A) in R^m. Together they form orthogonal complements: R^n = C(A^T) ⊕ N(A) and R^m = C(A) ⊕ N(A^T). The system Ax=b is consistent iff b ∈ C(A), i.e., b ⊥ N(A^T).

import numpy as np
A = np.array([[1,2,3],[4,5,6],[7,8,9.]])
# Null space via SVD (right singular vectors for zero singular values)
_, s, Vh = np.linalg.svd(A)
null_space = Vh[np.abs(s)<1e-10]  # or last rows if rank < n
print('rank:', np.sum(s>1e-10))  # 2
print('nullity:', 3-2)           # 1

The Frobenius norm ||A||_F = √(Σᵢⱼ aᵢⱼ²) = √(tr(A^TA)) = √(Σσᵢ²). Properties: it is unitarily invariant (||UAV||_F = ||A||_F for unitary U, V); it is the Hilbert-Schmidt norm for linear operators; it satisfies ||AB||_F ≤ ||A||_F·||B||_F (submultiplicativity). Relation to other norms: ||A||₂ (spectral) ≤ ||A||_F ≤ √min(m,n)·||A||₂. In optimization: regularizing with ||W||_F² = tr(W^TW) is L2 weight decay; gradient of ||W||_F² w.r.t. W = 2W. The Frobenius norm is the natural metric for comparing matrices as vectors in R^{mn} and measures total signal energy across all modes.

import numpy as np
A = np.random.randn(4, 5)
frob = np.linalg.norm(A, 'fro')
_, s, _ = np.linalg.svd(A)
print(frob, np.sqrt(np.sum(s**2)))  # same
print(frob, np.sqrt(np.trace(A.T@A)))  # same

The orthogonal projector onto the column space of A (full column rank) is P = A(A^TA)^{-1}A^T. Via QR (A=QR): P = QQ^T. Two key properties: (1) Idempotency: P² = P (projecting twice gives the same result). (2) Symmetry: P^T = P. The complementary projector I-P projects onto the orthogonal complement (null space of A^T). These are the only conditions: any symmetric idempotent matrix is an orthogonal projector onto its column space. Eigenvalues of P are 0 and 1 (ranks determine multiplicities). In linear regression, Ŷ = Py is the projection of y onto the column space of X (hat matrix P = X(X^TX)^{-1}X^T).

import numpy as np
A = np.random.randn(10, 3)
Q, _ = np.linalg.qr(A)
P = Q @ Q.T  # projector onto col(A)
print(np.allclose(P@P, P))   # True: idempotent
print(np.allclose(P, P.T))   # True: symmetric

The Kronecker product A⊗B is the concrete matrix representation of the tensor product when A and B are finite-dimensional matrices. The tensor product is an abstract algebraic construction: V⊗W is a vector space with basis {eᵢ⊗eⱼ} where {eᵢ} and {eⱼ} are bases of V and W. A linear map f:V→V' and g:W→W' extend to f⊗g:V⊗W→V'⊗W', represented as the Kronecker product in coordinates. In deep learning, higher-order tensors (e.g., weight tensors of dimension d₁×d₂×d₃) are elements of V₁⊗V₂⊗V₃. Tensor decompositions (Tucker, CP/PARAFAC) generalize matrix SVD to higher-order tensors for parameter compression.

The Perron-Frobenius theorem states: a square matrix with all positive entries has a unique largest real eigenvalue ρ (Perron root, equal to the spectral radius). Its corresponding eigenvectors have all positive entries. For non-negative irreducible matrices (all entries ≥ 0 and the matrix cannot be block-triangularized by permutation): ρ is a simple eigenvalue, and the corresponding eigenvector is entry-wise positive. Applications: Google PageRank algorithm (stochastic matrix; stationary distribution is the principal eigenvector), Markov chain stationary distributions, topic models, and nonnegative matrix factorization. The power iteration method converges to the Perron eigenvector for positive matrices.

Three important matrix norms in ML: (1) Spectral norm ||A||₂ = σ_max (largest singular value) = the induced ℓ² operator norm, measuring maximum stretching of unit vectors. (2) Frobenius norm ||A||_F = √(Σσᵢ²) = Hilbert-Schmidt norm, measuring total energy. (3) Nuclear norm ||A||_* = Σσᵢ, sum of all singular values, a convex proxy for rank. Inequalities: ||A||₂ ≤ ||A||_F ≤ √r·||A||₂ and ||A||_F ≤ ||A||_* ≤ √r·||A||_F, where r = rank. Unitarily invariant: all three are unchanged by left/right multiplication by unitary matrices. They are dual norms: ||A||_* = max_{||B||₂≤1} tr(A^TB) and ||A||_F = max_{||B||_F≤1} tr(A^TB).

import numpy as np
A = np.random.randn(5, 4)
_, s, _ = np.linalg.svd(A)
spectral = s[0]      # sigma_max
frobenius = np.sqrt(np.sum(s**2))
nuclear = np.sum(s)
print(spectral, frobenius, nuclear)
print(spectral <= frobenius <= nuclear)  # True

The condition number κ(A) = ||A||·||A^{-1}|| (for the spectral norm: σ_max/σ_min) measures how much relative error in the right-hand side b can amplify into relative error in the solution x: ||δx||/||x|| ≤ κ(A)·||δb||/||b||. A condition number of 10^k means up to k digits of precision can be lost. κ=1 for orthogonal matrices (perfect conditioning). For floating-point arithmetic (machine epsilon ε_mach ≈ 10^{-16}), solutions are unreliable when κ > 1/ε_mach ≈ 10^{16}. Preconditioning (replacing Ax=b with PAx=Pb for preconditioner P≈A^{-1}) reduces κ. Always check κ before trusting a linear solve.

import numpy as np
A = np.array([[1,0],[0,1e-10]])
print(np.linalg.cond(A))  # 1e10 - very ill-conditioned
b = np.array([1., 1e-10])
x = np.linalg.solve(A, b)
print(x)  # [1., 1.] but accuracy limited by cond

The vec operator stacks the columns of a matrix A (m×n) into a column vector of length mn: vec(A) = [a₁; a₂; ...; aₙ] where aᵢ is the i-th column. Key identity: vec(AXB) = (B^T⊗A)vec(X), converting a matrix equation AXB=C into the linear system (B^T⊗A)vec(X)=vec(C). This enables solving Sylvester equations AX+XB=C and Lyapunov equations AX+XA^T=Q by vectorization followed by Gaussian elimination on the mn×mn system. Also used in matrix Gaussian distributions where the covariance of vec(X) is a Kronecker product. In deep learning, gradient computations with respect to weight matrices are naturally expressed via vec and Kronecker products.

import numpy as np
A = np.array([[1,2],[3,4.]])
X = np.array([[1,0],[0,1.]])
B = np.array([[2,0],[1,3.]])
lhs = np.kron(B.T, A) @ A.flatten(order='F')
rhs = (A @ X @ B).flatten(order='F')
print(np.allclose(lhs, rhs))  # True

Cayley-Hamilton states that every square matrix satisfies its own characteristic polynomial: if p(λ) = det(A-λI) = (-1)^n(λ^n - tr(A)λ^{n-1} + ... + det(A)) is the characteristic polynomial, then p(A) = 0. Practical applications: (1) Computing A^{-1} from the characteristic polynomial: if det(A)≠0, p(A)=0 gives A^{-1} expressed as a polynomial in A of degree n-1. (2) Reducing high powers of A: A^n = linear combination of A^{n-1},...,A,I (using p(A)=0 to reduce degree). (3) Matrix functions f(A) via polynomial interpolation on eigenvalues. In theoretical computer science, Cayley-Hamilton underpins algorithmic speedups for matrix power computations.

import numpy as np
A = np.array([[1,2],[3,4.]])
# Char poly for 2x2: lambda^2 - tr(A)*lambda + det(A) = 0
# Verify: A^2 - tr(A)*A + det(A)*I = 0
result = A@A - np.trace(A)*A + np.linalg.det(A)*np.eye(2)
print(np.allclose(result, 0))  # True

Randomized SVD (Halko et al. 2011) computes a rank-k approximation of an m×n matrix A in O(mn log k + (m+n)k²) time, faster than full SVD's O(mn·min(m,n)) for k << min(m,n). Algorithm: (1) Draw random matrix Ω (n×(k+p), p=oversampling). (2) Form Y=AΩ (m×(k+p)). (3) Compute QR of Y to get Q (m×(k+p)). (4) Form B=Q^TA (small (k+p)×n matrix). (5) SVD of B gives approximate top-k singular values/vectors. The randomized range finder captures the column space in just 2 passes over A. Used in sklearn PCA (n_components < n_features) and large-scale recommendation systems where the full SVD is impractical.

import numpy as np
from sklearn.utils.extmath import randomized_svd
A = np.random.randn(1000, 500)
U, s, Vt = randomized_svd(A, n_components=10, random_state=42)
print(U.shape, s.shape, Vt.shape)  # (1000,10) (10,) (10,500)

A quadratic form Q(x) = x^TAx (A symmetric) can be diagonalized via eigendecomposition A=QΛQ^T: substituting x=Qy gives Q(x) = y^TΛy = Σλᵢyᵢ². In principal axes (eigenvector coordinates), the quadratic form has no cross terms. The level set {x: x^TAx = 1} for positive definite A is an ellipsoid with semi-axes 1/√λᵢ along eigenvectors qᵢ. Longest axis corresponds to smallest eigenvalue. This characterizes: covariance ellipses in Gaussian distributions, optimization landscape curvature (quadratic loss level sets), and Mahalanobis distance x^TΣ^{-1}x (which defines the metric for Gaussian distributions with covariance Σ).

import numpy as np
Sigma = np.array([[4,2],[2,3.]])
vals, vecs = np.linalg.eigh(Sigma)
print('Semi-axes:', 1/np.sqrt(vals))   # ellipse semi-axes
# Mahalanobis distance of x from origin:
x = np.array([1., 1.])
mahal = x @ np.linalg.inv(Sigma) @ x
print(mahal)

The QR algorithm iteratively computes QR decompositions of A: A_k = Q_kR_k; A_{k+1} = R_kQ_k = Q_k^TA_kQ_k. Each step is an orthogonal similarity transformation, preserving eigenvalues. For generic matrices, A_k converges to quasi-upper-triangular (Schur) form where eigenvalues appear on the diagonal. Practical QR algorithm: (1) Reduce to Hessenberg form (upper triangular + one sub-diagonal) in O(n³) to make each QR step O(n²). (2) Apply shifts (Francis implicit double-shift QZ) for quadratic convergence. Total cost O(n³). The QR algorithm is the standard eigenvalue solver in LAPACK (dgeev, dsyev) because it handles all cases (symmetric/non-symmetric, real/complex), is numerically backward stable, and converges cubically with good shifts.

The Lanczos algorithm builds a Krylov subspace K_m(A, v) = span{v, Av, A²v,...,A^{m-1}v} and finds the projection of A onto this subspace as a symmetric tridiagonal matrix T_m. Eigenvalues of T_m (Ritz values) approximate extreme eigenvalues of A with high accuracy after m << n iterations. Unlike power iteration (one eigenvalue at a time), Lanczos simultaneously captures the k largest and smallest eigenvalues. Cost per iteration: one matrix-vector multiply O(nnz(A)) + O(mn) orthogonalization. In practice, re-orthogonalization prevents loss of numerical orthogonality. Used in: principal eigenvalue computation for large sparse matrices, conjugate gradient (Lanczos is the theoretical foundation), and spectral graph methods.

Direct methods (Gaussian elimination, LU, Cholesky) produce an exact solution (up to floating-point errors) in a fixed number of operations: O(n³) for dense systems. They are preferred when n < ~10,000 and A is dense or when multiple right-hand sides need solving with the same A (factorize once, solve many times via triangular substitution in O(n²) each). Iterative methods (Conjugate Gradient, GMRES, MINRES) generate a sequence of approximations, converging in at most n steps but often much fewer for well-conditioned or structured (sparse/banded) A. They are preferred for very large sparse systems (n > 10^5) where direct methods cannot store the factor. Preconditioning (PA=I_approx) is critical for iterative method convergence speed.

The best rank-1 approximation of A in Frobenius norm is σ₁u₁v₁^T (outer product of first left and right singular vectors, scaled by first singular value). This follows from Eckart-Young: among all rank-1 matrices B, ||A-B||_F ≥ ||A - σ₁u₁v₁^T||_F = √(Σᵢ₌₂σᵢ²). Equivalently, σ₁ = max_{||u||=||v||=1} u^TAv (variational characterization): the singular value problem maximizes the Rayleigh quotient of the bilinear form. Power method for SVD alternates: v_{k+1} = A^T(Av_k)/||A^T(Av_k)||, converging at rate (σ₂/σ₁)²—fast when there is a large gap between top two singular values.

import numpy as np
A = np.random.randn(10, 8)
U, s, Vh = np.linalg.svd(A, full_matrices=False)
A_rank1 = s[0] * np.outer(U[:,0], Vh[0,:])
print(np.linalg.matrix_rank(A_rank1))  # 1
error = np.linalg.norm(A - A_rank1, 'fro')
print(np.allclose(error, np.sqrt(np.sum(s[1:]**2))))  # True

A Toeplitz matrix has constant diagonals: T_{ij} = t_{i-j}. It arises in convolution (a 1D convolution Toeplitz matrix multiplied by a vector = convolution), time series AR models, and signal processing. An n×n Toeplitz matrix is fully specified by its 2n-1 entries. Key fast algorithms: (1) Toeplitz matrix-vector product via FFT in O(n log n) (embed in a 2n×2n circulant matrix, use FFT). (2) Solving Tx=b using Levinson-Durbin in O(n²) vs O(n³) for dense LU. In deep learning, Toeplitz-structured weight matrices enable O(n log n) matrix multiplication, used in efficient convolutions and compression of fully-connected layers.

import numpy as np
from scipy.linalg import toeplitz
c = np.array([4,3,2,1.])
r = np.array([4,5,6,7.])
T = toeplitz(c, r)
print(T)  # constant diagonals
# Matrix-vector product via FFT (for large n):
x = np.random.randn(4)
y = T @ x  # naive O(n^2)

The matrix square root A^{1/2} satisfies (A^{1/2})² = A. For symmetric positive semi-definite A = QΛQ^T, the unique symmetric PSD square root is A^{1/2} = Q Λ^{1/2} Q^T (take elementwise square roots of eigenvalues). For general square matrices, square roots may not be unique or may not exist (e.g., a nilpotent matrix with no square root over reals). The Newton-Schulz iteration computes A^{1/2} iteratively. Applications: whitening transformation (covariance^{-1/2} decorrelates features), sampling from multivariate Gaussians (use Cholesky or symmetric square root), and computing matrix functions in operator theory and quantum information.

import numpy as np
A = np.array([[4,2],[2,3.]])
vals, Q = np.linalg.eigh(A)
A_sqrt = Q @ np.diag(np.sqrt(vals)) @ Q.T
print(np.allclose(A_sqrt @ A_sqrt, A))  # True
# Whitening:
X = np.random.randn(100,2)
W = np.linalg.inv(A_sqrt)
X_white = X @ W.T  # decorrelated

The generalized eigenvalue problem Av = λBv arises when B ≠ I. It appears in Fisher's Linear Discriminant Analysis (maximize between-class scatter relative to within-class scatter: S_B v = λ S_W v) and in quantum mechanics. Solved by: (1) Cholesky factorize B = LL^T (requires B PD); (2) Transform A' = L^{-1}A L^{-T}; (3) Solve standard eigenvalue problem A'u = λu; (4) Recover v = L^{-T}u. Or use scipy.linalg.eig(A, B) directly. For symmetric PD B, scipy.linalg.eigh(A, B) is preferred (guaranteed real eigenvalues, numerically stable). LDA in sklearn uses the generalized eigenvalue formulation internally.

import numpy as np
from scipy.linalg import eigh
A = np.random.randn(5,5); A = A.T@A + np.eye(5)
B = np.random.randn(5,5); B = B.T@B + 2*np.eye(5)
vals, vecs = eigh(A, B)  # generalized eigenproblem
print(vals)  # real eigenvalues

The Gershgorin circle theorem states: every eigenvalue of a square matrix A lies within at least one Gershgorin disk D(aᵢᵢ, Rᵢ), where center aᵢᵢ is the i-th diagonal entry and radius Rᵢ = Σⱼ≠ᵢ |aᵢⱼ| (sum of off-diagonal absolute values in row i). If k disks form an isolated set from the remaining n-k disks, exactly k eigenvalues lie in their union. This gives a quick, cheap (no computation beyond summing entries) bound on eigenvalue locations. Applications: verifying positive definiteness (if all disks have positive centers and small radii so all disks are in the right half), stability analysis of dynamical systems, and preconditioning analysis.

import numpy as np
A = np.array([[10,1,0],[1,8,1],[0,1,9.]])
# Gershgorin disks
for i in range(3):
    center = A[i,i]
    radius = np.sum(np.abs(A[i])) - np.abs(A[i,i])
    print(f'Disk {i}: center={center}, radius={radius}')
print('Eigenvalues:', np.linalg.eigvalsh(A))

The Courant-Fischer (minimax) theorem characterizes eigenvalues of a symmetric matrix A as: λₖ = min_{dim(S)=k} max_{x∈S,||x||=1} x^TAx = max_{dim(S)=n-k+1} min_{x∈S,||x||=1} x^TAx. Equivalently: λₖ = max_{dim(S)=k} min_{x∈S,||x||=1} x^TAx. This variational characterization has applications: (1) Proving the smallest eigenvalue is the minimum of the Rayleigh quotient x^TAx/||x||². (2) Interlacing theorem: eigenvalues of a principal submatrix interlace with those of A. (3) Weyl's inequality: |λₖ(A+E) - λₖ(A)| ≤ ||E||₂, bounding eigenvalue perturbations. Used in spectral graph theory and quantum mechanics.

import numpy as np
A = np.array([[3.,1.],[1.,2.]])
vals, vecs = np.linalg.eigh(A)
# Verify: lambda_1 = min Rayleigh quotient
theta = np.linspace(0,2*np.pi,1000)
RQ = np.array([np.array([np.cos(t),np.sin(t)])@A@np.array([np.cos(t),np.sin(t)]) for t in theta])
print(vals[0], RQ.min())  # should match

The polar decomposition of an invertible matrix A = UP, where U is orthogonal (unitary) and P is symmetric positive definite. Via SVD A = QΣV^T: set U = QV^T (orthogonal) and P = VΣV^T (SPD, eigenvalues = singular values of A). P is the symmetric positive definite factor (stretching) and U represents pure rotation/reflection. Uniqueness: for invertible A, the polar decomposition is unique. Applications: (1) Orthogonal Procrustes problem: find orthogonal U minimizing ||A-UB||_F—solution is U = QV^T from SVD of AB^T. (2) Gradient of loss w.r.t. orthogonal constraints. (3) Computing the nearest orthogonal matrix to a given matrix.

import numpy as np
A = np.random.randn(4,4)
U_svd, s, Vh = np.linalg.svd(A)
U = U_svd @ Vh  # orthogonal factor
P = Vh.T @ np.diag(s) @ Vh  # SPD factor
print(np.allclose(A, U@P))  # True
print(np.allclose(U@U.T, np.eye(4)))  # True

The Nystrom method approximates an n×n PSD kernel matrix K using m << n landmark points. Select m landmarks X_m (randomly or via k-means); compute K_mn (m×n cross-kernel) and K_mm (m×m kernel). Approximate K ≈ K_mn^T K_mm^{-1} K_mn (rank-m approximation). This reduces storage from O(n²) to O(mn) and computation from O(n³) to O(m²n) for kernel methods. The approximation error decreases as more landmarks are added. Used in large-scale Gaussian processes (O(m²n) training vs O(n³)), kernel SVM on millions of points, and scalable kernel PCA. The quality of landmark selection critically affects approximation quality.

import numpy as np
def nystrom_approx(X, m, gamma=1.0, seed=42):
    n=len(X); rng=np.random.default_rng(seed)
    idx=rng.choice(n,m,replace=False)
    Xm=X[idx]
    def rbf(A,B): sq=((A[:,None]-B[None,:])**2).sum(-1); return np.exp(-gamma*sq)
    Kmm=rbf(Xm,Xm); Kmn=rbf(Xm,X)
    Kmm_inv=np.linalg.pinv(Kmm+1e-8*np.eye(m))
    return Kmn.T@Kmm_inv@Kmn  # (n,n) approx

The rank of matrix A equals the maximum number of linearly independent columns (= maximum number of linearly independent rows). Columns c₁,...,cₙ are linearly dependent if there exists a non-zero vector x such that Xc = 0, i.e., x is in the null space of A^T. Rank = n - dim(null space). Full column rank (rank = n): all columns independent; Ax=0 implies x=0; A^TA is invertible; unique least-squares solution. Full row rank (rank = m): all rows independent; A has a right inverse; Ax=b always has a solution. Determining linear independence computationally: compute SVD and count singular values above a tolerance threshold. Gram-Schmidt orthogonalization naturally discards linearly dependent vectors.

import numpy as np
def is_linearly_independent(vectors, tol=1e-10):
    A = np.column_stack(vectors)
    rank = np.linalg.matrix_rank(A, tol=tol)
    return rank == A.shape[1]  # True if all cols independent
v1,v2,v3=np.array([1,0,0.]),np.array([0,1,0.]),np.array([1,1,0.])
print(is_linearly_independent([v1,v2,v3]))  # False: v3=v1+v2

Hadamard's inequality states: |det(A)| ≤ Πᵢ ||aᵢ||₂, where aᵢ are the columns of A. Equality holds if and only if the columns are mutually orthogonal. This upper-bounds the volume of the parallelepiped spanned by the columns. Application: Hadamard matrices achieve the maximum determinant for {±1} matrices of size n, giving |det| = n^{n/2}. They are used in Hadamard transforms (O(n log n)), compressed sensing, error-correcting codes, and random feature maps. The Hadamard inequality also underpins the Gram determinant identity: det(A^TA) = (volume of parallelepiped)², connecting determinants to geometric volume.

import numpy as np
A = np.random.randn(4,4)
hadamard_bound = np.prod(np.linalg.norm(A,axis=0))
print(abs(np.linalg.det(A)), '<=', hadamard_bound)  # det <= bound

Matrix completion recovers a low-rank matrix M from a subset of observed entries (e.g., movie ratings). Direct rank minimization is NP-hard (combinatorial). Nuclear norm minimization provides a convex relaxation: minimize ||M||_* subject to M_{ij} = X_{ij} for observed entries (i,j). Under mild conditions (incoherence of the row/column spaces and sufficient random observations O(rn log n)), nuclear norm minimization recovers the exact rank-r matrix M. Solved efficiently with proximal gradient: the proximity operator of the nuclear norm is matrix soft-thresholding on singular values. Applications: recommendation systems (Netflix Prize), image inpainting, and collaborative filtering. Candès & Recht (2009) established the theoretical guarantees.

import numpy as np
def matrix_soft_threshold(A, tau):
    U, s, Vh = np.linalg.svd(A, full_matrices=False)
    s_thresh = np.maximum(s - tau, 0)
    return U @ np.diag(s_thresh) @ Vh  # prox of nuclear norm

For a regular Markov chain with transition matrix P (rows sum to 1, all entries ≥ 0), the stationary distribution π satisfies πP = π (left eigenvector for eigenvalue 1). The spectral gap δ = 1 - |λ₂| (where |λ₂| is the second-largest absolute eigenvalue) determines mixing time: the chain converges to stationarity in O(log(1/ε)/δ) steps. Large spectral gap (close to 1) = fast mixing; small spectral gap = slow mixing (poorly connected graph). The power method P^t π₀ → π converges at geometric rate |λ₂|^t. Spectral gap analysis is used in MCMC convergence analysis and random walks on graphs (PageRank, spectral clustering).

import numpy as np
P = np.array([[0.7,0.3],[0.4,0.6]])
vals, vecs = np.linalg.eig(P.T)  # left eigenvectors
idx = np.argmax(vals.real)
pi = vecs[:,idx].real; pi /= pi.sum()
print('Stationary dist:', pi)
print('Spectral gap:', 1-abs(sorted(vals.real)[-2]))

K-FAC (Martens & Grosse 2015) approximates the Fisher Information Matrix block for each layer as a Kronecker product: F_l ≈ A_{l-1} ⊗ G_l, where A_{l-1} = E[aₜ⁻¹aₜ⁻¹^T] (activation covariance) and G_l = E[δₜδₜ^T] (pre-activation gradient covariance). Inverting F_l = A⊗G uses the identity (A⊗G)^{-1} = A^{-1}⊗G^{-1}, reducing O(d_l^3) (full block inversion) to O(n_l^3 + n_{l-1}^3) (two small inversions). K-FAC gives natural gradient updates that are invariant to reparameterization and can converge in 10-100× fewer iterations than SGD on convex problems, at 2-3× cost per step versus SGD.

The matrix Cauchy-Schwarz inequality states: ||AB||_F ≤ ||A||_F·||B||_F (submultiplicativity of Frobenius norm). For vectors: |a·b| ≤ ||a||₂·||b||₂ (standard CS). For matrices: tr(A^TB)² ≤ tr(A^TA)·tr(B^TB) = ||A||_F²·||B||_F². Also: tr(A^TB) ≤ Σσᵢ(A)σᵢ(B) (von Neumann trace inequality). Applications in ML: (1) Bounding gradient norms: ||∇L||_F ≤ ||X||_F·||δ||_F (sensitivity to input). (2) Proving generalization bounds via Rademacher complexity. (3) Attention score bounds: query-key dot products are bounded by their respective norms. (4) Proving convergence of matrix-valued updates in optimization.

import numpy as np
A = np.random.randn(5, 4)
B = np.random.randn(4, 3)
# ||AB||_F <= ||A||_F * ||B||_F
lhs = np.linalg.norm(A@B, 'fro')
rhs = np.linalg.norm(A,'fro') * np.linalg.norm(B,'fro')
print(lhs, '<=', rhs, ':', lhs<=rhs)  # True

The covariance matrix Σ (d×d) has Σᵢⱼ = Cov(Xᵢ,Xⱼ) = E[(Xᵢ-μᵢ)(Xⱼ-μⱼ)], with variances on the diagonal. It is scale-dependent. The correlation matrix R has Rᵢⱼ = Σᵢⱼ/(σᵢσⱼ), normalizing to unit variance, making it scale-invariant and comparable across features. R = D^{-1/2} Σ D^{-1/2} where D = diag(Σ). Both are symmetric positive semi-definite. In ML: PCA uses the covariance matrix (or its eigendecomposition via SVD); Gaussian graphical models use the precision matrix Σ^{-1}; feature selection uses the correlation matrix to identify redundant features. Whitening (normalizing by Σ^{-1/2}) produces identity covariance and decorrelates features.

import numpy as np
X = np.random.randn(100, 5)
X -= X.mean(0)  # center
Cov = X.T @ X / (len(X)-1)
stds = np.sqrt(np.diag(Cov))
Corr = Cov / np.outer(stds, stds)
print(np.diag(Corr))  # all ones