Matrix Operations

Chapter 1: LU Decomposition

We know the goal: factor A = LU so that solving Ax = b reduces to two cheap triangular solves. But what are L and U, and how do we find them?

Lower triangular means all entries above the diagonal are zero. Upper triangular means all entries below the diagonal are zero. The crucial property: solving a triangular system takes only O(n²) time — you just substitute one variable at a time, top to bottom (forward substitution for L) or bottom to top (back substitution for U).

The Two-Step Solve

Given A = LU and Ax = b:

Why is each triangular solve only O(n²)? Because the triangular structure means each equation involves one new unknown and all previously computed values. In forward substitution of an n×n lower triangular system, row i involves a dot product of length i − 1 plus one division. Summing 1 + 2 + … + (n−1) = n(n−1)/2, which is O(n²). Back substitution is symmetric.

Total for the two solves: O(n²). Compare to O(n³/3) for full Gaussian elimination. If you solve 100 systems with the same A, you do the O(n³/3) factorization once and then 100 × O(n²) triangular solves. For n = 1000, that is 3.3 × 10⁸ + 2 × 10⁸ instead of 3.3 × 10¹⁰ — a 10× speedup.

Computing L and U

LU decomposition is just Gaussian elimination with bookkeeping. When you subtract a multiple of row i from row j to create a zero below the pivot, that multiplier goes into L. The resulting upper triangular matrix is U. The diagonal of L is all ones (this is called a unit lower triangular matrix).

Let us walk through a 3×3 example by hand.

Implementation

python
import numpy as np

def lu_decompose(A):
    """Factor A = LU (no pivoting). Returns (L, U)."""
    n = A.shape[0]
    L = np.eye(n)           # start with identity
    U = A.astype(float).copy()

    for k in range(n):
        for i in range(k + 1, n):
            # multiplier: how much of row k to subtract from row i
            L[i, k] = U[i, k] / U[k, k]
            # eliminate: row i -= multiplier * row k
            U[i, k:] -= L[i, k] * U[k, k:]
    return L, U

def forward_sub(L, b):
    """Solve Ly = b for y (L is unit lower triangular)."""
    n = L.shape[0]
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i] - L[i, :i] @ y[:i]
    return y

def back_sub(U, y):
    """Solve Ux = y for x (U is upper triangular)."""
    n = U.shape[0]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
    return x

# Example
A = np.array([[2,3,1],[4,7,5],[6,18,22]], dtype=float)
b = np.array([1, 3, 5], dtype=float)
L, U = lu_decompose(A)
y = forward_sub(L, b)
x = back_sub(U, y)
print("x =", x)            # [ 1.75  -0.5   -0.5 ]
print("Ax =", A @ x)       # [1. 3. 5.] == b ✓

In-Place Storage

A clever implementation trick: L and U can be stored in the space of the original matrix A. Since L has 1s on the diagonal (which we do not need to store) and zeros above, while U has zeros below, they fit together like a jigsaw:

This is how LAPACK (and by extension NumPy/SciPy) stores LU factors. The function scipy.linalg.lu_factor returns a single n×n array plus a permutation vector — not separate L and U matrices. This halves memory usage, which matters for large matrices.

Worked Example: Solving Ax = b via LU

Let us solve the full problem. With A = LU already computed above, solve Ax = b where b = [1, 3, 5].

Notice that the forward and back substitution steps are short — O(n²). The hard work was in the LU factorization. If someone now gives us b = [10, 20, 30], we reuse the same L and U and just do the two substitution steps again.

Chapter 2: Pivoting

The LU decomposition from Chapter 1 has a fatal flaw. What happens when the pivot element u_kk is zero? Division by zero. The algorithm crashes. And even if the pivot is not exactly zero, a very small pivot amplifies rounding errors catastrophically.

The Instability Trap

Consider this 2×2 system using standard floating-point arithmetic with 3 significant digits:

The exact solution is x₁ ≈ 1.001, x₂ ≈ 0.999. Without pivoting, we eliminate using multiplier m = 1/0.001 = 1000. Row 2 becomes: 0·x₁ + (1 − 1000)x₂ = 2 − 1000, so −999 x₂ = −998. In 3-digit arithmetic, −999 rounds to −1000 and −998 rounds to −1000, giving x₂ = 1.000 and then x₁ = 0. The computed x₁ is completely wrong — off by 100%.

With partial pivoting, we swap rows so the largest element in the current column becomes the pivot. Swap rows 1 and 2:

Now the multiplier is m = 0.001/1 = 0.001. Row 2 becomes: 0·x₁ + 0.999 x₂ = 0.998. Even in 3-digit arithmetic: x₂ = 0.998/0.999 = 0.999, x₁ = 2 − 0.999 = 1.001. Both correct to 3 digits.

