Linear Algebra with NumPy & SciPy

Systems of linear equations, as pictures

A classic puzzle: "Two coffees and one muffin cost $7; one coffee and three muffins cost $11. What does each cost?" Calling the coffee price x and the muffin price y:

2x +  y =  7
 x + 3y = 11

Each equation is a straight line in the (x, y) plane: all the price pairs that satisfy it. The answer is the point that lies on both lines, where they cross. There are only three possible pictures for two lines:

Stack the coefficients into a matrix and the unknowns into a vector, and the whole system becomes one matrix-vector equation, using the multiplication from M47:

A = [[2, 1],     x = [x,     b = [ 7,        A @ x == b
     [1, 3]]          y]          11]

There's also a column picture: A x = x·(2, 1) + y·(1, 3), so solving means "find how much of each column to mix to reach b". That reading is the one that generalizes to least squares later.

Step through a matrix product below to refresh the "row dot column" rule that A @ x uses.

Remember: each linear equation is a line (or plane); the solution is where they all meet; A x = b packages the whole system.

solve, not inverse

For numbers, you solve 3x = 12 by dividing: x = 12 / 3. Matrices have an analogue of "divide by": the inverse A⁻¹, the matrix that undoes A (A⁻¹A = I). Then x = A⁻¹ b. Textbooks write that; working code almost never computes it. Use np.linalg.solve(A, b) instead, for three reasons:

• Speed: solve uses Gaussian elimination (the "subtract a multiple of one equation from another" method you may have done by hand), packaged as an LU factorization. Computing the full inverse does about three times as much work, then you still multiply.
• Accuracy: every extra floating-point operation adds rounding error. solve typically gives a smaller residual |A x - b|, and the difference becomes dramatic for nearly-singular matrices.
• Structure: the inverse of a sparse matrix is usually dense, which can turn megabytes into terabytes.

The condition number np.linalg.cond(A) measures how sensitive the answer is to tiny changes in A or b. Roughly, you lose about log10(cond) of your ~16 significant digits. A condition number of 10¹² leaves about 4 trustworthy digits, whichever method you use.

Remember: x = np.linalg.solve(A, b). Treat inv as a code smell, and check cond when results look odd.

Determinant and rank: does the matrix squash space?

Recall from M47 that a 2 × 2 matrix transforms the plane, sending the unit square (area 1) to a parallelogram whose sides are the matrix's columns. The determinant det(A) is the area of that parallelogram: the factor by which the transformation scales areas. For [[a, b], [c, d]] it is ad - bc. (In 3-D it scales volumes.)

• |det| > 1 expands areas; |det| < 1 shrinks them; rotations and shears have det = 1.
• A negative determinant means the transformation also flips orientation (a mirror image, like the reflected "F" in M47).
• det = 0 means the matrix squashes the plane onto a line (or a point). Many different inputs land on the same output, so you can't undo it: the matrix is singular (not invertible), and A x = b has either no solution or infinitely many. These are the parallel and same-line pictures above.

Rank is the number of genuinely independent directions the output fills: the dimension of the space you can reach. A 3 × 3 matrix with rank 3 fills all of space; rank 2 flattens it to a plane; rank 1 to a line. Equivalently, it's the number of linearly independent columns (columns that can't be made by mixing the others). Rank is the honest version of "is it singular?", because np.linalg.matrix_rank uses a tolerance, while a determinant can be tiny for perfectly healthy matrices (det(0.1 · I₁₀₀) = 10⁻¹⁰⁰).

Remember: det = area (volume) scale factor, sign = orientation flip, zero = squashed and non-invertible. Rank counts independent directions; prefer it to "det == 0" tests.

Norms: different ways to measure size

M47's length √(Σ xᵢ²) is one way to measure a vector's size. A norm is any such size measure that behaves sensibly (zero only for the zero vector, scales with the vector, and obeys the triangle inequality |v + w| ≤ |v| + |w|). Three norms do most of the work:

Remember: L2 for distance, L1 for robustness and sparsity, L∞ for worst-case error.

Projections and least squares

Projection: the closest point on a line

In M47 you saw the "shadow" of w on v. The shadow, as a vector, is the projection of w onto the line through v:

