NumPy Linear Algebra

Solving Ax = b

Three equations, three unknowns:

import numpy as np

A = np.array([[ 2.,  1., -1.],
              [-3., -1.,  2.],
              [-2.,  1.,  2.]])
b = np.array([8., -11., -3.])

x = np.linalg.solve(A, b)
print(x)
print(A @ x)
print(np.linalg.norm(A @ x - b))

$ python solve.py
[ 2.  3. -1.]
[  8. -11.  -3.]
1.9860273225978185e-15

The solution is x = 2, y = 3, z = -1. Multiplying back reproduces b as printed, yet the residual norm is about 2e-15, not zero — rounding error in the last bits. Floating-point solutions are checked with a tolerance (np.allclose), never with ==.

solve vs inv on an Ill-Conditioned Matrix

import numpy as np

n = 10
H = 1.0 / (np.arange(n)[:, None] + np.arange(n) + 1)
x_true = np.ones(n)
b = H @ x_true

x_solve = np.linalg.solve(H, b)
x_inv = np.linalg.inv(H) @ b

print(f"cond(H)      {np.linalg.cond(H):.2e}")
print(f"solve error  {np.abs(x_solve - x_true).max():.2e}")
print(f"inv error    {np.abs(x_inv - x_true).max():.2e}")

$ python solve_vs_inv.py
cond(H)      1.60e+13
solve error  1.69e-04
inv error    5.55e-03

H is the 10 × 10 Hilbert matrix, a textbook example of a nearly singular matrix, and its condition number is around 1e13. Both methods lose digits — the true answer is all ones — but the explicit inverse’s worst error is about 30 times larger than solve’s here. That gap is why solve is the default.

A Singular Matrix

import numpy as np

A = np.array([[1., 2.], [2., 4.]])
print(np.linalg.det(A), np.linalg.matrix_rank(A))
try:
    np.linalg.solve(A, [3., 6.])
except np.linalg.LinAlgError as e:
    print("LinAlgError:", e)

x, *_ = np.linalg.lstsq(A, [3., 6.], rcond=None)
print(x)

$ python singular.py
0.0 1
LinAlgError: Singular matrix
[0.6 1.2]

The second row is twice the first, so the rank is 1 and the determinant is zero. solve refuses. lstsq still answers: of the infinitely many solutions to x + 2y = 3, it returns the one with the smallest length.

Least Squares Regression

With more equations than unknowns, lstsq finds the best fit.

import numpy as np

rng = np.random.default_rng(42)
x = rng.uniform(0, 10, 200)
y = 3 * x + 2 + rng.normal(0, 1, 200)

X = np.column_stack([np.ones_like(x), x])
coef, resid, rank, sv = np.linalg.lstsq(X, y, rcond=None)

print(X.shape)
print("intercept, slope:", coef.round(4))
print("rank:", rank)

$ python lstsq.py
(200, 2)
intercept, slope: [1.6908 3.0659]
rank: 2

The column of ones gives the model an intercept. The slope lands very close to the true 3; the intercept (about 1.69 against a true 2) is further off, because with x spread over 0–10 the data pins down the slope much more tightly than the value at x = 0. A rank of 2 means both columns carry independent information.

Eigenvalues of a Covariance Matrix

import numpy as np

rng = np.random.default_rng(42)
data = rng.standard_normal((500, 3))
data[:, 1] = 0.8 * data[:, 0] + 0.3 * rng.standard_normal(500)

cov = np.cov(data, rowvar=False)
vals, vecs = np.linalg.eigh(cov)
order = np.argsort(vals)[::-1]
vals, vecs = vals[order], vecs[:, order]

print("eigenvalues:", vals.round(3))
print("share:", (vals / vals.sum()).round(3))
print("top direction:", vecs[:, 0].round(3))

$ python eigen.py
eigenvalues: [1.796 0.935 0.055]
share: [0.645 0.336 0.02 ]
top direction: [-0.775 -0.632 -0.007]

eigh is for symmetric matrices like covariances and returns eigenvalues in ascending order, so they’re re-sorted. Columns 0 and 1 are strongly correlated, so the largest eigenvalue carries most of the variance and its eigenvector weights those two columns and barely touches column 2. The sign of an eigenvector is arbitrary — -v is equally valid.

PCA with SVD

import numpy as np

rng = np.random.default_rng(42)
latent = rng.standard_normal((500, 2))
X = latent @ rng.standard_normal((2, 10)) + 0.1 * rng.standard_normal((500, 10))

Xc = X - X.mean(axis=0)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
ratio = S**2 / (S**2).sum()

print(U.shape, S.shape, Vt.shape)
print("variance ratio:", ratio[:4].round(4))
print("projected:", (Xc @ Vt[:2].T).shape)

$ python pca_svd.py
(500, 10) (10,) (10, 10)
variance ratio: [0.8086 0.1841 0.0011 0.0011]
projected: (500, 2)

The 10 columns were built from only 2 hidden factors plus a little noise, and the SVD finds that: the first two components hold almost all the variance and the rest is noise. Projecting onto the top two rows of Vt compresses each sample from 10 numbers to 2.

Norms, Rank and Determinant

import numpy as np

A = np.array([[2., 1.], [5., 3.]])
print(np.linalg.norm(A), np.linalg.norm(A, 2))
print(np.linalg.matrix_rank(A), np.linalg.cond(A))
print(np.linalg.det(A))
print(np.linalg.inv(A))

$ python norms.py
6.244997998398398 6.242943383865535
2 38.97434209415041
1.0000000000000009
[[ 3. -1.]
 [-5.  2.]]

The first norm is the Frobenius norm (square root of the sum of squares); norm(A, 2) is the largest singular value. The determinant is mathematically exactly 1 (2·3 − 1·5), but det works through an LU factorization in floating point, so it can print as a value a hair off 1. Never test a determinant with == 0.

Practice

1. Solve x + y = 10, x − y = 2.
2. Build the 12 × 12 Hilbert matrix and print its condition number. Roughly how many of float64’s ~16 significant digits would you expect to lose?
3. Fit y = a·x² + b with lstsq to x = [0, 1, 2, 3], y = [1, 3, 9, 19].

import numpy as np

print(np.linalg.solve([[1., 1.], [1., -1.]], [10., 2.]))

$ python practice_solve2.py
[6. 4.]

import numpy as np

n = 12
H = 1.0 / (np.arange(n)[:, None] + np.arange(n) + 1)
c = np.linalg.cond(H)
print(f"{c:.2e}  ~{np.log10(c):.0f} digits lost")

$ python practice_hilbert12.py
1.81e+16  ~16 digits lost

import numpy as np

x = np.array([0., 1., 2., 3.])
y = np.array([1., 3., 9., 19.])
X = np.column_stack([x**2, np.ones_like(x)])
coef, *_ = np.linalg.lstsq(X, y, rcond=None)
print(coef.round(6))

$ python practice_quadfit.py
[2. 1.]