The difference is dramatic. Without pivoting: x₁ = 0 (100% error). With pivoting: x₁ = 1.001 (0.01% error). Same matrix, same arithmetic precision — the only change was which row we used as the pivot. This is why every serious numerical library (LAPACK, NumPy, MATLAB) uses partial pivoting by default. It is not optional.

The mathematical explanation: the multiplier m = b/a_pivot appears in every subsequent entry. If a_pivot is tiny, m is huge, and the product m · (other entries) has many significant digits that must be subtracted from a number of comparable magnitude. In finite precision, this subtraction causes catastrophic cancellation — the leading digits cancel and you are left with rounding noise. Large pivots keep m small, keeping the arithmetic well-behaved.

The PA = LU Algorithm

At each step k, before eliminating below the pivot:

Worked Example: PA = LU with Partial Pivoting

Factor the matrix B = 1  2  34  5  67  8  0

Full vs. Partial Pivoting

Partial pivoting searches only the current column for the largest element. Complete (full) pivoting searches the entire remaining submatrix — finding the largest element in all remaining rows and columns. Complete pivoting requires column swaps too (PAQ = LU), which makes bookkeeping harder and adds O(n²) comparisons total. In practice, partial pivoting is almost always sufficient. The extra stability of complete pivoting is rarely needed, and the added cost is not justified.

There is a deep theorem behind this: with partial pivoting, all entries of L satisfy |l_ij| ≤ 1. This means that rounding errors in L are bounded. The entries of U can grow (by a factor of at most 2ⁿ⁻¹ in the worst case), but in practice the growth factor is almost always small — typically O(n^2/3) for random matrices. Over 60 years of numerical experience confirms that partial pivoting is robust for virtually all practical problems.

The Permutation Matrix P

A permutation matrix P is the identity matrix with its rows shuffled. Multiplying PA just reorders the rows of A. Key properties:

python
import numpy as np

def lu_partial_pivot(A):
    """PA = LU with partial pivoting. Returns (L, U, P)."""
    n = A.shape[0]
    U = A.astype(float).copy()
    L = np.zeros((n, n))
    P = np.eye(n)

    for k in range(n):
        # Find pivot: row with largest |value| in column k
        pivot_row = k + np.argmax(np.abs(U[k:, k]))

        # Swap rows in U, L, and P
        if pivot_row != k:
            U[[k, pivot_row]] = U[[pivot_row, k]]
            P[[k, pivot_row]] = P[[pivot_row, k]]
            if k > 0:
                L[[k, pivot_row], :k] = L[[pivot_row, k], :k]

        # Eliminate below pivot
        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]

    L += np.eye(n)  # add 1s on diagonal
    return L, U, P

# Solve PA = LU => Ax = b => PAx = Pb => LUx = Pb
A = np.array([[0.001, 1], [1, 1]])
b = np.array([1.0, 2.0])
L, U, P = lu_partial_pivot(A)
Pb = P @ b
y = forward_sub(L, Pb)  # reuse from Ch 1
x = back_sub(U, y)
print(x)  # [1.001001  0.998999] — correct!

Chapter 3: Matrix Inversion

A matrix A is invertible (or nonsingular) if there exists a matrix A⁻¹ such that A · A⁻¹ = A⁻¹ · A = I (the identity matrix). Geometrically, if A rotates and stretches space, then A⁻¹ undoes that rotation and stretching. Not all matrices are invertible — a matrix that squishes 3D space down to a plane has lost information and cannot be undone.

The condition for invertibility: det(A) ≠ 0. Equivalently, all pivots in LU decomposition are nonzero. Equivalently, the rows (or columns) of A are linearly independent.

Computing A⁻¹ via LU

The key idea: A⁻¹ is the matrix X satisfying AX = I. But AX = I is just n linear systems, one for each column of X:

where e_j is the j-th column of the identity matrix (all zeros except a 1 in position j). We already know how to solve Ax = b efficiently using LU. So: compute PA = LU once (cost: n³/3), then solve n systems (each costs n²). Total: n³/3 + n · n² = O(n³).

python
def invert_via_lu(A):
    """Compute A^{-1} using PA = LU decomposition."""
    n = A.shape[0]
    L, U, P = lu_partial_pivot(A)
    A_inv = np.zeros((n, n))

    for j in range(n):
        # Solve Ax_j = e_j  =>  LUx_j = P e_j
        e_j = np.zeros(n)
        e_j[j] = 1.0
        Pe = P @ e_j
        y = forward_sub(L, Pe)
        A_inv[:, j] = back_sub(U, y)

    return A_inv