proj = (w @ v) / (v @ v) * v          # a multiple of v
residual = w - proj                   # what's left over; perpendicular to v

The projection is the point on the line closest to w, and the leftover (residual) is perpendicular to the line. "Closest point = drop a perpendicular" is the key geometric fact behind least squares.

Least squares: when there are more equations than unknowns

Suppose you measured study hours and exam scores for 30 students and want a line score ≈ m·hours + c. That's 30 equations (one per student) but only 2 unknowns (m and c). Real data is noisy, so no line passes through every point: the system is overdetermined and has no exact solution. Least squares picks the line that minimizes the sum of squared vertical errors Σ (yᵢ - (m xᵢ + c))².

In matrix form, the design matrix X has a column of x-values and a column of 1s, and we want X @ [m, c] ≈ y. Geometrically, the columns of X span a flat plane inside 30-dimensional space; y sticks out of it; the best fit is the projection of y onto that plane. np.linalg.lstsq(X, y) computes it stably (via the SVD, below). You may see the "normal equations" (XᵀX) β = Xᵀy in textbooks; they give the same answer but square the condition number, so prefer lstsq.

Remember: projection = closest point, residual is perpendicular. Least squares projects y onto the column space of X; use np.linalg.lstsq.

Eigenvalues and eigenvectors: directions that only stretch

Apply a matrix to most vectors and they change direction. But some special directions are only stretched (or shrunk, or flipped): the output points along the same line as the input. Those are eigenvectors ("eigen" is German for "own" or "characteristic"), and the stretch factor is the eigenvalue λ:

A @ v == λ * v        # v is an eigenvector, λ its eigenvalue

For A = [[2, 1], [1, 2]], the direction (1, 1) is stretched by 3, and (1, -1) is left alone (stretch 1). Every other direction gets pulled toward (1, 1), the direction with the largest eigenvalue. The plot below shows it: gray arrows are unit vectors, colored arrows are where A sends them, and the eigen-directions are the only ones that don't turn.

Why care? Eigenvectors reveal what a matrix "does" in its most natural coordinates:

• Repeated application: applying A many times makes the largest-eigenvalue direction dominate. That's power iteration, and it's how PageRank ranks web pages.
• Stability: in a system that evolves by xₙ₊₁ = A xₙ, eigenvalues with |λ| > 1 mean something blows up; all |λ| < 1 means it settles down.
• Data: the eigenvectors of a covariance matrix are the principal directions of the data (PCA, below).

For symmetric matrices (A == A.T, like covariance matrices), the eigenvalues are real and the eigenvectors are perpendicular. Use np.linalg.eigh for them: it is faster, more accurate, and returns eigenvalues sorted ascending. A general matrix can have complex eigenvalues: a rotation, for example, turns every direction, so it has no real eigenvectors.

Remember: A v = λ v: eigenvectors keep their direction, eigenvalues are the stretch. Use eigh for symmetric matrices.

The SVD and low-rank approximation

Eigenvectors only work for square matrices, and not even all of those. The singular value decomposition (SVD) works for every matrix, any shape. It says that any linear transformation is, in the right coordinates, just three simple steps:

A = U @ diag(s) @ Vt
#   rotate  ·  stretch along axes  ·  rotate
#   U, Vt: orthogonal (pure rotations/reflections)   s: singular values, s[0] ≥ s[1] ≥ … ≥ 0

The SVD can also be read as a sum of simple pieces: A = s₀ u₀ v₀ᵀ + s₁ u₁ v₁ᵀ + …, where each uᵢ vᵢᵀ is a rank-1 matrix (an outer product, as in M47). The pieces are ordered by importance. Keeping only the first k gives the best possible rank-k approximation of A (the Eckart–Young theorem): no other rank-k matrix is closer. That is the idea behind image compression, noise removal, recommender systems, and PCA.

Storage tells the story. An m × n image needs m·n numbers; a rank-k approximation needs only k·(m + n + 1). For a 200 × 300 image and k = 20, that's 60,000 vs 10,020: about 6× smaller.

The singular-value plot is the diagnostic: a few big values followed by a long flat tail means the image (or dataset) is "really" low-rank plus noise. Cutting the tail compresses and denoises.

