Mastering Matrix Decompositions and SVD: LU, QR, Eigendecomposition, and Singular Value Decomposition

1. Why Factorize Matrices? Computational Efficiency & Numerical Stability

Solving a linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$ for an $N imes N$ matrix $\mathbf{A}$ using standard Gaussian elimination or explicit matrix inversion ($\mathbf{x} = \mathbf{A}^{-1} \mathbf{b}$) requires $O(N^3)$ floating-point operations. Worse, explicit matrix inversion is numerically unstable due to floating-point round-off amplification under IEEE 754 double-precision arithmetic. When a matrix is ill-conditioned, small measurement errors in input vector $\mathbf{b}$ cause massive catastrophic cancellation errors in output $\mathbf{x}$.

By factoring $\mathbf{A}$ into simpler structural forms (such as lower triangular $\mathbf{L}$, upper triangular $\mathbf{U}$, orthogonal $\mathbf{Q}$, or diagonal $\mathbf{\Sigma}$ matrices), solving $\mathbf{A} \mathbf{x} = \mathbf{b}$ reduces to forward and back-substitution in $O(N^2)$ time! Furthermore, factorizations guarantee numerical stability by maintaining tight condition bounds throughout computation.

2. LU Decomposition: Doolittle & Gaussian Elimination

LU Decomposition factors a square matrix $\mathbf{A}$ into a product of a Lower Triangular matrix $\mathbf{L}$ (with 1s on the main diagonal) and an Upper Triangular matrix $\mathbf{U}$:

2.1 Step-by-Step Worked Trace (2x2 Matrix)

Let's factor matrix $\mathbf{A} = egin{pmatrix} 2 & 4 \ 3 & 8 \end{pmatrix}$ manually:

Step 1: Set up matrix multiplication terms:
  u_11 = a_11 = 2
  u_12 = a_12 = 4

Step 2: Solve for L term l_21:
  l_21 * u_11 = a_21 => l_21 * 2 = 3 => l_21 = 1.5

Step 3: Solve for U term u_22:
  l_21 * u_12 + u_22 = a_22 => 1.5 * 4 + u_22 = 8 => 6 + u_22 = 8 => u_22 = 2

Result Matrices:
  L = [[1.0, 0.0], [1.5, 1.0]]
  U = [[2.0, 4.0], [0.0, 2.0]]

Verification L * U:
  [1.0*2.0 + 0,  1.0*4.0 + 0]       = [2.0, 4.0]
  [1.5*2.0 + 2,  1.5*4.0 + 1*2]     = [3.0 + 0, 6.0 + 2.0] = [3.0, 8.0] (Matches A!)

2.2 Python Implementation of Partial Pivoting LU Decomposition

To prevent division by zero when $a_{ii} = 0$, we use Partial Pivoting permutation matrix $\mathbf{P}$ such that $\mathbf{P} \mathbf{A} = \mathbf{L} \mathbf{U}$:

import numpy as np