A = np.array([[2,3,1],[4,7,5],[6,18,22]], dtype=float)
A_inv = invert_via_lu(A)
print(A @ A_inv)  # ≈ identity matrix

Worked Example: Inverting a 3×3 Matrix via LU

Using our familiar A = 2  3  14  7  56  18  22 with L and U from Chapter 1:

Assembling: A⁻¹ = −4    3    −0.53.625   −2.375   0.375−1.875   1.125   −0.125

You can verify: A · A⁻¹ = I (try the first row: 2(−4)+3(3.625)+1(−1.875) = −8+10.875−1.875 = 1 ✓).

When Do You Actually Need A⁻¹?

Almost never for solving Ax = b. But there are genuine use cases:

Matrix Inversion and Strassen

The Determinant for Free

Once you have A = LU (or PA = LU), computing the determinant is trivial:

The determinant of a triangular matrix is the product of its diagonal entries. L has all 1s on the diagonal, so det(L) = 1. U has the pivots on its diagonal. The sign depends on whether P involves an even or odd number of row swaps.

For our example: det(A) = det(U) = 2 · 1 · (−8) = −16. You can verify with the Sarrus rule for 3×3 matrices: 2(7·22 − 5·18) − 3(4·22 − 5·6) + 1(4·18 − 7·6) = 2(154−90) − 3(88−30) + 1(72−42) = 128 − 174 + 30 = −16 ✓.

This is how all serious software computes determinants. The naïve cofactor expansion has cost O(n!), which is utterly impractical for n > 15. LU gives det in O(n³/3) — the same cost as solving a system.

A beautiful result from CLRS: matrix inversion, matrix multiplication, and LU decomposition are all computationally equivalent up to constant factors. If you can multiply n×n matrices in O(n^ω) time (where ω < 3 for Strassen and beyond), then you can also invert and factor in O(n^ω) time. Currently the best known ω ≈ 2.371, so theoretically all three can be done in O(n^2.371). In practice, the constants are huge, so O(n³) LU is used for matrices smaller than about 10,000 × 10,000.

The proof of equivalence is elegant. To show "multiply ≤ invert": given two n×n matrices A, B, build the 3n×3n matrix M = I   A   00   I   B0   0   I. Its inverse has AB in the (1,3) block. So one matrix inversion gives us a matrix product. The other direction ("invert ≤ multiply") uses a divide-and-conquer argument, partitioning the matrix into blocks and recursing. This mutual reduction proves they have the same asymptotic complexity.

Chapter 4: Symmetric Positive Definite & Cholesky

A special class of matrices appears so frequently in science and engineering that it deserves its own decomposition. A matrix A is symmetric positive definite (SPD) if:

Where do SPD matrices come from? Everywhere:

Cholesky Decomposition

For an SPD matrix, LU decomposition can be simplified. Because A is symmetric, U turns out to be a scaled version of L^T. Specifically, A = LL^T where L is lower triangular with positive diagonal entries. This is the Cholesky decomposition.

The cost: n³/6 — half the cost of general LU. This is because symmetry means we only need to compute the lower triangle; the upper triangle is just the transpose.

Why SPD Matrices Are Special

SPD matrices have remarkable properties that make them the "nicest" matrices in linear algebra:

The no-pivoting property is the most practical advantage. For general matrices, you must pivot to avoid numerical catastrophe. For SPD matrices, the pivots are naturally positive and well-behaved — Cholesky is backward stable without any row swaps. This is why Cholesky is the default solver in optimization (positive definite Hessians), statistics (covariance matrices), and Kalman filtering (covariance updates).

Generating SPD Matrices

A common interview/code question: how do you create an SPD matrix for testing? The standard trick: take any matrix B with full column rank, then A = B^TB is guaranteed SPD.

python
# Generate a random SPD matrix
B = np.random.randn(10, 5)  # any 10×5 matrix
A = B.T @ B                    # 5×5, guaranteed SPD

# Alternative: start from eigenvalues
Q, _ = np.linalg.qr(np.random.randn(5, 5))  # random orthogonal
eigvals = np.array([10, 5, 2, 1, 0.1])   # choose condition number
A = Q @ np.diag(eigvals) @ Q.T                # κ = 10/0.1 = 100

The Algorithm

Derive L row by row. For each row i, we know that A = LL^T, so:

Solving for the entries of L:

python
def cholesky(A):
    """Cholesky decomposition A = LL^T for SPD matrix A."""
    n = A.shape[0]
    L = np.zeros((n, n))

    for i in range(n):
        for j in range(i + 1):
            s = np.sum(L[i, :j] * L[j, :j])
            if i == j:
                # Diagonal: square root
                L[i, j] = np.sqrt(A[i, i] - s)
            else:
                # Below diagonal
                L[i, j] = (A[i, j] - s) / L[j, j]
    return L

