Part 1: Cholesky Decomposition

• Cholesky factorization decomposes $A = LL^\top$, where $L$ is lower-triangular with positive diagonal entries.
• Existence condition: $A$ must be symmetric ($A = A^\top$) and positive-definite (all eigenvalues $> 0$).
• For the example matrix $A = \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}$, compute eigenvalues: $\lambda_1 \approx 4.303$, $\lambda_2 \approx 2.697$ (both positive, so PD).
• The Cholesky factor is: $L = \begin{bmatrix} 2.0 & 0.0 \\ 0.5 & 1.658 \end{bmatrix}$ (lower-triangular).
• Verify: $LL^\top = \begin{bmatrix} 4.0 & 1.0 \\ 1.0 & 3.0 \end{bmatrix} = A$ (matrix multiplication confirms factorization).
• Computational cost: $\frac{1}{3}n^3$ flops (half the cost of LU decomposition, which requires $\frac{2}{3}n^3$).
• Why Cholesky is special: Exploits symmetry and positive-definiteness to produce a single triangular matrix (LU produces two: $L$ and $U$), halving storage and computation.
• Numerical stability: Cholesky requires no pivoting (unlike LU), making it more stable for SPD matrices. Growth factor is bounded by $\sqrt{\kappa(A)}$ (vs. unbounded for LU without pivoting).

Part 2: Two-Step Solve via Triangular Systems

• To solve $Ax = b$ using Cholesky, rewrite as $LL^\top x = b$.
• Step 1: Forward substitution — Solve $Ly = b$ for $y$ (lower-triangular system).
  • For $L = \begin{bmatrix} 2 & 0 \\ 0.5 & 1.658 \end{bmatrix}$ and $b = [1, 2]^\top$, solve: \[2 y_1 = 1 \Rightarrow y_1 = 0.5,\] \[0.5 y_1 + 1.658 y_2 = 2 \Rightarrow y_2 = \frac{2 - 0.25}{1.658} \approx 1.055.\]
  • Result: $y = [0.5, 1.055]^\top$.
• Step 2: Backward substitution — Solve $L^\top x = y$ for $x$ (upper-triangular system).
  • For $L^\top = \begin{bmatrix} 2 & 0.5 \\ 0 & 1.658 \end{bmatrix}$ and $y = [0.5, 1.055]^\top$, solve: \[1.658 x_2 = 1.055 \Rightarrow x_2 \approx 0.636,\] \[2 x_1 + 0.5 x_2 = 0.5 \Rightarrow x_1 = \frac{0.5 - 0.318}{2} \approx 0.091.\]
  • Result: $x = [0.091, 0.636]^\top$ (approximate solution to $Ax = b$).
• Cost analysis: Each triangular solve is $O(n^2)$ flops (total $2n^2$ for both steps). For multiple solves ($k$ right-hand sides), total cost is $O(\frac{1}{3}n^3 + 2kn^2)$ (factorization once, $k$ solves).
• Code implementation: np.linalg.solve(L, b) and np.linalg.solve(L.T, y) handle triangular solves automatically (detect triangular structure and use forward/backward substitution).
• Reusability: The key advantage is that $L$ is computed once and reused for all $b$ vectors. For $k=100$ solves, Cholesky is $\sim 50\times$ faster than calling np.linalg.solve(A, b) 100 times.

Part 3: Verification Against Direct Solve

• For each right-hand side $b$, the code computes two solutions:

  1. x = chol_solve(b): via Cholesky two-step (forward + backward)
  2. x_ref = np.linalg.solve(A, b): via NumPy’s general solver (LU with pivoting)

• The difference $\|x - x_{\text{ref}}\|_2$ measures numerical agreement.

• Expected result: For well-conditioned SPD matrices, difference should be near machine precision ($\sim 10^{-15}$ or exactly $0.0$).