Remember: A = U diag(s) Vᵀ for any matrix: rotate, stretch, rotate. Truncating to the top k singular values gives the best rank-k approximation.

PCA from scratch via the SVD

Principal component analysis finds the directions along which data varies the most. Picture a cloud of points shaped like a tilted cigar: the first principal component points along the cigar, the second across it. Projecting onto the first few components keeps most of the information with fewer numbers, which is useful for visualization (squash 50 features into a 2-D scatter plot), compression, and speeding up models.

The recipe, with data X of shape (n_samples, n_features):

1. Center: subtract each column's mean, so the cloud sits at the origin (broadcasting from M48).
2. SVD: U, s, Vt = np.linalg.svd(Xc, full_matrices=False). The rows of Vt are the principal directions, already sorted by importance.
3. Variance explained: the variance along component i is s[i]**2 / (n - 1); divide by the total for the explained variance ratio.
4. Project: Z = Xc @ Vt[:k].T gives the new k-dimensional coordinates. Reconstruct approximately with Z @ Vt[:k] + mean.

(Equivalently, the principal directions are the eigenvectors of the covariance matrix XcᵀXc / (n - 1); the SVD gets them without forming that matrix, which is more accurate.)

Remember: center, SVD, keep the top rows of Vt; s²/(n-1) are the variances; project with Xc @ Vt[:k].T.

scipy.linalg: factorizations you can reuse

scipy.linalg contains everything in np.linalg plus more, always built against LAPACK. Its main extra value is factorizations: splitting a matrix into simpler pieces once, then reusing them.

• LU: A = P L U, where L is lower-triangular, U is upper-triangular, and P is a row permutation. This is Gaussian elimination written down. lu_factor/lu_solve reuse it for many right-hand sides.
• Cholesky: for symmetric positive-definite matrices (covariance matrices, many physics problems), A = L Lᵀ, about twice as fast as LU. Failure (LinAlgError) is also a cheap test that a matrix is not positive-definite.
• Triangular solves: a triangular system is solved by substitution in O(n²) (one unknown at a time, top to bottom). solve_triangular does it; solving LUx = b is two of these.

Remember: factor once, solve many times. LU for general square systems, Cholesky for symmetric positive-definite, triangular solves are O(n²).

Sparse matrices: storing only what's there

Many real matrices are almost entirely zeros. A friendship graph of 1 million people where each has about 200 friends has 10¹² possible entries but only 2 × 10⁸ ones: 99.98% zeros. Stored densely as float64, that's 8 TB; stored sparsely, a few gigabytes. Machine-learning feature matrices are similar: one-hot encoded categories and bag-of-words text produce rows with a handful of non-zeros among tens of thousands of columns. A sparse matrix stores only the non-zero entries and their positions.

Formats

• COO (coordinate): three parallel arrays row, col, data. Easiest to build; duplicates are summed on conversion.
• CSR (compressed sparse row): data and indices (column of each value), listed row by row, plus indptr, where row i's entries are data[indptr[i]:indptr[i+1]]. Fast row slicing and matrix-vector products. The workhorse.
• CSC (compressed sparse column): the same idea column by column. Fast column slicing; preferred by some solvers.
• Also LIL and DOK for incremental construction, DIA for banded matrices.

Why machine learning loves sparse matrices

• One-hot features: a "city" column with 5,000 possible values becomes 5,000 columns, each row with a single 1. scikit-learn's OneHotEncoder and CountVectorizer return CSR matrices by default.
• Text: a bag-of-words matrix over a 50,000-word vocabulary has a few hundred non-zeros per document.
• Graphs: an adjacency matrix has a 1 at [i, j] when there is an edge from i to j (see M26). A @ x spreads values along edges, which is one step of PageRank or a graph neural network layer; scipy.sparse.csgraph has shortest paths and connected components.

Remember: store only non-zeros. Build in COO, compute in CSR (rows) or CSC (columns). Sparse mat-vec costs O(nnz), not O(n²).

Practice

Use NumPy (and SciPy where asked). Tests compare with tolerances, and eigenvector/SVD signs are not assumed.

import numpy as np

