math-linear-algebra-computation

Solving systems of linear equations (Ax = b)

Linear Algebra Computation

Metadata

Overview

Computational linear algebra using NumPy and SciPy for matrix operations, decompositions, solving linear systems, and dimensionality reduction. Essential for scientific computing, machine learning, and numerical analysis.

When to Use This Skill

When NOT to Use This Skill

Prerequisites

Core Concepts

1. Matrix Operations Fundamentals

Efficient Matrix Multiplication

import numpy as np

# Basic matrix multiplication
A = np.random.rand(1000, 500)
B = np.random.rand(500, 300)

# Standard: A @ B or np.dot(A, B)
C = A @ B  # Preferred Python 3.5+

# Matrix-vector: use @ for clarity
x = np.random.rand(500)
y = A @ x

# Batch matrix multiplication
batch_A = np.random.rand(10, 100, 50)
batch_B = np.random.rand(10, 50, 30)
batch_C = np.einsum('ijk,ikl->ijl', batch_A, batch_B)
# Or: batch_C = batch_A @ batch_B (broadcasting)

# Outer product
u = np.random.rand(100)
v = np.random.rand(200)
outer = np.outer(u, v)  # 100x200 matrix

# Inner product (dot product)
inner = np.inner(u, u)  # Scalar

Matrix Norms and Properties

from numpy.linalg import norm

A = np.random.rand(100, 100)

# Frobenius norm (default)
fro_norm = norm(A)  # sqrt(sum of squared elements)

# 2-norm (spectral norm, largest singular value)
spectral_norm = norm(A, ord=2)

# Nuclear norm (sum of singular values)
nuclear_norm = norm(A, ord='nuc')

# Infinity norm (max row sum)
inf_norm = norm(A, ord=np.inf)

# 1-norm (max column sum)
one_norm = norm(A, ord=1)

# Condition number (sensitivity to perturbations)
cond = np.linalg.cond(A)  # ||A|| * ||A^-1||
# High condition number = ill-conditioned matrix

2. Solving Linear Systems

Direct Methods (Ax = b)

from scipy.linalg import solve, lu_factor, lu_solve, cho_factor, cho_solve

# General system: prefer solve() over inv()
A = np.random.rand(100, 100)
b = np.random.rand(100)

# NEVER do this (slow and numerically unstable):
# x = np.linalg.inv(A) @ b  # BAD!

# ALWAYS do this instead:
x = np.linalg.solve(A, b)  # GOOD!

# Multiple right-hand sides
B = np.random.rand(100, 5)
X = np.linalg.solve(A, B)

# LU decomposition for repeated solves
lu, piv = lu_factor(A)
x1 = lu_solve((lu, piv), b)
x2 = lu_solve((lu, piv), b + 0.1)  # Reuse factorization

# Symmetric positive definite: use Cholesky
A_spd = A @ A.T + np.eye(100)  # Force SPD
c, low = cho_factor(A_spd)
x_cho = cho_solve((c, low), b)

# Check solution quality
residual = norm(A @ x - b) / norm(b)
print(f"Relative residual: {residual:.2e}")

Least Squares Problems

from numpy.linalg import lstsq
from scipy.linalg import qr

# Overdetermined system (more equations than unknowns)
A = np.random.rand(200, 100)  # 200 equations, 100 unknowns
b = np.random.rand(200)

# Solve min ||Ax - b||^2
x, residuals, rank, s = lstsq(A, b, rcond=None)

# Underdetermined system (infinite solutions, find minimum norm)
A_under = np.random.rand(50, 100)  # 50 equations, 100 unknowns
b_under = np.random.rand(50)
x_min_norm = lstsq(A_under, b_under, rcond=None)[0]

# Weighted least squares
W = np.diag(np.random.rand(200))  # Weight matrix
A_weighted = W @ A
b_weighted = W @ b
x_weighted = lstsq(A_weighted, b_weighted, rcond=None)[0]

# QR-based least squares (more stable for ill-conditioned A)
Q, R = qr(A, mode='economic')
x_qr = solve(R, Q.T @ b)

3. Matrix Decompositions

Singular Value Decomposition (SVD)

from numpy.linalg import svd
from scipy.sparse.linalg import svds

# Full SVD: A = U @ S @ V^T
A = np.random.rand(100, 50)
U, s, Vt = svd(A, full_matrices=False)

# U: 100x50 (left singular vectors)
# s: 50 (singular values, sorted descending)
# Vt: 50x50 (right singular vectors transposed)