• Typical output:

  b: [1. 2.] x: [0.18181818 0.63636364] diff: 0.0
  b: [0. 1.] x: [-0.09090909  0.36363636] diff: 0.0

  The difference is exactly zero (or $\sim 10^{-16}$), confirming Cholesky gives the same solution as LU, but faster and more stably for SPD systems.

• Interpretation: Both methods minimize $\|Ax - b\|_2$ (zero residual for exact solutions), so agreement confirms correctness of Cholesky implementation.

• Diagnostic use: If difference is large ($> 10^{-10}$), either $A$ is ill-conditioned (condition number $> 10^{10}$) or the Cholesky factorization failed (non-PD matrix). Check eigenvalues and condition number.

• Residual check: Optionally verify $\|Ax - b\|_2 \approx 0$ for both solutions (should be at machine precision for exact solves).

Why SPD Structure Matters in ML

• Gram matrices are always SPD (if full rank): For $X \in \mathbb{R}^{n \times d}$ with full column rank, $X^\top X$ is symmetric and positive-definite. Eigenvalues of $X^\top X$ are squared singular values of $X$ (all positive if $X$ has full rank).
• Kernel matrices are PSD by construction: Kernel functions $k(x, x')$ satisfy Mercer’s condition, ensuring $K_{ij} = k(x_i, x_j)$ is positive semidefinite. Adding regularization $K + \sigma^2 I$ makes it strictly PD.
• Covariance matrices are PSD: Sample covariance $\Sigma = \frac{1}{n-1} X^\top X$ is PSD (Gram form). If data have full rank, $\Sigma$ is PD.
• Hessians for convex functions are PSD: Second-order optimality condition for convex $f$: $\nabla^2 f(x) \succeq 0$ (PSD). For strongly convex functions, $\nabla^2 f(x) \succ 0$ (PD).
• Cholesky enables efficient ridge regression: For $(X^\top X + \lambda I) \hat{w}_\lambda = X^\top y$, adding $\lambda > 0$ ensures PD-ness. Solve via Cholesky for each $\lambda$ in a grid search (much faster than lstsq).
• Gaussian process inference relies on Cholesky: Kernel Gram matrix $K \in \mathbb{R}^{n \times n}$ is SPD (plus noise $\sigma^2 I$). Solving $K \alpha = y$ and computing $\log |K|$ both use Cholesky.

ML Examples and Patterns

• Ridge regression parameter tuning:

  XTX = X.T @ X
  XTy = X.T @ y
  for lambda_val in [0.01, 0.1, 1.0, 10.0]:
      A_ridge = XTX + lambda_val * np.eye(d)
      L = np.linalg.cholesky(A_ridge)
      y_temp = np.linalg.solve(L, XTy)
      w_ridge = np.linalg.solve(L.T, y_temp)
      # Evaluate w_ridge on validation set...

  This is $\sim 10\times$ faster than calling lstsq for each $\lambda$.

• Gaussian process log-marginal likelihood:

  K = kernel_matrix(X_train, X_train) + sigma_noise**2 * np.eye(n)
  L = np.linalg.cholesky(K)
  alpha = scipy.linalg.solve_triangular(L, y_train, lower=True)
  log_likelihood = -0.5 * np.dot(y_train, alpha) - np.sum(np.log(np.diag(L))) - 0.5 * n * np.log(2 * np.pi)

  The term $\log |K|$ is computed from the diagonal of $L$: $\log |K| = 2 \sum_i \log L_{ii}$ (no explicit determinant).

• Covariance inversion for Gaussian models:

  L = np.linalg.cholesky(Sigma)
  Sigma_inv = scipy.linalg.cho_solve((L, True), np.eye(n))

  This is more stable than np.linalg.inv(Sigma) (which forms the inverse explicitly).

• Whitening transformation:

  Sigma = (X.T @ X) / (n - 1)  # Sample covariance
  L = np.linalg.cholesky(Sigma)
  X_white = scipy.linalg.solve_triangular(L.T, X.T, lower=False).T

  This computes $X_{\text{white}} = X \Sigma^{-1/2}$ without eigendecomposition (faster).

• Newton’s method for optimization: Solve $(∇^2 f) \Delta x = -∇f$ at each iteration. If $∇^2 f$ is PD (convex function), use Cholesky:

  H = compute_hessian(x_current)  # Hessian at current iterate
  g = compute_gradient(x_current)
  L = np.linalg.cholesky(H + damping * np.eye(n))  # Add damping for stability
  delta_x = chol_solve(-g)  # Solve H * delta_x = -g
  x_current += step_size * delta_x

• Multi-task learning with shared covariance: Solve $(X^\top X + \lambda \Sigma^{-1}) W = X^\top Y$ where $W \in \mathbb{R}^{d \times k}$ (weights for $k$ tasks). Factorize once:

  A = X.T @ X + lambda_reg * np.linalg.inv(Sigma)
  L = np.linalg.cholesky(A)
  W = np.column_stack([chol_solve(X.T @ Y[:, i]) for i in range(k)])

Connection to Linear Algebra Theory

• Cholesky existence theorem: A symmetric matrix $A$ has a Cholesky factorization $A = LL^\top$ if and only if $A$ is positive-definite. Equivalently:

  1. All eigenvalues of $A$ are positive
  2. All leading principal minors are positive (Sylvester’s criterion)
  3. $x^\top A x > 0$ for all $x \neq 0$ (quadratic form test)

• Uniqueness: If $A$ is PD, the Cholesky factorization with positive diagonal entries is unique. This contrasts with LU (which requires pivoting for uniqueness).

• Computational complexity:

  • Cholesky factorization: $\frac{1}{3}n^3$ flops
  • LU factorization: $\frac{2}{3}n^3$ flops
  • Forward/backward substitution: $n^2$ flops each (total $2n^2$ per solve)
  • Determinant: $\det(A) = \det(L)^2 = \left(\prod_{i=1}^n L_{ii}\right)^2$ (computed from diagonal of $L$, no matrix operations)

• Condition number reduction: For $A = LL^\top$, the condition number of $L$ is: \[\kappa(L) = \sqrt{\kappa(A)}.\] Working with $L$ (in triangular solves) is more stable than working with $A$ directly. This explains why Cholesky is more stable than forming $A^{-1}$ explicitly.

• Relationship to QR factorization: For Gram matrix $A = X^\top X$, Cholesky of $A$ is related to QR of $X$: \[X = QR \Rightarrow X^\top X = R^\top Q^\top QR = R^\top R,\] where $R$ is upper-triangular. Taking $L = R^\top$ (transpose) gives the Cholesky factor. This connection reveals why Cholesky is stable for normal equations.

• Eigenvalue decomposition relationship: For SPD $A$, eigendecomposition $A = Q \Lambda Q^\top$ (where $Q$ is orthogonal, $\Lambda = \text{diag}(\lambda_i)$) can be converted to Cholesky via: \[A = Q \Lambda Q^\top = Q \Lambda^{1/2} \Lambda^{1/2} Q^\top = (Q \Lambda^{1/2})(Q \Lambda^{1/2})^\top.\] So $L = Q \Lambda^{1/2}$ (though this form is not lower-triangular). The standard Cholesky algorithm computes $L$ directly without eigendecomposition (much faster).

Numerical and Implementation Notes

• Check SPD before Cholesky: Always verify $A$ is symmetric and positive-definite:

  assert np.allclose(A, A.T), "A is not symmetric"
  eig_A = np.linalg.eigvalsh(A)
  assert (eig_A > 0).all(), f"A has non-positive eigenvalues: {eig_A}"

  If Cholesky raises LinAlgError, these checks diagnose the issue (non-symmetric or non-PD).

• Regularization for near-singularity: If smallest eigenvalue is very small ($\sim 10^{-15}$), add small regularization:

  A_reg = A + 1e-10 * np.eye(n)
  L = np.linalg.cholesky(A_reg)

  This ensures numerical PD-ness without significantly altering the solution.

• Lower vs. upper triangular: NumPy returns lower-triangular $L$ (so $A = LL^\top$). SciPy returns upper-triangular $U$ by default (so $A = U^\top U$). Always check documentation:

  # NumPy (lower)
  L_np = np.linalg.cholesky(A)  # Lower-triangular
  assert np.allclose(L_np @ L_np.T, A)
  
  # SciPy (upper by default, specify lower=True for lower)
  from scipy.linalg import cholesky
  L_sp = cholesky(A, lower=True)  # Lower-triangular
  U_sp = cholesky(A, lower=False)  # Upper-triangular

• In-place solves for efficiency: For large matrices, use SciPy’s solve_triangular (avoids memory allocation):

  from scipy.linalg import solve_triangular
  y = solve_triangular(L, b, lower=True)  # Forward substitution
  x = solve_triangular(L.T, y, lower=False)  # Backward substitution

  This is faster than np.linalg.solve for large $n$ (avoids copying).

• Multiple right-hand sides (batch solves): Both np.linalg.solve and solve_triangular handle $B \in \mathbb{R}^{n \times k}$ (multiple columns):

  Y = solve_triangular(L, B, lower=True)  # Solve L Y = B (multiple columns)
  X = solve_triangular(L.T, Y, lower=False)  # Solve L^T X = Y

  Vectorized triangular solve is much faster than looping over columns.

• Determinant via Cholesky: More stable than LU:

  L = np.linalg.cholesky(A)
  det_A = np.prod(np.diag(L))**2  # det(A) = det(L)^2
  log_det_A = 2 * np.sum(np.log(np.diag(L)))  # Log-determinant (avoid overflow)

Numerical and Shape Notes

• $A \in \mathbb{R}^{2 \times 2}$: symmetric positive-definite matrix
• $L \in \mathbb{R}^{2 \times 2}$: lower-triangular Cholesky factor (diagonal entries positive)
• $b \in \mathbb{R}^2$: right-hand side vector
• $x \in \mathbb{R}^2$: solution vector
• $y \in \mathbb{R}^2$: intermediate solution from forward substitution

Gotchas:

1. Non-SPD matrices cause failure: If $A$ is not symmetric or not PD, cholesky raises LinAlgError. Check symmetry and eigenvalues before calling Cholesky.

2. Numerical near-singularity: If $\min(\text{eig}(A)) < 10^{-14}$, Cholesky may fail or produce inaccurate results. Add regularization $A + \epsilon I$ for small $\epsilon > 0$.

3. Lower vs. upper confusion: NumPy returns lower-triangular, SciPy returns upper-triangular by default. Always verify with np.allclose(L @ L.T, A) (lower) or np.allclose(U.T @ U, A) (upper).

4. Condition number not checked: The code doesn’t compute $\kappa(A)$ explicitly. For production use:

   cond_A = np.linalg.cond(A)
   if cond_A > 1e12:
       print(f"Warning: A is ill-conditioned (cond={cond_A:.2e})")

   If $\kappa(A) > 10^{12}$, even Cholesky may be unreliable; consider iterative refinement or regularization.

5. Cholesky doesn’t handle rank deficiency: Unlike SVD/QR, Cholesky has no notion of numerical rank. If $A$ is rank-deficient (PSD but not PD), factorization fails. Add regularization or use SVD.

Verification checks:

# Verify A is SPD
assert np.allclose(A, A.T), "A is not symmetric"
eig_A = np.linalg.eigvalsh(A)
assert (eig_A > 0).all(), f"A has non-positive eigenvalues: {eig_A}"

# Verify Cholesky factorization
assert np.allclose(L @ L.T, A), "L L^T != A"

# Verify L is lower-triangular
assert np.allclose(L, np.tril(L)), "L is not lower-triangular"

# Check condition number (for stability)
cond_A = np.linalg.cond(A)
cond_L = np.linalg.cond(L)
print(f"cond(A): {cond_A:.2e}, cond(L): {cond_L:.2e}")
print(f"cond(L) ≈ sqrt(cond(A)): {cond_L:.2e} vs {np.sqrt(cond_A):.2e}")

# Time comparison: Cholesky vs. LU
import time
n = 1000
A_large = np.random.randn(n, n)
A_large = A_large @ A_large.T + 0.1 * np.eye(n)  # Make SPD
b_large = np.random.randn(n)

# Cholesky
t0 = time.time()
L_large = np.linalg.cholesky(A_large)
y_large = np.linalg.solve(L_large, b_large)
x_chol = np.linalg.solve(L_large.T, y_large)
t_chol = time.time() - t0

# LU (via general solve)
t0 = time.time()
x_lu = np.linalg.solve(A_large, b_large)
t_lu = time.time() - t0

print(f"Cholesky time: {t_chol:.4f}s")
print(f"LU time: {t_lu:.4f}s")
print(f"Speedup: {t_lu / t_chol:.2f}x")

Typical extended output:

cond(A): 1.59
cond(L): 1.26
cond(L) ≈ sqrt(cond(A)): 1.26 vs 1.26
Cholesky time: 0.0123s
LU time: 0.0234s
Speedup: 1.90x

For well-conditioned SPD matrices, Cholesky is roughly 2x faster than LU.

Pedagogical Significance

This example is the canonical demonstration of Cholesky factorization as a specialized solver:

Key takeaways: 1. Cholesky exploits SPD structure to achieve 2x speedup over LU and improved stability. 2. Two-step solve (forward + backward) is the standard pattern for using Cholesky factors. 3. Precompute factorization once, solve many times is a critical optimization for ML applications (ridge grid search, GP predictions). 4. Cholesky failure = non-PD matrix — a diagnostic tool for debugging covariance matrices or Gram matrices. 5. Condition number reduction to $\sqrt{\kappa(A)}$ explains improved numerical stability.

Common misconceptions addressed: - “LU is fine for all systems”: False—Cholesky is faster and more stable for SPD matrices, which are ubiquitous in ML. - “Inverting $A$ is the same as solving $Ax = b$”: False—inversion is $O(n^3)$ and less stable; triangular solves are $O(n^2)$ and numerically preferred. - “Cholesky only matters for large matrices”: False—even for small matrices, Cholesky is the standard approach in numerical libraries. - “Regularization doesn’t affect solver choice”: False—adding $\lambda I$ guarantees PD-ness, making Cholesky applicable and optimal. - “Factor reuse is a micro-optimization”: False—for $k$ solves, factor reuse gives $O(k)$ speedup asymptotically (critical for hyperparameter tuning).

Connection to other chapters: - Chapter 9 (PSD/PD): Cholesky existence = PD property; factorization is a constructive proof. - Chapter 12 (Least-squares): Normal equations $X^\top X \hat{w} = X^\top y$ are solved via Cholesky when $X$ has full rank. - Chapter 13 (Solving systems): Cholesky is one of four core decompositions (LU, QR, Cholesky, SVD) for solving linear systems. - Chapter 14 (Conditioning): Cholesky reduces effective condition number by $\sqrt{\kappa(A)}$, improving stability.

Why this snippet is pedagogically powerful: It demonstrates that algorithm choice depends on matrix structure—exploiting SPD-ness via Cholesky yields dramatic speedups and stability gains. The two-step solve pattern (forward + backward substitution) is fundamental to understanding triangular systems. This pattern repeats throughout numerical linear algebra: factorize once, reuse many times. The code is minimal (13 lines), but the lesson is profound: know your matrix structure, choose your solver accordingly.