def solve_system(A, b):
    A = np.asarray(A, dtype=float)
    if A.ndim != 2 or A.shape[0] != A.shape[1]:
        raise ValueError(f"A must be square, got shape {A.shape}")
    if np.linalg.matrix_rank(A) < A.shape[0]:
        raise ValueError("A is singular: no unique solution")
    return np.linalg.solve(A, b)

def residual(A, x, b):
    return float(np.linalg.norm(A @ x - b))

import numpy as np

def fit_poly(x, y, degree):
    x = np.asarray(x, dtype=float)
    X = np.column_stack([x ** p for p in range(degree + 1)])
    coeffs, *_ = np.linalg.lstsq(X, y, rcond=None)
    return coeffs

def predict(coeffs, x):
    x = np.asarray(x, dtype=float)
    return sum(c * x ** p for p, c in enumerate(coeffs))

import numpy as np

def power_iteration(A, iters=500, tol=1e-12):
    A = np.asarray(A, dtype=float)
    v = np.ones(A.shape[0]) / np.sqrt(A.shape[0])
    lam = v @ A @ v
    for _ in range(iters):
        w = A @ v
        v = w / np.linalg.norm(w)
        lam_new = v @ A @ v
        if abs(lam_new - lam) < tol:
            lam = lam_new
            break
        lam = lam_new
    return float(lam), v

import numpy as np

def low_rank(A, k):
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    return (U[:, :k] * s[:k]) @ Vt[:k]

def storage_ratio(shape, k):
    m, n = shape
    return k * (m + n + 1) / (m * n)

def rank_for_energy(A, fraction):
    s = np.linalg.svd(A, compute_uv=False)
    energy = np.cumsum(s ** 2) / (s ** 2).sum()
    return int(np.argmax(energy >= fraction - 1e-12)) + 1

import numpy as np

def pca(X, k):
    X = np.asarray(X, dtype=float)
    mean = X.mean(axis=0)
    Xc = X - mean
    _, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    var = s ** 2 / (X.shape[0] - 1)
    components = Vt[:k]
    return components, var[:k] / var.sum(), Xc @ components.T, mean

import numpy as np

def to_csr(dense):
    dense = np.asarray(dense)
    data, indices, indptr = [], [], [0]
    for row in dense:
        cols = np.nonzero(row)[0]
        data.extend(row[cols])
        indices.extend(cols)
        indptr.append(len(data))
    return np.array(data, dtype=float), np.array(indices, dtype=int), np.array(indptr, dtype=int)

def csr_matvec(data, indices, indptr, x):
    n_rows = len(indptr) - 1
    out = np.zeros(n_rows)
    for i in range(n_rows):
        s, e = indptr[i], indptr[i + 1]
        out[i] = data[s:e] @ x[indices[s:e]]
    return out

Mini-project: a PCA-based compressor

pixi init pcazip && cd pcazip
pixi add python=3.14 numpy matplotlib pytest
pixi run python -m pcazip photo.png --energy 0.95      # writes photo.pcaz.npz and report.png
pixi run pytest -q

Quiz

Recap

• Linear systems are intersecting lines/planes; write them as A x = b and use np.linalg.solve, not inv. Check cond for sensitivity.
• Determinant = area/volume scale factor (sign = flip, zero = squashed). Rank = independent directions.
• L2, L1, and L∞ norms measure size differently: distance, robustness/sparsity, worst case.
• Least squares projects y onto the column space of X; np.linalg.lstsq.
• Eigenvectors are only stretched: A v = λ v. Use eigh for symmetric matrices.
• SVD: rotate, stretch, rotate, for any matrix; truncation gives optimal low-rank approximations; PCA is SVD of centered data.
• Reuse scipy.linalg factorizations (LU, Cholesky). Store mostly-zero data in scipy.sparse CSR/CSC.

Going further

• numpy.linalg reference
• SciPy tutorial: linear algebra
• SciPy tutorial: sparse arrays
• 3Blue1Brown: Essence of Linear Algebra (determinants, eigenvectors)
• Gilbert Strang, Introduction to Linear Algebra, and his free MIT OpenCourseWare 18.06 lectures
• Next: M50 · SciPy for Science & Engineering