# Example: covariance matrix (always SPD)
A = np.array([[4, 2, 2],
              [2, 5, 3],
              [2, 3, 6]], dtype=float)
L = cholesky(A)
print("L =\n", L)
print("LL^T =\n", L @ L.T)  # == A ✓

# Solve Ax = b using Cholesky
b = np.array([1.0, 2.0, 3.0])
y = forward_sub(L, b)
x = back_sub(L.T, y)
print("x =", x)

Worked Example: 3×3 Cholesky by Hand

Decompose A = 4  2  22  5  32  3  6

Notice every diagonal entry of L is 2 — because the diagonal of A is [4, 5, 6] and the accumulated sum of squares at each step leaves exactly 4 under each square root. This is not a coincidence for this matrix but illustrates how the algorithm peels off one layer of "variance" at a time.

Cholesky in NumPy and SciPy

In practice, you almost never implement Cholesky yourself. The standard tools:

python
import numpy as np
from scipy import linalg

A = np.array([[4,2,2],[2,5,3],[2,3,6]], dtype=float)

# NumPy: returns lower triangular L
L = np.linalg.cholesky(A)

# SciPy: returns upper triangular by default!
U = linalg.cholesky(A)       # upper: A = U^T U
L2 = linalg.cholesky(A, lower=True)  # lower: A = L L^T

# Solve Ax = b using Cholesky (SciPy one-liner)
b = np.array([1.0, 2.0, 3.0])
x = linalg.cho_solve(linalg.cho_factor(A), b)
print("x =", x)
print("Ax =", A @ x)  # == b ✓

Chapter 5: Least Squares

Not every system Ax = b has a solution. If you have 100 data points and want to fit a line (2 parameters), you have 100 equations and 2 unknowns — the system is overdetermined. There is no x that satisfies all 100 equations exactly. But you can find the x that comes closest.

The Least Squares Problem

Given an m×n matrix A (m > n) and a vector b ∈ R^m, find x ∈ Rⁿ that minimizes the squared residual:

Each term (a_i^Tx − b_i)² is the squared error between the model's prediction a_i^Tx and the observed value b_i. Minimizing the sum of squared errors is the principle behind linear regression — the foundation of statistics and ML.

The Normal Equations

To minimize ||Ax − b||², take the derivative with respect to x and set it to zero:

These are the normal equations. The matrix A^TA is n×n (small!), symmetric, and positive semidefinite. If A has full column rank (columns are independent), then A^TA is positive definite — perfect for Cholesky. The solution pipeline:

Total cost: O(mn²) — dominated by forming A^TA. For 10,000 data points and 10 features, that is about 10⁶ operations. Fast.

Where Does ||Ax − b||² Come From?

Why minimize the squared error rather than, say, the absolute error |Ax − b|? Three reasons: (1) The squared error is differentiable everywhere, so calculus gives us a clean closed-form solution. Absolute error has a kink at zero, making optimization harder. (2) Under Gaussian noise assumptions, minimizing squared error gives the maximum likelihood estimate — the statistically optimal estimate. (3) The resulting normal equations are a symmetric positive (semi)definite system, which we can solve efficiently with Cholesky.

If the noise is not Gaussian (e.g., heavy-tailed outliers), least squares can be fooled. A single extreme outlier pulls the fit toward it. Robust alternatives include least absolute deviation (L1 regression, solved by linear programming — Ch 29) and Huber loss (a hybrid that is quadratic for small residuals and linear for large ones).

Geometric Interpretation

The column space of A is the set of all vectors Ax for any x. When Ax = b has no solution, b is not in the column space. The least squares solution finds the point in the column space closest to b — the orthogonal projection. The residual r = b − Ax̂ is perpendicular to the column space, which is exactly what A^Tr = 0 says (the normal equations).

Think of it this way: you are standing at point b in a high-dimensional room. The column space of A is a flat surface (a subspace). The least squares solution is the shadow of b cast straight down onto that surface. The normal equations guarantee the shadow is directly below you — the residual (the vertical drop) is perpendicular to the surface.

Worked Example: Linear Regression by Hand

Fit y = c₀ + c₁t to data: (1,2), (2,3), (3,6), (4,8).

Polynomial Fitting

To fit a polynomial of degree d to m data points (t₁,y₁), …, (t_m,y_m), construct the Vandermonde matrix:

The coefficient vector x = [c₀, c₁, …, c_d] gives the polynomial p(t) = c₀ + c₁t + c₂t² + … + c_dt^d that best fits the data in the least-squares sense.