# Reconstruct original matrix
S = np.diag(s)
A_reconstructed = U @ S @ Vt
assert np.allclose(A, A_reconstructed)

# Low-rank approximation (truncated SVD)
k = 10  # Keep top 10 components
A_approx = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

# Approximation error
error = norm(A - A_approx, 'fro') / norm(A, 'fro')
print(f"Relative error with k={k}: {error:.4f}")

# Sparse SVD (for large matrices, top k singular values only)
A_large = np.random.rand(1000, 500)
k = 20
U_k, s_k, Vt_k = svds(A_large, k=k)
# Note: svds returns singular values in ASCENDING order
# Reverse them for consistency
idx = np.argsort(s_k)[::-1]
s_k = s_k[idx]
U_k = U_k[:, idx]
Vt_k = Vt_k[idx, :]

# Matrix pseudoinverse via SVD
tolerance = 1e-10
s_inv = np.array([1/si if si > tolerance else 0 for si in s])
A_pinv = Vt.T @ np.diag(s_inv) @ U.T

Eigenvalue Decomposition

from numpy.linalg import eig, eigh
from scipy.linalg import eig as scipy_eig

# General matrix (possibly complex eigenvalues)
A = np.random.rand(100, 100)
eigenvalues, eigenvectors = eig(A)

# A @ v = λ * v
# eigenvectors[:, i] corresponds to eigenvalues[i]

# Symmetric/Hermitian matrices (real eigenvalues, orthogonal eigenvectors)
A_sym = (A + A.T) / 2
eigenvalues_sym, eigenvectors_sym = eigh(A_sym)
# eigenvalues are sorted in ascending order

# Verify decomposition
Lambda = np.diag(eigenvalues_sym)
assert np.allclose(A_sym, eigenvectors_sym @ Lambda @ eigenvectors_sym.T)

# Generalized eigenvalue problem: A @ v = λ * B @ v
B = np.random.rand(100, 100)
B_spd = B @ B.T + np.eye(100)
eigenvalues_gen, eigenvectors_gen = scipy_eig(A_sym, B_spd)

# Power iteration for dominant eigenvalue
def power_iteration(A, num_iterations=100):
    v = np.random.rand(A.shape[1])
    for _ in range(num_iterations):
        v = A @ v
        v = v / norm(v)
    eigenvalue = (v @ A @ v) / (v @ v)
    return eigenvalue, v

lambda_max, v_max = power_iteration(A_sym)

QR Decomposition

from numpy.linalg import qr

# A = Q @ R (Q orthogonal, R upper triangular)
A = np.random.rand(100, 50)

# Full QR (Q is 100x100)
Q_full, R_full = qr(A, mode='complete')

# Economic QR (Q is 100x50, preferred)
Q, R = qr(A, mode='reduced')

# Verify orthogonality
assert np.allclose(Q.T @ Q, np.eye(50))
assert np.allclose(A, Q @ R)

# Applications:
# 1. Least squares (already shown)
# 2. Gram-Schmidt orthogonalization
def gram_schmidt(A):
    Q, R = qr(A)
    return Q

# 3. QR algorithm for eigenvalues
def qr_eigenvalue_algorithm(A, iterations=100):
    A_k = A.copy()
    for _ in range(iterations):
        Q, R = qr(A_k)
        A_k = R @ Q
    # Diagonal of A_k approximates eigenvalues
    return np.diag(A_k)

Cholesky Decomposition

from numpy.linalg import cholesky
from scipy.linalg import cholesky as scipy_cholesky

# A = L @ L^T (A must be symmetric positive definite)
A_spd = np.random.rand(100, 50)
A_spd = A_spd @ A_spd.T + np.eye(100)  # Force SPD

# Lower triangular
L = cholesky(A_spd)
assert np.allclose(A_spd, L @ L.T)

# Upper triangular (SciPy)
U = scipy_cholesky(A_spd, lower=False)
assert np.allclose(A_spd, U.T @ U)

# Solving Ax = b with Cholesky
# A = L @ L^T => L @ (L^T @ x) = b
# Solve L @ y = b, then L^T @ x = y
from scipy.linalg import solve_triangular

b = np.random.rand(100)
y = solve_triangular(L, b, lower=True)
x = solve_triangular(L.T, y, lower=False)
assert np.allclose(A_spd @ x, b)

# Updating Cholesky factorization (rank-1 update)
def cholesky_rank1_update(L, u):
    """Update L when A' = A + u @ u^T"""
    from scipy.linalg.lapack import dchud
    L_updated, info = dchud(L.T, u)
    return L_updated.T

