Iterative Linear Solvers

A large risk model's covariance matrix might have tens of thousands of assets, making $Ax=b$ a system with tens of thousands of equations. Gaussian elimination solves this exactly, but its cost grows like $n^3$ — for $n=50{,}000$, that's on the order of $10^{14}$ operations, far too slow. Iterative solvers take an entirely different approach: start from a rough guess for $x$, and repeatedly refine it using cheap operations (just multiplying by $A$, not inverting anything), stopping once the guess is close enough. For huge, sparse systems, this trade of exactness for speed is not just acceptable — it's the only practical option.

An analogy: tuning a radio dial by feel

Gaussian elimination is like a mechanic taking the radio apart to compute the exact dial position for a station from the circuit diagram — precise, but a lot of work per station. An iterative solver is like turning the dial, listening to how close you are to a clear signal, and nudging it a bit closer, then repeating — each nudge is cheap, and after enough small adjustments you're close enough that further tuning wouldn't be audible. You never take the radio apart, and for a radio with a thousand stations to find, this is dramatically faster overall even though it never gives you the mathematically exact dial position.

The math, one piece at a time

The simplest iterative scheme, Jacobi iteration, updates every component of the guess $x$ simultaneously using the current guess's other components:

$$x_i^{(k+1)} = \frac{1}{A_{ii}}\left(b_i - \sum_{j \neq i} A_{ij}\,x_j^{(k)}\right) .$$

In words: solve equation $i$ for $x_i$ as if every other variable's current (possibly still-wrong) guess were correct, then update all $n$ variables this way at once using last iteration's values, and repeat. Each iteration costs only a matrix-vector multiply — $O(n^2)$, or far less if $A$ is sparse (mostly zeros, common in real covariance and network structures), versus $O(n^3)$ for a full direct solve. Gauss-Seidel improves on this by using already-updated values within the same sweep instead of waiting for the next iteration, typically converging faster for the same cost per step. More advanced methods like conjugate gradient (for symmetric positive-definite $A$, the common case for covariance matrices) choose each step's direction and size more cleverly, often converging in far fewer iterations than simple sweeps.

Convergence isn't automatic: Jacobi and Gauss-Seidel are guaranteed to converge for certain well-behaved matrices (notably diagonally dominant ones, where each diagonal entry exceeds the sum of the absolute values of the rest of its row), but can diverge or converge painfully slowly otherwise. This is the trade quants make deliberately: iterative solvers gain speed and memory efficiency by giving up the direct method's unconditional guarantee.

Worked example 1: Jacobi iteration by hand

Solve $\begin{pmatrix}4&1\\1&3\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}=\begin{pmatrix}9\\8\end{pmatrix}$ (true answer: $x_1=2.2, x_2=1.9$ — check: $4(2.2)+1.9=10.7$... let's instead solve directly: from equations $4x_1+x_2=9$ and $x_1+3x_2=8$, substituting gives $x_1=2.09, x_2=1.964$ approximately). Starting guess $x^{(0)} = (0,0)$. Iteration 1: $x_1^{(1)} = \frac{9-1(0)}{4}=2.25$, $x_2^{(1)}=\frac{8-1(0)}{3}=2.667$. Iteration 2: $x_1^{(2)}=\frac{9-1(2.667)}{4}=\frac{6.333}{4}=1.583$, $x_2^{(2)}=\frac{8-1(2.25)}{3}=\frac{5.75}{3}=1.917$. Iteration 3: $x_1^{(3)}=\frac{9-1.917}{4}=1.771$, $x_2^{(3)}=\frac{8-1.583}{3}=2.139$. The guesses are visibly oscillating in but converging toward roughly $(2.09, 1.96)$ over further iterations, getting closer each sweep without ever needing to invert the matrix.

Worked example 2: when the cost saving matters

For a covariance matrix with $n=10{,}000$ assets (sparse, with most cross-asset entries near zero after thresholding), a direct solve costs on the order of $n^3 = 10^{12}$ operations. An iterative solver like conjugate gradient, needing perhaps 50-100 iterations to converge to acceptable precision, each costing roughly $O(n \times \text{nonzeros per row})$ — say 200 nonzero entries per row after sparsification — costs on the order of $100 \times 10{,}000 \times 200 = 2 \times 10^8$ operations: roughly four orders of magnitude cheaper than the direct solve, which is the entire reason large-scale risk systems use iterative solvers rather than Gaussian elimination.

What this means in practice

Iterative solvers are what makes large-scale portfolio optimization, big sparse risk-factor models, and machine-learning-adjacent linear algebra (like solving huge regularized regression systems) computationally feasible at all — direct methods simply don't scale to the sizes involved. The tradeoff to keep in mind is that an iterative solver returns an approximate answer with a chosen tolerance, not an exact one, and its convergence speed depends heavily on the matrix's condition number: a poorly-conditioned system can force far more iterations, sometimes eroding the speed advantage that motivated using an iterative method in the first place.