Linear Algebra Operations: NumPy linalg

NumPy's linear algebra module (numpy.linalg) provides high-performance implementations of essential matrix operations: multiplication, inversion, decompositions (SVD, QR, Cholesky), eigenvalue solving, and linear equation solving. These operations are crucial for machine learning, scientific simulation, signal processing, and numerical optimization. Understanding which operations to use and their computational complexity helps you write fast algorithms for data-heavy problems.

Matrix Operations: Multiplication and Transposition

Matrix multiplication is the workhorse of linear algebra. NumPy delegates to BLAS (Basic Linear Algebra Subprograms), which is optimized for your CPU.

import numpy as np
import timeit

# Create matrices
A = np.random.randn(1000, 500)   # 1000x500
B = np.random.randn(500, 1000)   # 500x1000

# Matrix multiplication (1000x500) @ (500x1000) → 1000x1000
C = A @ B  # Modern Python syntax; equivalent to np.dot(A, B)
print(f"Result shape: {C.shape}")  # (1000, 1000)

# Alternative: np.dot() or np.matmul()
C_dot = np.dot(A, B)
C_matmul = np.matmul(A, B)
assert np.allclose(C, C_dot)
assert np.allclose(C, C_matmul)

# Batch matrix multiplication (NumPy 1.16+)
A_batch = np.random.randn(32, 10, 20)   # 32 batches of 10x20 matrices
B_batch = np.random.randn(32, 20, 15)   # 32 batches of 20x15 matrices
C_batch = A_batch @ B_batch             # (32, 10, 15)
print(f"Batch result shape: {C_batch.shape}")

# Transpose: efficient view (no copy)
A_T = A.T  # equivalent to np.transpose(A)
print(f"Transpose shape: {A_T.shape}")

# Hermitian transpose (conjugate transpose, for complex matrices)
C_complex = np.random.randn(5, 5) + 1j * np.random.randn(5, 5)
C_hermitian = C_complex.conj().T  # or np.conj(np.transpose(C_complex))

BLAS implementations are highly optimized; matrix multiplication on large matrices executes near peak CPU bandwidth. Prefer @ or np.dot() over explicit element-wise operations.

Solving Linear Systems

The most common linear algebra problem: solve Ax = b for unknown x. Use np.linalg.solve() instead of manually inverting A.

import numpy as np

# Problem: solve Ax = b
A = np.array([[3, 1],
              [1, 2]], dtype=float)
b = np.array([9, 8], dtype=float)

# Solution 1: np.linalg.solve() — preferred, fast and numerically stable
x = np.linalg.solve(A, b)
print(f"Solution: x = {x}")  # [2, 3]

# Verify: Ax should equal b
result = A @ x
assert np.allclose(result, b)

# Solution 2: Manual inversion (slow and unstable for ill-conditioned matrices)
A_inv = np.linalg.inv(A)
x_inv = A_inv @ b
assert np.allclose(x, x_inv)  # Same result, but slower and less stable

# Overdetermined system (more equations than unknowns): least-squares
A_overdetermined = np.random.randn(100, 5)  # 100 equations, 5 unknowns
b_overdetermined = np.random.randn(100)
x_lstsq, residuals, rank, s = np.linalg.lstsq(A_overdetermined, b_overdetermined, rcond=None)
print(f"Least-squares solution shape: {x_lstsq.shape}")  # (5,)
print(f"Residuals (sum of squared errors): {residuals}")

Never use np.linalg.inv(A) @ b to solve Ax = b—it's slower and numerically unstable. Always use np.linalg.solve().

Matrix Decompositions

Decompositions express a matrix as a product of simpler matrices, enabling fast algorithms for solving, inversion, and rank computation.

Singular Value Decomposition (SVD)

SVD factors A = U @ Σ @ V.T, where U and V are orthogonal, Σ is diagonal. Essential for low-rank approximation, dimensionality reduction, and computing rank.

import numpy as np

# Data matrix (e.g., images or samples)
A = np.random.randn(100, 50)  # 100 samples, 50 features

# Full SVD
U, S, Vt = np.linalg.svd(A, full_matrices=True)
print(f"U shape: {U.shape}")  # (100, 100)
print(f"S shape: {S.shape}")  # (50,) — singular values
print(f"Vt shape: {Vt.shape}")  # (50, 50)

# Reconstruct A
A_reconstructed = U @ np.diag(S[:min(len(U), len(Vt))]) @ Vt[:len(S)]
# (simplified; exact reconstruction requires careful indexing)

# Low-rank approximation (keep top k singular values)
k = 10
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
A_approx = U_k @ np.diag(S_k) @ Vt_k  # (100, 50) with rank k
print(f"Low-rank approximation shape: {A_approx.shape}")

# Compute matrix rank (number of non-zero singular values)
rank = np.sum(S > 1e-10)
print(f"Matrix rank: {rank}")

# Condition number (ratio of largest to smallest singular value)
condition_number = S[0] / S[-1]
print(f"Condition number: {condition_number:.2e}")  # Ill-conditioned if large

QR Decomposition

QR factors A = Q @ R, where Q is orthogonal, R is upper triangular. Useful for solving least-squares and orthogonalizing columns.

import numpy as np