4. Dimensionality Reduction

Principal Component Analysis (PCA)

from sklearn.decomposition import PCA

# Data matrix: n samples x p features
X = np.random.rand(1000, 100)

# Center the data
X_centered = X - X.mean(axis=0)

# Method 1: SVD-based PCA (preferred for stability)
U, s, Vt = svd(X_centered, full_matrices=False)

# Principal components (eigenvectors of covariance)
components = Vt  # p x p

# Explained variance
explained_variance = (s ** 2) / (X.shape[0] - 1)
explained_variance_ratio = explained_variance / explained_variance.sum()

# Project data onto top k components
k = 10
X_reduced = U[:, :k] @ np.diag(s[:k])

# Reconstruct approximation
X_reconstructed = X_reduced @ Vt[:k, :] + X.mean(axis=0)

# Method 2: Eigendecomposition of covariance matrix
cov_matrix = np.cov(X_centered.T)  # p x p
eigenvalues, eigenvectors = eigh(cov_matrix)

# Sort in descending order
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# Same principal components (up to sign)
assert np.allclose(np.abs(components[0]), np.abs(eigenvectors[:, 0]))

# Scikit-learn PCA (recommended for production)
pca = PCA(n_components=10)
X_pca = pca.fit_transform(X)
print(f"Explained variance ratio: {pca.explained_variance_ratio_.sum():.4f}")

Randomized SVD (Fast Approximation)

from sklearn.utils.extmath import randomized_svd

# For very large matrices, randomized SVD is much faster
A_large = np.random.rand(10000, 5000)
k = 50

# Randomized SVD
U_rand, s_rand, Vt_rand = randomized_svd(A_large, n_components=k, random_state=42)

# Compare with truncated SVD
from scipy.sparse.linalg import svds
U_exact, s_exact, Vt_exact = svds(A_large, k=k)

# Randomized is approximate but much faster
# Error typically < 1% for k << min(m, n)

5. Sparse Matrices

Sparse Matrix Operations

from scipy.sparse import csr_matrix, csc_matrix, lil_matrix
from scipy.sparse.linalg import spsolve, gmres, cg

# Create sparse matrix (row-based for fast row operations)
row = np.array([0, 0, 1, 2, 2, 2])
col = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6])
A_sparse = csr_matrix((data, (row, col)), shape=(3, 3))

# Dense equivalent
A_dense = A_sparse.toarray()

# Matrix-vector multiplication
x = np.array([1, 2, 3])
y = A_sparse @ x  # Fast for sparse matrices

# Solving sparse linear systems
n = 1000
# Create sparse positive definite matrix
from scipy.sparse import diags
A_sparse_spd = diags([1, -2, 1], [-1, 0, 1], shape=(n, n)).tocsr()
A_sparse_spd = A_sparse_spd @ A_sparse_spd.T + diags([10], [0], shape=(n, n))

b = np.random.rand(n)

# Direct solver (LU-based)
x_direct = spsolve(A_sparse_spd, b)

# Iterative solvers (better for very large systems)
# Conjugate gradient (for SPD matrices)
x_cg, info = cg(A_sparse_spd, b, tol=1e-8)

# GMRES (for general matrices)
x_gmres, info = gmres(A_sparse_spd, b, tol=1e-8)

# Check convergence
if info == 0:
    print("Converged successfully")
elif info > 0:
    print(f"Did not converge after {info} iterations")

Patterns and Best Practices

Pattern 1: Never Invert Matrices

# WRONG: Explicitly computing inverse
A = np.random.rand(100, 100)
b = np.random.rand(100)
x = np.linalg.inv(A) @ b  # Slow and unstable!

# RIGHT: Use solve
x = np.linalg.solve(A, b)  # Fast and stable!

# WRONG: Computing (A^T A)^-1 A^T b (normal equations)
A = np.random.rand(200, 100)
b = np.random.rand(200)
x = np.linalg.inv(A.T @ A) @ A.T @ b  # Numerically unstable!

# RIGHT: Use lstsq
x = np.linalg.lstsq(A, b, rcond=None)[0]  # Stable QR-based!

Pattern 2: Check Conditioning

def solve_with_conditioning_check(A, b, threshold=1e12):
    cond = np.linalg.cond(A)
    if cond > threshold:
        print(f"Warning: Matrix is ill-conditioned (cond={cond:.2e})")
        print("Consider regularization or preconditioning")

    x = np.linalg.solve(A, b)

    # Verify solution
    residual = norm(A @ x - b) / norm(b)
    if residual > 1e-6:
        print(f"Warning: Large residual {residual:.2e}")

    return x

