NumPy Linear Algebra

Real-World Scenario

A data scientist building a recommendation system uses matrix factorization to decompose a 50,000 user × 10,000 movie rating matrix into latent factors. A quantitative researcher solves a portfolio optimization problem — finding optimal asset weights that minimize variance for a target return. Both reduce to linear algebra. np.linalg provides the tools without requiring a PhD in numerical methods.

Solving Linear Systems

The classic problem: given Ax = b, find x. Appears everywhere from curve fitting to portfolio optimization.

import numpy as np

# A 3×3 system of equations:
# 2x + y - z = 8
# -3x - y + 2z = -11
# -2x + y + 2z = -3
A = np.array([
    [ 2,  1, -1],
    [-3, -1,  2],
    [-2,  1,  2]
], dtype=float)

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

# Solve Ax = b — faster and more stable than computing inv(A) @ b
x = np.linalg.solve(A, b)
print(x)          # [2. 3. -1.]

# Verify: Ax should equal b
residual = np.linalg.norm(A @ x - b)
print(f"Residual: {residual:.2e}")  # ~0.0 — exact solution

Least Squares Regression

When the system is overdetermined (more equations than unknowns — the common case in ML), use least squares to find the best-fit solution.

import numpy as np

rng = np.random.default_rng(42)
n = 200

# Generate data: y = 3x + 2 + noise
x = rng.uniform(0, 10, n)
y = 3 * x + 2 + rng.normal(0, 1, n)

# Design matrix: column of ones for intercept, column of x for slope
X = np.column_stack([np.ones(n), x])   # (200, 2)

# Least squares: minimizes ||Xw - y||^2
# Returns coefficients, residuals, rank, singular values
coeffs, residuals, rank, sv = np.linalg.lstsq(X, y, rcond=None)
intercept, slope = coeffs
print(f"Intercept: {intercept:.4f}")  # ~2.0
print(f"Slope:     {slope:.4f}")      # ~3.0

# The normal equations approach (equivalent but less stable)
# w = (X^T X)^{-1} X^T y
w_normal = np.linalg.solve(X.T @ X, X.T @ y)
print(f"Normal eq: {w_normal}")  # same result

Eigenvalues and Eigenvectors

Eigendecomposition is at the heart of PCA, spectral clustering, and Markov chains.

import numpy as np

# Covariance matrix of a 3-feature dataset
rng = np.random.default_rng(42)
data = rng.standard_normal((100, 3))
data[:, 1] = data[:, 0] * 0.8 + rng.normal(0, 0.3, 100)  # correlated features

cov = np.cov(data.T)   # (3, 3) covariance matrix

# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eigh(cov)  # eigh for symmetric matrices

# Sort by descending eigenvalue (largest = most variance)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

print("Eigenvalues:", eigenvalues.round(3))
print("Variance explained:", (eigenvalues / eigenvalues.sum()).round(3))
# First component likely explains >60% of variance due to correlation

PCA from Scratch using SVD

import numpy as np

rng = np.random.default_rng(42)

# Dataset: 500 samples, 50 features
X = rng.standard_normal((500, 50))
# Introduce correlation structure
X[:, :10] = X[:, :10] @ rng.standard_normal((10, 10))

# Step 1: center the data
X_centered = X - X.mean(axis=0)

# Step 2: SVD decomposition — X = U @ S_diag @ Vt
# U: (n, n), S: (k,) singular values, Vt: (k, p) principal components
U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
# U: (500, 50), S: (50,), Vt: (50, 50)

# Step 3: project onto first 2 principal components
n_components = 2
X_pca = X_centered @ Vt[:n_components].T   # (500, 2)

# Variance explained by each component
variance_ratio = (S ** 2) / (S ** 2).sum()
print("Variance explained by PC1:", f"{variance_ratio[0]:.1%}")
print("Variance explained by PC2:", f"{variance_ratio[1]:.1%}")
print("Projected shape:", X_pca.shape)  # (500, 2)

Matrix Norms and Condition Number

import numpy as np

A = np.array([[1., 2.], [3., 4.]])

# Frobenius norm — square root of sum of squared elements
print(np.linalg.norm(A, 'fro'))  # 5.477

# Spectral norm (largest singular value)
print(np.linalg.norm(A, 2))      # 5.465

# Condition number — ratio of largest to smallest singular value
# High condition number → ill-conditioned → numerical instability in solve()
print(np.linalg.cond(A))         # 14.93 — acceptable

# Rank — number of linearly independent rows/columns
print(np.linalg.matrix_rank(A))  # 2 — full rank

Determinant and Inverse

import numpy as np

A = np.array([[2., 1.], [5., 3.]])

# Determinant — zero means singular (non-invertible)
det = np.linalg.det(A)
print(f"det(A) = {det:.2f}")  # 1.0

# Matrix inverse — use solve() instead for systems of equations
A_inv = np.linalg.inv(A)
print(A_inv)
# [[ 3. -1.]
#  [-5.  2.]]

# Verify: A @ A_inv should equal identity
print(np.allclose(A @ A_inv, np.eye(2)))  # True

QR Decomposition

Useful for numerically stable least squares and Gram-Schmidt orthogonalization.

import numpy as np

# Design matrix for regression
rng = np.random.default_rng(0)
X = rng.standard_normal((100, 5))
y = rng.standard_normal(100)

# QR decomposition: X = Q @ R
# Q: orthonormal columns, R: upper triangular
Q, R = np.linalg.qr(X)

# Solve least squares via QR: w = R^{-1} Q^T y
# More numerically stable than (X^T X)^{-1} X^T y for large systems
w = np.linalg.solve(R, Q.T @ y)
print("Coefficients:", w.round(4))

# Verify against lstsq
w_lstsq, *_ = np.linalg.lstsq(X, y, rcond=None)
print("Match:", np.allclose(w, w_lstsq))  # True