python
def least_squares_poly(t, y, degree):
    """Fit polynomial of given degree to data (t, y)."""
    # Build Vandermonde matrix
    A = np.column_stack([t**k for k in range(degree + 1)])

    # Normal equations: (A^T A) x = A^T y
    ATA = A.T @ A
    ATy = A.T @ y

    # Solve via Cholesky
    L = cholesky(ATA)
    z = forward_sub(L, ATy)
    coeffs = back_sub(L.T, z)
    return coeffs

# Example: fit a quadratic to noisy data
t = np.linspace(0, 5, 50)
y_true = 1 + 2*t - 0.5*t**2
y_noisy = y_true + np.random.randn(50) * 0.5

coeffs = least_squares_poly(t, y_noisy, degree=2)
print("Coefficients:", coeffs)
# ≈ [1.0, 2.0, -0.5] — recovers the true polynomial

The Overfitting Trap

In the simulation above, try increasing the polynomial degree to 8 or higher with only 10 data points. The curve passes through (or very close to) every point — the residuals drop to zero. But the curve oscillates wildly between the points. This is overfitting: the polynomial is fitting the noise, not the underlying pattern.

Overfitting happens when the model has more parameters than the data can constrain. A degree-d polynomial has d+1 parameters. With m data points, you need m >> d+1 to get a reliable fit. The rule of thumb in statistics is at least 10-20 observations per parameter.

This is directly connected to the condition number. As the polynomial degree approaches the number of data points, the Vandermonde matrix becomes increasingly ill-conditioned. The normal equations A^TA become nearly singular, and tiny changes in the data produce wild swings in the coefficients. Regularization (ridge regression, from Chapter 6) controls overfitting by penalizing large coefficients.

Chapter 6: Numerical Stability

Everything we have done so far assumes exact arithmetic. Real computers use floating-point numbers with roughly 16 significant digits (64-bit doubles). Rounding errors are tiny — but they can be amplified by the structure of the matrix. Understanding when and why this happens is the difference between numerical code that works and code that silently produces garbage.

Condition Number

The condition number κ(A) measures how sensitive Ax = b is to perturbations. If you perturb b by a tiny amount δb, the solution changes by δx, and:

If κ(A) = 10^k, then you lose roughly k digits of accuracy. With 16-digit doubles and κ = 10¹², you get only 4 correct digits.

Think of it this way: the condition number tells you how many digits of your input get "burned up" by the problem. You start with 16 digits of precision. The problem burns κ of them. You are left with 16 − log₁₀(κ) correct digits. If κ uses up all 16, you get zero correct digits — the output is pure noise.

The condition number is defined as:

where σ_max and σ_min are the largest and smallest singular values of A. A matrix with κ ≈ 1 is well-conditioned; a matrix with κ >> 1 is ill-conditioned.

Some common condition numbers to calibrate your intuition:

QR Decomposition

For least squares, the normal equations A^TAx = A^Tb square the condition number: κ(A^TA) = κ(A)². If A has κ = 10⁶, then A^TA has κ = 10¹², and you lose 12 digits — nearly all your precision.

The fix: QR decomposition. Factor A = QR where Q has orthonormal columns (Q^TQ = I) and R is upper triangular. Then:

(using the fact that orthogonal transformations preserve norms). So we just solve the triangular system Rx = Q^Tb. No squaring of the condition number. QR costs about twice as much as Cholesky for least squares (2mn² vs mn²), but it is far more numerically stable.

Worked Example: Condition Number Disaster

Consider the 2×2 matrix:

This matrix is invertible (det = 0.0001), but its condition number is κ ≈ 40,001. Solving Ax = [2, 2.0001] gives x = [1, 1]. But solving Ax = [2, 2.0002] (perturbing b by 0.00005%) gives x = [0, 2] — a 50% change in x₁. A tiny perturbation in b caused a massive change in x. This is what ill-conditioning looks like.

Now form A^TA for least squares:

This matrix has κ ≈ 1.6 × 10⁹ — the condition number squared. With 16-digit doubles, you get about 7 correct digits. With Cholesky on A^TA, the answer might be okay but barely. With QR directly on A, you get about 12 correct digits. That is the difference between "usable" and "pushing the limits."

SVD: The Gold Standard

The singular value decomposition A = UΣV^T decomposes any matrix into a rotation U, a scaling Σ (diagonal), and another rotation V^T. It is the most stable decomposition, handles rank-deficient matrices gracefully, and directly reveals the condition number (σ_max/σ_min). The cost is O(mn²) — similar to QR. NumPy's np.linalg.lstsq uses SVD internally.

SVD is the Swiss Army knife of numerical linear algebra. It solves least squares, computes pseudoinverses (for non-square or rank-deficient matrices), performs PCA (principal component analysis), and is the basis for low-rank approximation (the best rank-k approximation of a matrix is given by truncating its SVD — the Eckart-Young theorem).