A = np.random.randn(10, 5)
Q, R = np.linalg.qr(A)

print(f"Q shape: {Q.shape}")  # (10, 5) or (10, 10) depending on mode
print(f"R shape: {R.shape}")  # (5, 5)

# Q is orthogonal (Q.T @ Q ≈ I)
assert np.allclose(Q.T @ Q, np.eye(Q.shape[1]))

# Reconstruct A
A_reconstructed = Q @ R
assert np.allclose(A, A_reconstructed)

Cholesky Decomposition

For symmetric positive-definite matrices, Cholesky factors A = L @ L.T where L is lower triangular. Faster than LU for this special case.

import numpy as np

# Symmetric positive-definite matrix (e.g., covariance)
n = 5
A = np.random.randn(n, n)
A = A @ A.T  # A is now symmetric positive-definite
print(f"Eigenvalues (all positive): {np.linalg.eigvalsh(A)}")

# Cholesky decomposition
L = np.linalg.cholesky(A)
print(f"L shape: {L.shape}")  # (5, 5) lower triangular

# Reconstruct A
A_reconstructed = L @ L.T
assert np.allclose(A, A_reconstructed)

# Solve Ax = b using Cholesky (faster than general solve for SPD)
b = np.random.randn(n)
# Method 1: General solver
x1 = np.linalg.solve(A, b)
# Method 2: Using Cholesky (illustrative; np.linalg.solve uses Cholesky internally)
y = np.linalg.solve(L, b)  # solve Ly = b
x2 = np.linalg.solve(L.T, y)  # solve L.T x = y
assert np.allclose(x1, x2)

Eigenvalues and Eigenvectors

Eigendecomposition finds vectors v and scalars λ such that Av = λv. Essential for understanding matrix behavior, stability analysis, and principal component analysis.

import numpy as np

# Sample: symmetric covariance matrix
A = np.array([[4, 1],
              [1, 3]], dtype=float)

# Eigendecomposition
eigenvalues, eigenvectors = np.linalg.eig(A)
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# Verify: A @ v = λ @ v
v0 = eigenvectors[:, 0]
lam0 = eigenvalues[0]
assert np.allclose(A @ v0, lam0 * v0)

# For symmetric matrices, use eigvalsh (more stable, returns sorted values)
eigenvalues_sym, eigenvectors_sym = np.linalg.eigh(A)
print(f"Symmetric eigenvalues (sorted): {eigenvalues_sym}")

# Principal component analysis (PCA) using SVD
X = np.random.randn(100, 50)  # 100 samples, 50 features
U, S, Vt = np.linalg.svd(X, full_matrices=False)
# Vt rows are principal components (eigenvectors of X.T @ X)
# S**2 / (n-1) are the variances (eigenvalues)

Matrix Norms and Condition Numbers

Matrix norms measure size; condition numbers measure sensitivity to perturbations (numerical stability).

import numpy as np

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

# Frobenius norm (Euclidean norm of all elements)
frobenius = np.linalg.norm(A, 'fro')
print(f"Frobenius norm: {frobenius:.4f}")

# Spectral norm (largest singular value)
spectral = np.linalg.norm(A, 2)  # or np.linalg.norm(A) — default is spectral
print(f"Spectral norm: {spectral:.4f}")

# Nuclear norm (sum of singular values)
nuclear = np.linalg.norm(A, 'nuc')
print(f"Nuclear norm: {nuclear:.4f}")

# Condition number (spectral norm)
condition_number = np.linalg.cond(A)
print(f"Condition number: {condition_number:.2e}")
# Small (close to 1): well-conditioned, stable solve
# Large (> 1e10): ill-conditioned, solve is numerically unstable

Key Takeaways

• Use np.linalg.solve() for Ax = b; never invert and multiply—it's slower and unstable.
• SVD reveals rank, enables low-rank approximation, and is the basis for PCA.
• QR and Cholesky decompositions speed up specific problem classes.
• Condition number predicts numerical stability: large values indicate ill-conditioned systems.
• BLAS backend makes NumPy matrix operations competitive with compiled languages.

Frequently Asked Questions

What is the time complexity of matrix multiplication?

(m x n) @ (n x p) is O(mnp). For square matrices (m = n = p), it's O(n^3). Modern BLAS uses algorithms like Strassen for better constants but same asymptotic complexity.

Should I use @ or np.dot() for matrix multiplication?

Both are equivalent for 2D arrays. @ is more explicit and modern (Python 3.5+); use it. np.dot() also handles vector dot products and higher-dimensional broadcasting, so choose based on clarity.

When should I use lstsq instead of solve?

Use lstsq for overdetermined systems (more equations than unknowns) or when A is rank-deficient. Use solve for square, full-rank systems.

Is SVD always the best for low-rank approximation?

SVD is optimal but O(n^3). For interactive systems, consider randomized SVD or truncated methods. For streaming data, incremental methods are better.

How do I know if my matrix is ill-conditioned?

Compute condition number: np.linalg.cond(A). Values > 1e10 indicate problems. Also, eigenvalues spread across many orders of magnitude suggest ill-conditioning.

Further Reading

• NumPy Linear Algebra Reference
• Matrix Decompositions Explained
• Eigenvalue Solvers and Applications
• Understanding Condition Numbers