Appendix D — numpy / scipy.linalg Cheatsheet

This is the dense lookup card: every numpy or scipy operation the book uses, the mathematics it computes, and a one-line example you can paste and run. The book's philosophy is that code demonstrates the math (it never replaces the proof or the hand computation), and this card is the dictionary that lets you translate a formula into a call. Each example assumes the imports import numpy as np and from scipy import linalg as sla, and that A, B are square matrices and x, b are vectors, defined as needed.

A reminder that governs the whole card: math indexes from 1, numpy from 0 (so $a_{ij}$ is A[i-1, j-1]), @ is matrix multiplication while * is elementwise, and never test equality with == — use np.allclose. With that, the dictionary.

D.1 Creating and inspecting arrays

Common Pitfall — A 1-D array of shape (n,) is neither a row nor a column; it is a flat vector that numpy reshapes contextually. A @ x treats x as a column and returns a 1-D result. If you genuinely need a column matrix of shape (n,1), write x.reshape(-1, 1) or x[:, None]. Most of the book is happier with flat 1-D vectors.

D.2 The core operations

A = np.array([[2.,1.],[1.,3.]]); x = np.array([1.,2.])
A @ x            # -> array([4., 7.])     matrix-vector product A x
A.T @ A          # -> the symmetric matrix  A^T A  (appears in least squares)

D.3 Solving systems and inverting (Ch. 4, 9)

A = np.array([[2.,1.],[1.,3.]]); b = np.array([1.,2.])
np.linalg.solve(A, b)          # -> array([0.2, 0.6])   solves A x = b
np.linalg.det(A)               # -> 5.0
np.linalg.matrix_power(np.array([[1,1],[0,1]]), 3)   # -> [[1,3],[0,1]]  the shear cubed

Computational Note — Prefer np.linalg.solve(A, b) to np.linalg.inv(A) @ b. They are mathematically equal, but solve is faster and numerically more accurate — forming the inverse explicitly wastes work and amplifies rounding error. Compute the inverse only when you genuinely need the matrix $A^{-1}$ itself, not merely a solution.

D.4 Eigenvalues: eig versus eigh (Ch. 23–27)

S = np.array([[2.,1.],[1.,2.]])      # symmetric
w, V = np.linalg.eigh(S)             # w -> [1., 3.]  (REAL, ASCENDING); V columns orthonormal
R = np.array([[0.,-1.],[1.,0.]])     # a 90-degree rotation
np.linalg.eigvals(R)                 # -> [0.+1.j, 0.-1.j]  COMPLEX pair: no real fixed direction

D.5 The factorizations (Ch. 10, 20, 28, 30)

A = np.array([[3.,0.],[0.,-2.],[0.,0.]])
U, s, Vt = np.linalg.svd(A)          # s -> [3., 2.]  singular values, DESCENDING
# reconstruct: build a full m-by-n Sigma, then U @ Sigma @ Vt == A
Sigma = np.zeros_like(A); Sigma[:len(s),:len(s)] = np.diag(s)
np.allclose(U @ Sigma @ Vt, A)       # -> True   (note Vt is V^T already, not V)

Q, R = np.linalg.qr(np.array([[1.,1.],[1.,0.],[0.,1.]]))
np.allclose(Q.T @ Q, np.eye(2))      # -> True   Q has orthonormal columns (3x2 reduced by default)

P, L, U = sla.lu(np.array([[0.,2.],[1.,1.]]))
np.allclose(P @ L @ U, [[0.,2.],[1.,1.]])        # -> True   convention is A = P @ L @ U

D.6 Norms, conditioning, and least squares (Ch. 18, 19, 38)

np.linalg.norm([3., 4.])                 # -> 5.0
np.linalg.cond(np.array([[1.,0.],[0.,1e-3]]))    # -> 1000.0   (ill-conditioned)

# Least squares: fit a line y = c0 + c1 t through 3 points
M = np.array([[1.,0.],[1.,1.],[1.,2.]]); y = np.array([1.,2.,2.])
coef, *_ = np.linalg.lstsq(M, y, rcond=None)     # coef -> [1.1667, 0.5]

Computational Note — np.linalg.norm(A) on a matrix defaults to the Frobenius norm (treat all entries as one long vector), whereas on a vector it defaults to the Euclidean 2-norm. To get the matrix's operator 2-norm — its largest singular value, the gain of the worst-case input direction (Ch. 30) — you must pass the explicit order, np.linalg.norm(A, 2). The two matrix norms are different numbers; name the one you mean. And lstsq returns a tuple (solution, residuals, rank, singular values), so unpack it — coef, *_ = np.linalg.lstsq(...) — and pass rcond=None to silence a version-dependent warning and get the documented behavior.

D.7 A one-screen summary

import numpy as np
from scipy import linalg as sla

A @ B                      # matrix product (composition)     |  A * B  is elementwise!
A.T,  A.conj().T           # transpose A^T  ;  conjugate transpose A^*
np.linalg.solve(A, b)      # solve A x = b   (preferred over inv(A) @ b)
np.linalg.inv(A)           # inverse A^{-1}
np.linalg.det(A)           # determinant
np.linalg.matrix_rank(A)   # rank
np.linalg.matrix_power(A,k)# A^k
np.linalg.eig(A)           # general eigen-decomposition (unsorted, maybe complex)
np.linalg.eigh(S)          # symmetric/Hermitian: REAL, ascending, orthonormal
np.linalg.svd(A)           # U, s, Vt   (s descending; Vt is V^T)
np.linalg.qr(A)            # Q, R       (reduced by default; signs vary)
np.linalg.lstsq(A,b,rcond=None)   # least squares  min ||A x - b||
np.linalg.norm(x)          # ||x||_2 ;  norm(A) is Frobenius ;  norm(A,2) is sigma_1
np.linalg.cond(A)          # condition number sigma_1/sigma_n
sla.lu(A)                  # P, L, U   with A = P @ L @ U
sla.expm(A)                # matrix exponential e^A
sla.null_space(A)          # orthonormal basis of N(A)
sla.cholesky(A, lower=True)# lower L with A = L L^T   (scipy default is UPPER)

Pin this to the wall. The math comes first in this book; numpy is how you watch it come true. When a result surprises you, the cause is almost always a convention on this card — a default order, a transpose already taken, a sign you cannot pin down, or * where you meant @.