The tradeoff: SVD is about 10× more expensive than Cholesky for a well-conditioned least squares problem. Use Cholesky when you know the matrix is well-conditioned. Use QR for moderate ill-conditioning. Use SVD when you do not trust the data or when the matrix may be rank-deficient.

Practical Decision Tree for Least Squares

Regularization: Taming Ill-Conditioning

When κ(A) is large, you can improve the condition number by adding regularization. Instead of solving A^TAx = A^Tb, solve (A^TA + λI)x = A^Tb. This is ridge regression (also called Tikhonov regularization).

Why does it work? Adding λI to A^TA increases every eigenvalue by λ. If the smallest eigenvalue was σ_min², it becomes σ_min² + λ. The condition number drops from σ_max²/σ_min² to (σ_max² + λ)/(σ_min² + λ). For large λ, this approaches 1.

For example, with λ = 0.1 and a matrix that has σ_min² = 10⁻¹⁰ and σ_max² = 1, the condition number drops from 10¹⁰ to (1 + 0.1)/(10⁻¹⁰ + 0.1) ≈ 11. That is a nine-order-of-magnitude improvement in conditioning, at the cost of a small bias in the solution.

The tradeoff: regularization introduces bias. The solution is no longer the true least squares answer — it is shrunk toward zero. But the reduction in variance (from better conditioning) often more than compensates. This is the bias-variance tradeoff, and ridge regression is one of the most important tools in statistical machine learning.

python
# Three ways to solve least squares in NumPy
import numpy as np

A = np.random.randn(100, 5)
b = np.random.randn(100)

# Method 1: Normal equations (fastest, least stable)
x1 = np.linalg.solve(A.T @ A, A.T @ b)

# Method 2: QR decomposition (good balance)
Q, R = np.linalg.qr(A)
x2 = np.linalg.solve(R, Q.T @ b)

# Method 3: SVD (most stable, recommended)
x3, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)

# Check condition number
print("κ(A) =", np.linalg.cond(A))
print("κ(AᵀA) =", np.linalg.cond(A.T @ A))  # ≈ κ(A)²

python
# Ridge regression: regularized least squares
def ridge_regression(X, y, lam):
    """Solve (X^T X + λI) w = X^T y"""
    n = X.shape[1]
    ATA = X.T @ X + lam * np.eye(n)
    ATy = X.T @ y
    return np.linalg.solve(ATA, ATy)

# Compare condition numbers
X = np.random.randn(100, 10)
# Make it ill-conditioned by adding near-duplicate column
X[:, 9] = X[:, 0] + 1e-6 * np.random.randn(100)

print("κ(X^T X)        =", np.linalg.cond(X.T @ X))
print("κ(X^T X + 0.1I) =", np.linalg.cond(X.T @ X + 0.1*np.eye(10)))
# Condition drops by orders of magnitude!

Chapter 7: Applications

Matrix operations are not abstract mathematics — they are the computational backbone of modern technology. Every field that involves data, physics, or optimization eventually reduces to solving linear systems or computing decompositions. Here is where each technique from this chapter shows up in the real world.

Machine Learning: Linear Regression

The most direct application. Given a dataset of m samples with n features and a target variable, linear regression finds weights w that minimize ||Xw − y||². This is exactly the least squares problem from Chapter 5. Scikit-learn's LinearRegression solves the normal equations (or uses SVD for numerical stability). Every ML engineer uses this daily, even when the model is a neural network — because the loss function still bottoms out at a system of linear equations in the final layer.

Computer Graphics: Transformations

Every vertex in a 3D scene is transformed by a 4×4 matrix: rotation, translation, scaling, projection. A typical rendering pipeline applies 3-5 matrix multiplications per vertex, and a single frame has millions of vertices. GPUs are essentially massively parallel matrix multiplication engines. The model-view-projection matrix M = P · V · M multiplies three 4×4 matrices — and the inverse M⁻¹ is needed for ray casting and picking.

Physics: Finite Element Methods

Simulating how a bridge bends, how heat flows through a wall, or how an airplane wing vibrates all reduce to solving Kx = f, where K is a large sparse SPD stiffness matrix and f is the applied force. Cholesky decomposition (or iterative solvers for very large systems) is the workhorse. Modern structural simulations involve matrices with millions of rows — but they are sparse (mostly zeros), so specialized sparse Cholesky implementations handle them efficiently.

Signal Processing: Fourier and Wavelets

The Discrete Fourier Transform is a matrix-vector product y = Fx where F is the Fourier matrix. The FFT is an O(n log n) factorization of this matrix into sparse factors — essentially an LU-like decomposition that exploits the structure of F. Wavelets, cosine transforms, and Hadamard transforms are all matrix factorizations.