Pattern 3: Memory-Efficient Matrix Operations

# Avoid creating intermediate large matrices
n = 10000

# WRONG: Creates intermediate n x n matrix
A = np.random.rand(n, n)
B = np.random.rand(n, n)
result = (A @ B) @ A.T  # Two n x n intermediates!

# BETTER: Use einsum or explicit ordering
result = A @ (B @ A.T)  # Only one intermediate

# BEST: For complex expressions, use einsum
# result = np.einsum('ij,jk,ki->i', A, B, A)  # No intermediates!

Pattern 4: Numerical Stability

# Computing variance: two-pass algorithm
X = np.random.rand(1000, 100)

# WRONG: Naive variance (numerically unstable)
mean = X.mean(axis=0)
variance_naive = ((X - mean) ** 2).mean(axis=0)

# RIGHT: Use built-in (Welford's algorithm)
variance_stable = np.var(X, axis=0)

# For covariance: use np.cov, not manual computation
cov_stable = np.cov(X.T)

Quick Reference

Function Selection Guide

# Solving Ax = b:
# - General: np.linalg.solve(A, b)
# - SPD: scipy.linalg.cho_solve(cho_factor(A), b)
# - Sparse: scipy.sparse.linalg.spsolve(A, b)
# - Overdetermined: np.linalg.lstsq(A, b)

# Matrix decomposition:
# - SVD: np.linalg.svd(A) or scipy.sparse.linalg.svds(A, k)
# - Eigenvalues: np.linalg.eigh(A) for symmetric
# - QR: np.linalg.qr(A)
# - Cholesky: np.linalg.cholesky(A) for SPD

# Matrix properties:
# - Norm: np.linalg.norm(A, ord=...)
# - Condition: np.linalg.cond(A)
# - Rank: np.linalg.matrix_rank(A)
# - Determinant: np.linalg.det(A)

Performance Tips

  1. Use @ instead of np.dot() for clarity and speed
  2. Prefer solve() over inv() (10x faster, more accurate)
  3. Use eigh() for symmetric matrices (2x faster than eig())
  4. Set full_matrices=False in SVD for efficiency
  5. Use sparse matrices when >90% zeros
  6. Vectorize operations, avoid Python loops
  7. Use einsum() for complex tensor operations

Anti-Patterns

Anti-Pattern 1: Matrix Inversion

# NEVER do this
x = np.linalg.inv(A) @ b

# ALWAYS do this
x = np.linalg.solve(A, b)

Anti-Pattern 2: Ignoring Conditioning

# WRONG: Blindly solving without checking
x = np.linalg.solve(A, b)

# RIGHT: Check condition number first
cond = np.linalg.cond(A)
if cond > 1e12:
    # Add regularization
    A_reg = A + 1e-6 * np.eye(A.shape[0])
    x = np.linalg.solve(A_reg, b)

Anti-Pattern 3: Inefficient Operations

# WRONG: Loop over rows/columns
result = np.zeros(n)
for i in range(n):
    result[i] = A[i, :] @ x

# RIGHT: Vectorize
result = A @ x

Anti-Pattern 4: Premature Densification

# WRONG: Converting sparse to dense unnecessarily
A_sparse = csr_matrix(...)
A_dense = A_sparse.toarray()  # Memory explosion!
x = np.linalg.solve(A_dense, b)

# RIGHT: Use sparse solvers
x = spsolve(A_sparse, b)

Troubleshooting

Issue: Singular Matrix Error

# Cause: Matrix is not invertible (det = 0 or near 0)
# Solution 1: Check rank
rank = np.linalg.matrix_rank(A)
if rank < A.shape[0]:
    print("Matrix is rank-deficient")

# Solution 2: Add regularization
A_reg = A + lambda_reg * np.eye(A.shape[0])

# Solution 3: Use pseudoinverse
x = np.linalg.pinv(A) @ b

Issue: Slow Performance

# Check if matrix is sparse
sparsity = 1 - np.count_nonzero(A) / A.size
if sparsity > 0.9:
    print("Use sparse matrix format!")
    A_sparse = csr_matrix(A)

Issue: Numerical Instability

# Check condition number
cond = np.linalg.cond(A)
if cond > 1 / np.finfo(float).eps:
    print("Matrix is numerically singular!")
    # Use more stable algorithm or regularization

Related Skills

Learning Resources