Linear Systems: Solving Systems of Linear Equations

1. The Problem $A\mathbf{x} = \mathbf{b}$

Given invertible $A \in \mathbb{R}^{n\times n}$, the unique solution is $\mathbf{x} = A^{-1}\mathbf{b}$. But computing $A^{-1}$ directly is expensive ($O(n^3)$) and numerically unstable — in practice this is always solved via matrix decomposition instead.

2. LU Decomposition

with $L$ lower triangular (unit diagonal) and $U$ upper triangular. Solving $A\mathbf{x}=\mathbf{b}$ becomes two easy triangular steps:

1. Solve $L\mathbf{y} = \mathbf{b}$ (forward substitution)
2. Solve $U\mathbf{x} = \mathbf{y}$ (backward substitution)

3. Cholesky Decomposition (when $A$ is symmetric positive definite)

Roughly twice as fast as general LU because it exploits symmetry. This is the standard decomposition for solving Markowitz portfolio optimization, since the covariance matrix $\Sigma$ is always symmetric and (ideally) PD.

4. Overdetermined Systems: Least Squares

When $A \in \mathbb{R}^{m\times n}$ with $m > n$ (more observations than variables — the typical regression case), the system $A\mathbf{x}=\mathbf{b}$ generally has no exact solution. We instead seek $\mathbf{x}$ that minimizes the residual:

Differentiating and setting to zero gives the normal equations:

This is exactly the OLS formula in linear regression (statistics/02-regression.md). In practice, solving via $A^TA$ can be numerically unstable when $A$ is nearly singular (multicollinearity) — it’s preferable to solve directly via QR decomposition or SVD rather than explicitly computing $(A^TA)^{-1}$.

5. Underdetermined Systems: Minimum-Norm Solution

When $m < n$ (fewer observations than variables — common when the number of factors exceeds the number of historical data points), the system has infinitely many solutions. The one with smallest $L_2$ norm is:

6. Python Implementation

import numpy as np
from scipy.linalg import lu, cho_factor, cho_solve

def solve_via_lu(A: np.ndarray, b: np.ndarray) -> np.ndarray:
    return np.linalg.solve(A, b)  # numpy internally selects LU for square systems

def solve_via_cholesky(Sigma: np.ndarray, b: np.ndarray) -> np.ndarray:
    """Efficiently solve Sigma x = b when Sigma is a covariance matrix (PD)."""
    c, low = cho_factor(Sigma)
    return cho_solve((c, low), b)

def ols_via_qr(X: np.ndarray, y: np.ndarray) -> np.ndarray:
    """OLS estimation that is more numerically stable than (A^T A)^-1 A^T y under multicollinearity."""
    Q, R = np.linalg.qr(X)
    return np.linalg.solve(R, Q.T @ y)

7. Application: Unconstrained Minimum-Variance Portfolio Weights

The unconstrained minimum-variance portfolio optimization problem reduces to solving a linear system:

In practice, $\Sigma^{-1}$ is never computed directly — instead solve $\Sigma \mathbf{z} = \mathbf{1}$ via Cholesky and then normalize $\mathbf{z}$ — both faster and numerically more stable.