PageRank and Recommender Systems

Google's PageRank finds the dominant eigenvector of the web's link matrix — a matrix operation. Netflix's recommendation engine uses matrix factorization: approximate the user-item rating matrix as a product of two low-rank matrices. Both are solved by iterative methods that reduce to repeated matrix-vector products and linear solves.

Neural Networks: Backpropagation

Every layer of a neural network computes y = Wx + b (a matrix-vector product) followed by a nonlinearity. Backpropagation computes gradients by multiplying Jacobian matrices — chain rule is matrix multiplication. The reason GPUs accelerate deep learning is that training is dominated by matrix multiplications. Modern transformers (GPT, etc.) spend most of their compute on the attention matrix Q K^T/√d and the projection matrices WQ, WK, WV, WO.

Kalman Filters and Robotics

The Kalman filter — the workhorse of state estimation in robotics, navigation, and aerospace — performs a matrix inversion (or LU solve) at every timestep. The innovation covariance S = H P H^T + R must be inverted to compute the Kalman gain K = P H^T S⁻¹. For a robot tracking 100 landmarks in SLAM, the state vector is 200+ dimensions, and the Cholesky decomposition of the covariance matrix dominates the computation. Efficient matrix operations are the difference between real-time and too-slow.

Economics: Input-Output Models

Leontief's input-output model describes an economy where each industry consumes outputs of other industries. The consumption matrix C says how much of industry j's output is needed per unit of industry i's output. If d is the external demand vector, the total production x must satisfy x = Cx + d, or (I − C)x = d. This is a linear system, solved by LU decomposition. Leontief won the 1973 Nobel Prize in Economics for this model — and the computational technique behind it is exactly what we learned in Chapter 1.

Cryptography and Number Theory

Lattice-based cryptography (the basis for post-quantum encryption schemes like CRYSTALS-Kyber) uses matrix operations over finite fields. Finding short vectors in a lattice — the hard problem underlying the security — is related to LU decomposition and the Gram-Schmidt process (which is QR decomposition). The LLL algorithm for lattice basis reduction is essentially iterated QR factorization.

Computational Biology

In genomics, principal component analysis (PCA) on a genotype matrix (millions of SNPs × thousands of individuals) reveals population structure. PCA is computed via SVD — a matrix decomposition. Similarly, gene expression analysis uses least squares to fit linear models relating expression levels to experimental conditions. The entire field of bioinformatics is built on the matrix operations from this chapter.

Finance: Portfolio Optimization

Markowitz mean-variance optimization finds the portfolio weights w that minimize variance w^TΣw subject to a target return. The covariance matrix Σ is SPD, so the optimal weights come from solving a linear system using Cholesky. Risk modeling, option pricing via finite differences, and credit scoring all reduce to matrix operations.

Control Systems

Designing autopilots, robot controllers, and industrial regulators requires solving Riccati equations — matrix equations of the form A^TXA − X − A^TXB(R + B^TXB)⁻¹B^TXA + Q = 0. Each step of the iterative solver involves an LU solve (for the matrix inverse) and matrix multiplications. The linear-quadratic regulator (LQR) — the foundation of optimal control — is solved by iterating this equation to convergence. Every drone, self-driving car, and rocket engine controller depends on these matrix operations running fast enough for real-time control.

The Computational Bottleneck

A recurring theme: the most expensive part of many algorithms is a matrix operation. Interior-point methods for linear programming (Ch 29) solve a system of linear equations per iteration. The Hessian solve in Newton's method for optimization is an LU or Cholesky decomposition. Even randomized algorithms like stochastic gradient descent can be seen as approximate matrix operations. Understanding the cost of these operations tells you the fundamental speed limit of many algorithms.

Chapter 8: Interview Arsenal

Matrix operations appear in system design interviews (ML pipelines, recommendation engines), coding interviews (implement linear regression from scratch), and theory interviews (complexity of matrix operations). Here is your cheat sheet.

Decomposition Comparison

Complexity Cheat Sheet

Coding Drills

Implement each of these from scratch (no numpy.linalg). They are common in ML and systems interviews.

python
# Drill 1: LU decomposition with partial pivoting
# Input: n×n matrix A
# Output: L, U, perm (permutation vector)

def lu_pivot(A):
    n = len(A)
    U = [row[:] for row in A]  # deep copy, pure Python
    L = [[0.0]*n for _ in range(n)]
    perm = list(range(n))

    for k in range(n):
        # Find pivot
        max_val, max_row = abs(U[k][k]), k
        for i in range(k+1, n):
            if abs(U[i][k]) > max_val:
                max_val, max_row = abs(U[i][k]), i
        # Swap
        U[k], U[max_row] = U[max_row], U[k]
        L[k], L[max_row] = L[max_row], L[k]
        perm[k], perm[max_row] = perm[max_row], perm[k]
        # Eliminate
        for i in range(k+1, n):
            L[i][k] = U[i][k] / U[k][k]
            for j in range(k, n):
                U[i][j] -= L[i][k] * U[k][j]
    for i in range(n):
        L[i][i] = 1.0
    return L, U, perm