def lu_decomposition(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()
    P = np.eye(n)

    for i in range(n):
        # Pivot selection (find max element in column i)
        max_row = i + np.argmax(np.abs(U[i:, i]))
        if i != max_row:
            U[[i, max_row]] = U[[max_row, i]]
            P[[i, max_row]] = P[[max_row, i]]
            if i > 0:
                L[[i, max_row], :i] = L[[max_row, i], :i]

        for j in range(i + 1, n):
            factor = U[j, i] / U[i, i]
            L[j, i] = factor
            U[j, i:] -= factor * U[i, i:]

    return P, L, U

3. QR Decomposition: Gram-Schmidt & Householder Reflections

QR Decomposition factors an $M imes N$ matrix $\mathbf{A}$ into an Orthogonal matrix $\mathbf{Q} \in \mathbb{R}^{M imes M}$ (where $\mathbf{Q}^T \mathbf{Q} = \mathbf{I}$) and an Upper Triangular matrix $\mathbf{R} \in \mathbb{R}^{M imes N}$:

Because $\mathbf{Q}$ is orthogonal, its inverse is simply its transpose ($\mathbf{Q}^{-1} = \mathbf{Q}^T$). Solving least-squares problems $\mathbf{A} \mathbf{x} pprox \mathbf{b}$ simplifies to $\mathbf{R} \mathbf{x} = \mathbf{Q}^T \mathbf{b}$!

3.1 Gram-Schmidt Orthogonalization Process

Given column vectors $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n$ of matrix $\mathbf{A}$, Gram-Schmidt computes orthogonal basis vectors $\mathbf{u}_1, \mathbf{u}_2, \dots$ and normalizes them into unit vectors $\mathbf{q}_1, \mathbf{q}_2, \dots$:

def qr_gram_schmidt(A):
    m, n = A.shape
    Q = np.zeros((m, n))
    R = np.zeros((n, n))

    for k in range(n):
        v = A[:, k].astype(float)
        for j in range(k):
            R[j, k] = np.dot(Q[:, j], A[:, k])
            v -= R[j, k] * Q[:, j]
        R[k, k] = np.linalg.norm(v)
        Q[:, k] = v / R[k, k]

    return Q, R

4. Eigendecomposition: Eigenvalues and Eigenvectors

For a square $N imes N$ matrix $\mathbf{A}$, an eigenvector $\mathbf{v}$ is a non-zero vector that changes only in scale (not direction) when linear transformation $\mathbf{A}$ is applied:

where scalar $\lambda$ is the corresponding eigenvalue. If $\mathbf{A}$ has $N$ linearly independent eigenvectors, we assemble them into matrix $\mathbf{V} = [\mathbf{v}_1, \dots, \mathbf{v}_n]$ and diagonal matrix $\mathbf{\Lambda} = ext{diag}(\lambda_1, \dots, \lambda_n)$ to factor $\mathbf{A}$:

5. Singular Value Decomposition (SVD): Geometric Intuition

While Eigendecomposition is strictly limited to square matrices, Singular Value Decomposition (SVD) applies to any arbitrary $M imes N$ matrix. SVD states that any linear transformation can be broken down into three fundamental geometric operations: a rotation, a scaling, and a second rotation!

Geometrically, matrix $\mathbf{A}$ transforms a unit sphere in $\mathbb{R}^N$ into an hyper-ellipsoid in $\mathbb{R}^M$. The right singular vectors $\mathbf{v}_i$ form the orthogonal directions in the domain that map to the principal axes of the ellipsoid. The singular values $\sigma_i$ measure the exact lengths of these ellipsoid axes, while the left singular vectors $\mathbf{u}_i$ define the direction of these axes in the codomain space!

5.1 Mathematical Breakdown of SVD Components

• Left Singular Vectors ($\mathbf{U} \in \mathbb{R}^{M imes M}$): Orthogonal eigenvectors of $\mathbf{A} \mathbf{A}^T$. Form the orthonormal basis for output space.
• Singular Values Matrix ($\mathbf{\Sigma} \in \mathbb{R}^{M imes N}$): Rectangular diagonal matrix containing non-negative singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r \ge 0$. Note that $\sigma_i = \sqrt{\lambda_i(\mathbf{A}^T \mathbf{A})}$.
• Right Singular Vectors ($\mathbf{V} \in \mathbb{R}^{N imes N}$): Orthogonal eigenvectors of $\mathbf{A}^T \mathbf{A}$. Form the orthonormal basis for input domain space.

6. Numerical Step-by-Step Worked Trace: Computing SVD of a 2x2 Matrix

Let's manually compute the full SVD of matrix $\mathbf{A} = egin{pmatrix} 3 & 0 \ 4 & 5 \end{pmatrix}$:

Step 1: Compute A^T * A
  A^T = [[3, 4], [0, 5]]
  A^T * A = [[3, 4], [0, 5]] * [[3, 0], [4, 5]] = [[25, 20], [20, 25]]

Step 2: Find Eigenvalues of A^T * A
  det(A^T*A - lambda*I) = (25 - lambda)^2 - 400 = 0
  (25 - lambda) = +/- 20
  lambda_1 = 45,  lambda_2 = 5

Step 3: Calculate Singular Values sigma_i = sqrt(lambda_i)
  sigma_1 = sqrt(45) ≈ 6.7082
  sigma_2 = sqrt(5)  ≈ 2.2361
  Sigma = [[6.7082, 0.0], [0.0, 2.2361]]

Step 4: Find Eigenvectors of A^T * A (Matrix V)
  For lambda_1 = 45:
    [[25-45, 20], [20, 25-45]] * v1 = [[-20, 20], [20, -20]] * v1 = 0 => v1 = [1/sqrt(2), 1/sqrt(2)]
  For lambda_2 = 5:
    [[20, 20], [20, 20]] * v2 = 0 => v2 = [-1/sqrt(2), 1/sqrt(2)]
  V = [[0.7071, -0.7071], [0.7071, 0.7071]]

Step 5: Compute Left Singular Vectors u_i = (1 / sigma_i) * A * v_i
  u_1 = (1 / 6.7082) * [[3, 0], [4, 5]] * [0.7071, 0.7071]^T = [0.3162, 0.9487]^T
  u_2 = (1 / 2.2361) * [[3, 0], [4, 5]] * [-0.7071, 0.7071]^T = [-0.9487, 0.3162]^T
  U = [[0.3162, -0.9487], [0.9487, 0.3162]]

Final Check U * Sigma * V^T matches original matrix A!

7. Low-Rank Matrix Approximation & Image Compression

According to the Eckart-Young-Mirsky Theorem, the optimal rank-$k$ approximation $\mathbf{A}_k$ of matrix $\mathbf{A}$ (minimizing Frobenius norm error $\|\mathbf{A} - \mathbf{A}_k\|_F$) is obtained by keeping only the top-$k$ largest singular values and discarding the rest:

# Image Compression using Truncated SVD in Python
import numpy as np

def compress_matrix_svd(A, k_components):
    U, s, Vt = np.linalg.svd(A, full_matrices=False)
    # Truncate to top k components
    U_k = U[:, :k_components]
    s_k = s[:k_components]
    Vt_k = Vt[:k_components, :]
    
    # Reconstruct rank-k matrix
    A_compressed = np.dot(U_k * s_k, Vt_k)
    
    original_size = A.shape[0] * A.shape[1]
    compressed_size = k_components * (A.shape[0] + A.shape[1] + 1)
    compression_ratio = original_size / compressed_size
    
    return A_compressed, compression_ratio

For a 2000x2000 image, storing 50 singular components requires storing $50 imes (2000 + 2000 + 1) = 200,050$ numbers instead of 4,000,000—achieving a **20x compression ratio** with minimal visual loss!

8. Principal Component Analysis (PCA) via SVD

Principal Component Analysis (PCA) is the workhorse of unsupervised dimensionality reduction. Given a mean-centered data matrix $\mathbf{X} \in \mathbb{R}^{N imes D}$ ($N$ samples, $D$ features), the sample covariance matrix is $\mathbf{C} = rac{1}{N-1} \mathbf{X}^T \mathbf{X}$.

Instead of explicitly constructing the $D imes D$ covariance matrix $\mathbf{C}$ (which is slow for huge feature dimensions), we compute SVD directly on $\mathbf{X}$:

The right singular vectors $\mathbf{V}$ are exact principal component directions, and the singular values give the variance explained along each axis: $ ext{Var}(PC_i) = rac{\sigma_i^2}{N-1}$.

9. Computing the Moore-Penrose Pseudoinverse ($\mathbf{A}^+$)

When matrix $\mathbf{A}$ is non-square ($M e N$) or singular, its standard inverse $\mathbf{A}^{-1}$ does not exist. The Moore-Penrose Pseudoinverse $\mathbf{A}^+$ generalizes matrix inversion using SVD:

where $\mathbf{\Sigma}^+$ is formed by taking the reciprocal of each non-zero singular value ($1/\sigma_i$) and transposing the matrix shape. The pseudoinverse provides the exact minimum-norm solution to any overdetermined or underdetermined linear system $\mathbf{A} \mathbf{x} = \mathbf{b}$!

10. Randomized SVD for Huge Datasets

Standard full SVD algorithms (like Golub-Reinsch) require $O(M N \min(M, N))$ operations, making them unfeasible for massive datasets (e.g. $100,000 imes 100,000$ matrices). **Randomized SVD** uses random projection matrices $\mathbf{\Omega}$ to sample the range space of $\mathbf{A}$, projecting it to a small matrix of target rank $k$:

def randomized_svd(A, k, n_iter=2):
    m, n = A.shape
    # Generate Gaussian random matrix
    Omega = np.random.normal(size=(n, k + 10))
    
    # Form sample matrix Y and perform Power Iteration to amplify dominant spectra
    Y = A @ Omega
    for _ in range(n_iter):
        Y = A @ (A.T @ Y)
        
    # Orthonormalize Y to find Q
    Q, _ = np.linalg.qr(Y)
    
    # Project A into low-dimensional space
    B = Q.T @ A
    
    # Compute standard SVD on small matrix B
    U_hat, s, Vt = np.linalg.svd(B, full_matrices=False)
    U = Q @ U_hat
    
    return U[:, :k], s[:k], Vt[:k, :]

Randomized SVD computes top-$k$ singular components in $O(M N \log k)$ time—speeding up computation by over 50x for large sparse matrices!

11. Recommender Systems: SVD Collaborative Filtering

In recommender systems (such as the famous Netflix Prize algorithm), a User-Item rating matrix $\mathbf{R} \in \mathbb{R}^{U imes I}$ is partially observed. SVD factorizes $\mathbf{R}$ into user preference vectors $\mathbf{P}_u$ and item feature vectors $\mathbf{Q}_i$:

By optimizing user and item bias parameters along with latent vectors using Stochastic Gradient Descent (SGD), models predict unobserved ratings to provide personalized recommendations!

12. Cholesky Decomposition for Positive Definite Matrices

When matrix $\mathbf{A}$ is Symmetric Positive-Definite ($\mathbf{x}^T \mathbf{A} \mathbf{x} > 0$), Cholesky Decomposition factors $\mathbf{A}$ into a single lower triangular matrix $\mathbf{L}$:

Cholesky decomposition requires only $N^3 / 3$ FLOPs—twice as fast as standard LU decomposition! It is extensively used in Gaussian Process Regression, Kalman Filtering, and Monte Carlo risk simulations.

13. Non-Negative Matrix Factorization (NMF)

Unlike SVD, which allows negative entries in $\mathbf{U}$ and $\mathbf{V}$, Non-Negative Matrix Factorization (NMF) enforces strictly non-negative constraints ($\mathbf{W} \ge 0, \mathbf{H} \ge 0$):

In document topic modeling or audio spectrogram separation, non-negativity ensures pure additive parts-based representations, making topics easily interpretable by human analysts.

14. Polar Decomposition: Rotation and Stretch Separation

Every square matrix $\mathbf{A} \in \mathbb{R}^{N imes N}$ can be factored into a Polar Decomposition $\mathbf{A} = \mathbf{U}_p \mathbf{P}_s$, where $\mathbf{U}_p$ is an orthogonal rotation matrix and $\mathbf{P}_s = \sqrt{\mathbf{A}^T \mathbf{A}}$ is a symmetric positive semi-definite stretch matrix. Polar decomposition can be computed directly from SVD by setting $\mathbf{U}_p = \mathbf{U} \mathbf{V}^T$ and $\mathbf{P}_s = \mathbf{V} \mathbf{\Sigma} \mathbf{V}^T$, providing applications in 3D computer graphics physics engines!

15. SVD in Quantum Information (Schmidt Decomposition)

In quantum mechanics, a bipartite pure state $|\psi angle \in \mathcal{H}_A \otimes \mathcal{H}_B$ can be expanded using SVD into its Schmidt Decomposition:

where coefficients $\sqrt{\lambda_i}$ are singular values of the state coefficient matrix. SVD directly quantifies quantum entanglement entropy through von Neumann entropy $S = -\sum \lambda_i \log_2 \lambda_i$!

16. Tensor Decompositions: CP and Tucker Factorizations

When data extends beyond 2D matrices into multi-dimensional arrays (tensors $\mathcal{X} \in \mathbb{R}^{I imes J imes K}$), SVD generalizes into two primary tensor decompositions:

• CP (Canonical Polyadic) Decomposition: Factorizes a tensor into a sum of rank-1 vector outer products: $\mathcal{X} pprox \sum_{r=1}^R \mathbf{a}_r \circ \mathbf{b}_r \circ \mathbf{c}_r$.
• Tucker Decomposition (Higher-Order SVD): Factorizes a tensor into a small core tensor $\mathcal{G}$ multiplied by factor matrices along each mode: $\mathcal{X} pprox \mathcal{G} imes_1 \mathbf{A} imes_2 \mathbf{B} imes_3 \mathbf{C}$.

17. Proof of Eckart-Young-Mirsky Low-Rank Approximation Theorem

Let $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T = \sum_{i=1}^r \sigma_i \mathbf{u}_i \mathbf{v}_i^T$ be the full rank-$r$ SVD. We want to prove that the rank-$k$ matrix $\mathbf{A}_k = \sum_{i=1}^k \sigma_i \mathbf{u}_i \mathbf{v}_i^T$ minimizes $\|\mathbf{A} - \mathbf{B}\|_F^2$ over all matrices $\mathbf{B}$ of rank at most $k$.

By unitarily invariant properties of the Frobenius norm:

Since singular values are ordered $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r$, discarding the smallest $r-k$ singular values minimizes the residual sum of squares error, completing the proof!

18. Golub-Kahan Bidiagonalization SVD Algorithm

The standard LAPACK algorithm for computing SVD does not work directly on dense matrices. Instead, it proceeds in two distinct phases:

1. Bidiagonalization: Applies Householder reflections from left and right to reduce $\mathbf{A}$ into an upper bidiagonal matrix $\mathbf{B}$.
2. Implicit QR Steps: Computes the singular values of bidiagonal matrix $\mathbf{B}$ using implicit QR steps in $O(N^2)$ time.

19. Power Iteration and Rayleigh Quotient Iteration in Python

To find the dominant eigenvalue $\lambda_{\max}$ and eigenvector $\mathbf{v}_1$ of a symmetric matrix $\mathbf{A}$, we use Power Iteration:

def power_iteration(A, num_iterations=100):
    b_k = np.random.rand(A.shape[1])
    for _ in range(num_iterations):
        # Calculate matrix-vector product
        b_k1 = A @ b_k
        # Calculate norm
        b_k = b_k1 / np.linalg.norm(b_k1)
    
    # Rayleigh quotient for eigenvalue estimate
    eigenvalue = (b_k.T @ A @ b_k) / (b_k.T @ b_k)
    return eigenvalue, b_k

20. Matrix Completion via Soft Thresholding Singular Values

In matrix completion problems with missing entries (such as medical imaging or genomics), low-rank matrix recovery solves:

The Singular Value Thresholding (SVT) algorithm iteratively updates $\mathbf{X}_k = \mathcal{D}_{ au}(\mathbf{Y}_k)$ where $\mathcal{D}_{ au}(\mathbf{Y}) = \mathbf{U} ext{diag}(\max(\sigma_i - au, 0)) \mathbf{V}^T$, providing exact low-rank reconstruction!

21. Latent Semantic Indexing (LSI) Search Engine Matching

In search engines, Latent Semantic Indexing (LSI) uses SVD to build a latent concept space. A search query vector $\mathbf{q} \in \mathbb{R}^V$ is projected into $k$-dimensional latent space via $\mathbf{q}_k = \mathbf{q}^T \mathbf{U}_k \mathbf{\Sigma}_k^{-1}$. Document relevance is calculated via cosine similarity against stored document vectors $\mathbf{V}_k$:

22. Homography Estimation and DLT Solvers in Computer Vision

In 3D computer vision and image stitching, finding the $3 imes 3$ Homography matrix $\mathbf{H}$ that maps points $\mathbf{x}_i \leftrightarrow \mathbf{x}'_i$ requires solving an overdetermined homogeneous system $\mathbf{A} \mathbf{h} = \mathbf{0}$ subject to $\|\mathbf{h}\| = 1$.

The Direct Linear Transform (DLT) algorithm computes SVD $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$. The optimal solution vector $\mathbf{h}$ is the last column of $\mathbf{V}$ corresponding to the smallest singular value $\sigma_{\min}$, providing maximum numerical stability against camera noise!

23. Generalized Singular Value Decomposition (GSVD)

Given two matrices $\mathbf{A} \in \mathbb{R}^{M imes N}$ and $\mathbf{B} \in \mathbb{R}^{P imes N}$ sharing the same number of columns, Generalized SVD (GSVD) decomposes both matrices simultaneously using a single non-singular matrix $\mathbf{X}$ and orthogonal matrices $\mathbf{U}$ and $\mathbf{V}$:

GSVD is the core mathematical engine behind comparative gene expression profiling across biological species and joint signal detection in multi-user MIMO telecommunication arrays!

24. PyTorch Custom Autograd SVD Layer Function

Below is a custom PyTorch Autograd layer that applies Truncated SVD during the forward pass while enabling exact gradient propagation during backpropagation:

class TruncatedSVDLayer(torch.autograd.Function):
    @staticmethod
    def forward(ctx, X, k):
        U, S, Vh = torch.linalg.svd(X, full_matrices=False)
        U_k = U[:, :k]
        S_k = S[:k]
        Vh_k = Vh[:k, :]
        
        ctx.save_for_backward(U_k, S_k, Vh_k)
        return U_k @ torch.diag(S_k) @ Vh_k

    @staticmethod
    def backward(ctx, grad_output):
        U_k, S_k, Vh_k = ctx.saved_tensors
        # Backward gradient propagation for low-rank matrix
        grad_input = grad_output @ Vh_k.T @ torch.diag(1.0 / (S_k + 1e-6)) @ U_k.T
        return grad_input, None

25. Eigenfaces Facial Recognition Pipeline in Python

The classic Eigenfaces face recognition algorithm applies SVD to a database of face images to compute principal face components:

def compute_eigenfaces(face_images_matrix, n_components=30):
    # face_images_matrix shape: (N_samples, H * W)
    mean_face = np.mean(face_images_matrix, axis=0)
    centered_faces = face_images_matrix - mean_face
    
    # Compute SVD on centered face matrix
    U, S, Vt = np.linalg.svd(centered_faces, full_matrices=False)
    
    eigenfaces = Vt[:n_components, :]
    weights = centered_faces @ eigenfaces.T
    
    return mean_face, eigenfaces, weights

26. Dynamic Mode Decomposition (DMD) in Fluid Dynamics

In aerospace engineering and fluid dynamics modeling, Dynamic Mode Decomposition (DMD) extracts coherent spatial-temporal structures from experimental flow field snapshots. Given snapshot matrices $\mathbf{X}_1$ and $\mathbf{X}_2$, DMD computes SVD $\mathbf{X}_1 = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$ to construct a reduced-order system operator $ ilde{\mathbf{A}} = \mathbf{U}^T \mathbf{X}_2 \mathbf{V} \mathbf{\Sigma}^{-1}$, capturing fluid turbulence frequencies and growth rates!

27. SVD Digital Watermarking for Copyright Forensics

Digital watermarking embeds copyright metadata directly into singular values of image matrices. By modifying singular values $\mathbf{\Sigma}' = \mathbf{\Sigma} + lpha \mathbf{\Sigma}_w$, the watermark remains invisible to human vision while surviving compression, spatial scaling, and adversarial noise attacks.

28. Householder Reflection Transformation Matrices

A Householder Reflection matrix $\mathbf{H} = \mathbf{I} - 2 \mathbf{v} \mathbf{v}^T / (\mathbf{v}^T \mathbf{v})$ reflects any vector $\mathbf{x}$ across a hyperplane orthogonal to vector $\mathbf{v}$. Householder reflections are used in QR decomposition and bidiagonalization because they are perfectly orthogonal ($\mathbf{H}^T \mathbf{H} = \mathbf{I}$) and preserve vector norms with zero round-off drift!

29. Audio Subspace Noise Reduction via SVD

In audio signal processing, a noisy signal matrix $\mathbf{Y} = \mathbf{X} + \mathbf{N}$ constructed from Hankel time-delay frames is decomposed via SVD. Signals reside in dominant singular components, while white noise is uniformly spread across small singular values. Zeroing small singular values recovers clean speech waveforms!

30. Complete Python Performance Benchmarking Suite

Below is a benchmarking script comparing standard LAPACK SVD against Randomized SVD and Gram-Schmidt QR decomposition across matrix sizes up to $5000 imes 5000$:

import time
import numpy as np

def benchmark_matrix_decompositions():
    sizes = [500, 1000, 2000]
    print(f"{'Size':<8} | {'Exact SVD (s)':<15} | {'Rand SVD k=50 (s)':<18} | {'LU Dec (s)':<12}")
    print("-" * 60)

    for n in sizes:
        A = np.random.randn(n, n)
        
        # 1. Exact SVD
        t0 = time.time()
        np.linalg.svd(A, full_matrices=False)
        t_svd = time.time() - t0
        
        # 2. Randomized SVD
        t0 = time.time()
        randomized_svd(A, k=50)
        t_rsvd = time.time() - t0
        
        # 3. LU Decomposition
        t0 = time.time()
        lu_decomposition(A)
        t_lu = time.time() - t0
        
        print(f"{n:<8} | {t_svd:<15.4f} | {t_rsvd:<18.4f} | {t_lu:<12.4f}")

benchmark_matrix_decompositions()

31. Developer Pitfall Box

32. Production Engineering Summary Checklist

• Fast Linear Systems: Use LU decomposition when solving $\mathbf{A} \mathbf{x} = \mathbf{b}$ for multiple right-hand side vectors $\mathbf{b}$.
• Least-Squares Stability: Prefer QR decomposition over normal equations $(\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{b}$ for linear regression.
• Dimensionality Reduction: Use Truncated SVD or Randomized SVD for PCA on high-dimensional sparse matrices.
• Avoid Direct Inversion: Never compute np.linalg.inv(A); use np.linalg.solve(A, b) or pseudoinverse np.linalg.pinv(A).

33. Developer FAQ

Written by Professor Pixel · CodingPancake Systems Architecture Series