python
# Drill 2: Linear regression from scratch
# Input: X (m×n features), y (m targets)
# Output: weights w (n,)

def linear_regression(X, y):
    """Solve normal equations: (X^T X) w = X^T y"""
    XTX = [[0.0]*len(X[0]) for _ in range(len(X[0]))]
    XTy = [0.0]*len(X[0])
    m, n = len(X), len(X[0])

    # Form X^T X (n×n) and X^T y (n,)
    for i in range(n):
        for j in range(n):
            for k in range(m):
                XTX[i][j] += X[k][i] * X[k][j]
        for k in range(m):
            XTy[i] += X[k][i] * y[k]

    # Solve using LU
    L, U, perm = lu_pivot(XTX)
    Pb = [XTy[perm[i]] for i in range(n)]
    # Forward sub
    fwd = [0.0]*n
    for i in range(n):
        fwd[i] = Pb[i] - sum(L[i][k]*fwd[k] for k in range(i))
    # Back sub
    w = [0.0]*n
    for i in range(n-1, -1, -1):
        w[i] = (fwd[i] - sum(U[i][k]*w[k] for k in range(i+1,n))) / U[i][i]
    return w

python
# Drill 3: Cholesky decomposition (no numpy)
import math

def cholesky_pure(A):
    """Cholesky A = LL^T for SPD matrix. Pure Python."""
    n = len(A)
    L = [[0.0]*n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = sum(L[i][k] * L[j][k] for k in range(j))
            if i == j:
                val = A[i][i] - s
                if val <= 0:
                    raise ValueError("Not positive definite")
                L[i][j] = math.sqrt(val)
            else:
                L[i][j] = (A[i][j] - s) / L[j][j]
    return L

# Test with SPD matrix
A = [[4,2,2],[2,5,3],[2,3,6]]
L = cholesky_pure(A)
# Verify: L[0] = [2, 0, 0], L[1] = [1, 2, 0], L[2] = [1, 1, 2]

python
# Drill 4: Solve triangular system (interview classic)

def solve_upper_triangular(U, b):
    """Solve Ux = b for upper triangular U. Pure Python."""
    n = len(b)
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - sum(U[i][j]*x[j] for j in range(i+1, n))) / U[i][i]
    return x

# Drill 5: Compute determinant via LU
def det_via_lu(A):
    """det(A) = det(L) * det(U) * det(P).
    det(L) = 1 (unit diagonal). det(U) = product of diagonal.
    det(P) = (-1)^(number of swaps)."""
    L, U, perm = lu_pivot(A)
    # Count swaps in permutation
    n = len(perm)
    visited = [False] * n
    swaps = 0
    for i in range(n):
        if visited[i]: continue
        j, cycle_len = i, 0
        while not visited[j]:
            visited[j] = True
            j = perm[j]; cycle_len += 1
        swaps += cycle_len - 1

    det_U = 1.0
    for i in range(n):
        det_U *= U[i][i]
    return ((-1)**swaps) * det_U

Key Formulas to Memorize

Common Interview Pitfalls

Interview Questions & Answers

Complexity Genealogy

The relationships between matrix operations form a hierarchy that every computer scientist should know:

Chapter 9: Connections

Matrix operations are not an isolated topic — they are the computational substrate that other algorithms build upon. Here is how Chapter 28 connects to the rest of CLRS, other micro-lessons, and the broader world.

Within CLRS

Within This Site

The Decomposition Hierarchy

Here is how the decompositions relate — think of it as a decision tree:

Limitations

This chapter covers dense matrix operations — matrices where most entries are nonzero. Many real-world matrices (graphs, finite elements, sparse data) have mostly zero entries. Sparse solvers exploit this structure for massive speedups (storing only nonzero entries, using iterative methods like conjugate gradient instead of direct factorization). For truly large systems (millions of unknowns), iterative methods replace direct decompositions entirely.

The numerical stability discussion here barely scratches the surface. A full treatment requires understanding floating-point arithmetic (IEEE 754), backward error analysis, and the nuances of different pivot strategies. See Trefethen & Bau, Numerical Linear Algebra, for the definitive treatment.

What We Did Not Cover

Summary: The Big Ideas

If you remember nothing else from this chapter, remember these five principles:

Further Reading

"The purpose of computing is insight, not numbers."
— Richard Hamming