Brownian Bridge™

Setup

Direct solvers — LU factorisation (Module 2) and Cholesky (Module 2) — cost $O(n^3)$ flops and require $O(n^2)$ storage. For a Crank-Nicolson finite difference grid with $N = 10^4$ spatial nodes, forming and factorising the full stiffness matrix is $10^{12}$ flops and 800 MB of storage — impractical. For a 2-factor interest rate PDE on a $200 \times 200$ grid, the system is $n = 4 \times 10^4$ — far beyond direct methods.

Iterative solvers exploit the structure of the system: they require only matrix-vector products $Ap$ , never the factorisation. For sparse matrices (finite difference stencils, sparse covariance models), each $Ap$ costs $O(n)$ rather than $O(n^2)$ . The cost of the solve is $O(n \cdot k_\text{iter})$ where $k_\text{iter}$ depends on the condition number.

Where iterative solvers appear in quant finance:

PDE pricing (finite differences). Black-Scholes and Heston PDEs on dense grids produce large tridiagonal or banded systems at each time step. For 1D grids the Thomas algorithm (a direct tridiagonal solver, $O(n)$ ) is standard; for 2D grids, alternating direction implicit (ADI) splits into 1D tridiagonals. For 3D or coupled systems, iterative solvers (PCG) are used.
Yield curve calibration. Bootstrapping a swap curve with $n$ pillars gives an $n \times n$ lower-triangular system (direct, trivial). Fitting a regularised Nelson-Siegel-Svensson model to $m \gg n$ market quotes leads to a normal-equations solve $(A^\top A + \lambda I)x = A^\top b$ — either direct (via Cholesky) or iterative (CG) depending on $n$ .
Portfolio optimisation. Mean-variance with $n = 500$ assets and $k$ linear constraints gives a KKT system of dimension $n + k \approx 600$ . Cholesky is fine here. But factor-model covariance $\Sigma = F F^\top + D$ (low-rank plus diagonal) inverts analytically via Woodbury — no solver needed.
Large-scale calibration of stochastic vol models. Iterative refinement is used when the Jacobian is large and only approximately known (numerical differentiation), making direct factorisation unreliable.

Mathematical setting. We focus on three cases:

SPD system $Ax = b$ with $A$ sparse and symmetric positive definite: solved by Conjugate Gradient (CG).
Non-symmetric system: GMRES (Generalised Minimal RESidual).
Overdetermined LS $\min \|Ax - b\|_2$ : LSQR.

Conventions. $A \in \mathbb{R}^{n \times n}$ SPD unless stated. $r_k = b - Ax_k$ is the residual at step $k$ . $e_k = x^* - x_k$ is the error. $\|v\|_A = \sqrt{v^\top A v}$ is the A-norm (energy norm).

INSIGHT

Financial Insight. In production calibration at a bank, iterative solvers are used not just for speed but for robustness: they can be stopped early (before convergence) when the residual is below market-quote precision. Solving to $10^{-10}$ accuracy when quotes are only known to $10^{-4}$ wastes computation. Early stopping is equivalent to a form of regularisation — an important practical insight.

Theory

1. Quadratic Minimisation and the CG Motivation

For $A$ SPD, solving $Ax = b$ is equivalent to minimising the strictly convex quadratic: $f(x) = \frac{1}{2} x^\top A x - b^\top x.$ The gradient is $\nabla f(x) = Ax - b = -r(x)$ . The minimum is at $x^* = A^{-1}b$ .

Steepest descent. The simplest iterative approach: move in the direction of the residual $r_k = b - Ax_k$ by the optimal step size:

DEFINITION

Definition 5.1 (Steepest Descent). Given $x_0$ , iterate: $\alpha_k = \frac{r_k^\top r_k}{r_k^\top A r_k}, \qquad x_{k+1} = x_k + \alpha_k r_k.$ The step $\alpha_k$ is chosen to minimise $f(x_k + \alpha r_k)$ over $\alpha \in \mathbb{R}$ .

Convergence of steepest descent:

THEOREM

Theorem 5.1. Let $\kappa = \kappa_2(A) = \lambda_\text{max}/\lambda_\text{min}$ . Steepest descent satisfies: $\|e_{k+1}\|_A \leq \left(\frac{\kappa - 1}{\kappa + 1}\right)^2 \|e_k\|_A.$ For ill-conditioned problems ( $\kappa \gg 1$ ), the convergence factor approaches 1 — very slow.

The fundamental problem with steepest descent: consecutive search directions $r_{k+1} \perp r_k$ in the Euclidean inner product, but they are not $A$ -orthogonal. The method revisits the same subspaces, causing "zigzag" convergence along narrow valleys of $f$ when eigenvalues are spread widely.

2. Conjugate Gradient Method

The CG method eliminates the redundant search directions by enforcing A-conjugacy:

DEFINITION

Definition 5.2 (A-conjugate directions). Vectors $p, q$ are A-conjugate (A-orthogonal) if $p^\top A q = 0$ . A set $\{p_0, \ldots, p_{k-1}\}$ is A-conjugate if $p_i^\top A p_j = 0$ for all $i \neq j$ .

THEOREM

Theorem 5.2 (CG Algorithm — Hestenes-Stiefel, 1952). Starting from $x_0$ , $r_0 = b - Ax_0$ , $p_0 = r_0$ , iterate: $\alpha_k = \frac{r_k^\top r_k}{p_k^\top A p_k}, \quad x_{k+1} = x_k + \alpha_k p_k, \quad r_{k+1} = r_k - \alpha_k A p_k,$ $\beta_k = \frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}, \quad p_{k+1} = r_{k+1} + \beta_k p_k.$

Key properties (proofs by induction on the Krylov subspace structure):

$r_i \perp r_j$ for $i \neq j$ (residuals are mutually orthogonal in the Euclidean inner product).
$p_i^\top A p_j = 0$ for $i \neq j$ (search directions are A-conjugate).
$x_k$ minimises $f(x) = \frac{1}{2}x^\top Ax - b^\top x$ over the Krylov subspace $\mathcal{K}_k = \text{span}\{r_0, Ar_0, \ldots, A^{k-1}r_0\}$ .
In exact arithmetic, CG terminates in at most $n$ steps with the exact solution.

CG convergence rate:

THEOREM

Theorem 5.3 (CG convergence). The CG iterates satisfy: $\|e_k\|_A \leq 2\left(\frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1}\right)^k \|e_0\|_A.$ For $\kappa = 100$ : CG convergence factor $\approx 0.82$ ; steepest descent factor $\approx 0.96$ .

Comparison: CG depends on $\sqrt{\kappa}$ rather than $\kappa$ . For $\kappa = 10^4$ (typical calibration matrix):

Steepest descent: $\approx 10^4$ iterations to reduce error by $10^{-4}$ .
CG: $\approx 100$ iterations (100× faster).

EXAMPLE

Example 5.1 (Tridiagonal FD system). A Crank-Nicolson price step for a $n = 500$ node grid gives a tridiagonal SPD matrix $A$ with $\kappa \approx 4n^2/\pi^2 \approx 10^5$ (for the standard 1D heat equation stencil). CG converges in $\approx \sqrt{10^5} \approx 316$ iterations to machine precision — competitive with the Thomas algorithm ( $O(n) = O(500)$ steps). For 2D ADI, the subproblems are 1D tridiagonals, so Thomas algorithm dominates.

3. Preconditioning

The convergence rate of CG depends entirely on $\kappa_2(A)$ . Preconditioning replaces the original system with a better-conditioned equivalent.

DEFINITION

Definition 5.3 (Preconditioned CG). Given a preconditioner $M \approx A$ (symmetric positive definite, easy to invert), solve the equivalent system: $M^{-1/2} A M^{-1/2} \hat x = M^{-1/2} b, \quad \hat x = M^{1/2} x.$ The effective condition number is $\kappa_2(M^{-1}A)$ . If $M = A$ , $\kappa = 1$ and CG converges in one step.

Implementation. In practice, PCG avoids computing $M^{-1/2}$ explicitly. The residual preconditioned step replaces $r_k$ with $z_k = M^{-1} r_k$ , and the inner products $r_k^\top r_k$ become $r_k^\top z_k$ .

Common preconditioners:

Preconditioner	Construction cost	Application cost	Best for
Diagonal (Jacobi)	$O(n)$	$O(n)$	Diagonally dominant $A$
Incomplete Cholesky (IC0)	$O(\text{nnz})$	$O(\text{nnz})$	SPD sparse matrices
SSOR	$O(1)$	$O(\text{nnz})$	Banded / finite difference
AMG (algebraic multigrid)	$O(n \log n)$	$O(n)$	Elliptic PDE systems

REMARK

Remark. For yield curve calibration, the normal equations $(A^\top A + \lambda I)x = A^\top b$ with Tikhonov regularisation have an effective condition number $\kappa_2(A^\top A + \lambda I) \approx \sigma_1^2 / \lambda$ when $\lambda$ dominates the small singular values. Diagonal preconditioning $M = \lambda I$ (trivially $M^{-1} = I/\lambda$ ) reduces $\kappa$ by the regularisation. This is one way to see why regularisation aids iterative convergence.

4. Thomas Algorithm: Direct O(n) for Tridiagonal Systems

Before covering LSQR, it is important to note the most practically used "iterative-like" solver in quant finance: the Thomas algorithm (tridiagonal matrix algorithm), which is in fact a direct solver.

DEFINITION

Definition 5.4 (Thomas algorithm). For the tridiagonal system $a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i, \quad i = 1, \ldots, n,$ the Thomas algorithm performs Gaussian elimination in $O(n)$ flops via forward sweep (eliminating $a_i$ ) and backward substitution.

Stability. The Thomas algorithm is numerically stable when the matrix is diagonally dominant ( $|b_i| > |a_i| + |c_i|$ ) — which holds for all standard Crank-Nicolson and Crank-Nicolson-like FD time-stepping schemes. No pivoting is needed.

Cost: $O(8n)$ flops — orders of magnitude cheaper than LU ( $O(n^3/3)$ ) for 1D grids. For a grid with $n = 10^4$ nodes, Thomas costs $8 \times 10^4$ flops vs LU's $3 \times 10^{11}$ — a factor of $4 \times 10^6$ speedup.

5. LSQR: Iterative Least Squares

For large overdetermined LS problems $\min \|Ax - b\|_2$ where $A \in \mathbb{R}^{m \times n}$ is sparse, forming $A^\top A$ explicitly destroys sparsity and squares the condition number (Module 4). LSQR (Paige-Saunders, 1982) solves the LS problem iteratively using only matrix-vector products with $A$ and $A^\top$ .

THEOREM

Theorem 5.4 (LSQR). LSQR applies the Lanczos bidiagonalisation to build a sequence of iterates $x_k$ minimising $\|Ax - b\|_2$ over $\mathcal{K}_k(A^\top A, A^\top b)$ . It satisfies:

$x_k$ minimises the LS residual over the $k$ -dimensional Krylov subspace.
Convergence rate: $O\!\left(\left(\frac{\sqrt{\kappa_2(A)}-1}{\sqrt{\kappa_2(A)}+1}\right)^k\right)$ — same dependence as CG on $\sqrt{\kappa_2(A)}$ (not $\kappa_2(A)^2$ as for normal equations via CG).
Numerically equivalent to CG applied to $A^\top A x = A^\top b$ , but avoids forming $A^\top A$ .

When to use LSQR vs direct methods:

$n < 10^3$ : use direct (SciPy lstsq via LAPACK DGELSD).
$n \geq 10^3$ , $A$ sparse: use LSQR (SciPy lsqr).
$n \geq 10^4$ : use LSQR with a good preconditioner on $A^\top A$ .

6. Stopping Criteria

WARNING

Warning. Iterative solvers require explicit stopping criteria. The standard choices are:

Absolute residual: stop when $\|r_k\|_2 \leq \varepsilon_\text{abs}$ .
Relative residual: stop when $\|r_k\|_2 / \|b\|_2 \leq \varepsilon_\text{rel}$ .
Maximum iterations: always set a cap to avoid infinite loops on divergent problems.

For calibration, set $\varepsilon_\text{rel}$ to match quote precision (e.g., $10^{-4}$ for bid-ask), not machine precision ( $10^{-16}$ ). Over-solving wastes time and does not improve the hedge.

Validation

The companion notebook verifies:

Steepest descent vs CG convergence curves for a $50 \times 50$ SPD matrix with known $\kappa$ .
CG terminates in $\leq n$ steps on an $n \times n$ SPD system in exact arithmetic (within floating-point tolerance).
A-conjugacy: $p_i^\top A p_j = 0$ for all $i \neq j$ in the CG sequence.
Preconditioning: diagonal preconditioner reduces iteration count on a poorly scaled matrix.
Thomas algorithm: $O(n)$ tridiagonal solve matches scipy.linalg.solve_banded to machine precision.
LSQR convergence: matches the direct LS solution on a sparse overdetermined system.

PRACTICE

Hand exercise. Consider the $2 \times 2$ SPD system $A = \begin{pmatrix} 4 & 1 \\ 1 & 3 \end{pmatrix}$ , $b = \begin{pmatrix} 1 \\ 2 \end{pmatrix}$ , starting from $x_0 = (0, 0)^\top$ .

(a) Compute $r_0 = b - Ax_0$ . What is $p_0$ ?

(b) Compute $\alpha_0 = r_0^\top r_0 / (r_0^\top A r_0)$ .

(d) For a $2 \times 2$ SPD system, CG converges in at most 2 steps. Is $r_1 \perp r_0$ ?

(Answer: $r_0 = (1,2)^\top$ ; $r_0^\top r_0 = 5$ ; $Ar_0 = (6, 7)^\top$ ; $r_0^\top Ar_0 = 20$ ; $\alpha_0 = 1/4$ ; $x_1 = (1/4, 1/2)^\top$ ; $r_1 = (1 - 6/4, 2 - 7/4)^\top = (-1/2, 1/4)^\top$ ; $r_1 \cdot r_0 = -1/2 + 1/2 = 0$ ✓)

Limitations

WARNING

Warning: CG requires $A$ to be SPD. If $A$ is only SPSD (rank-deficient), CG breaks down when the search direction $p_k$ becomes $A$ -singular ( $p_k^\top A p_k = 0$ ). For indefinite or non-symmetric $A$ , use MINRES (symmetric indefinite) or GMRES (non-symmetric). GMRES stores all previous residuals and may require restarts for large problems.

WARNING

Warning: floating-point loss of orthogonality. In finite-precision arithmetic, the theoretical A-conjugacy $p_i^\top Ap_j = 0$ degrades as $k$ grows. For ill-conditioned problems, "ghost" eigenvalues (Ritz values that have already converged) pollute the Krylov subspace, causing CG to stagnate. Remedies: periodic restart, re-orthogonalisation (Lanczos with full re-ortho), or switching to a robust direct solver once the problem is small enough.

WARNING

Warning: Thomas algorithm fails with exact zeros on the diagonal. If any pivot $b_i - a_i c_{i-1} / b_{i-1}$ becomes zero during the forward sweep, the algorithm fails. For FD schemes from parabolic PDEs with positive coefficients, this cannot happen — but for schemes with non-standard boundary conditions, check diagonal dominance before applying Thomas.

Practical failure modes:

Slow convergence: $\kappa$ too large without a preconditioner. Symptom: residual stagnates after initial rapid drop. Fix: apply IC0 or block Jacobi preconditioning.
Divergence in GMRES: restart frequency too low, or the matrix is nearly singular. Increase restart frequency or switch to a direct solve on the ill-conditioned subproblem.
LSQR damping: using LSQR with a damping parameter $\lambda > 0$ solves the Tikhonov-regularised LS — equivalent to the direct Tikhonov formula but without forming $A^\top A$ .

Interview Angle

PRACTICE

L1 (Junior Quant Developer).

Expected depth: algorithm description, when to use which solver.

"What is the conjugate gradient method and when would you use it over LU factorisation?" Answer: CG is an iterative solver for SPD $Ax = b$ that only requires matrix-vector products — cost $O(n \cdot k)$ vs $O(n^3)$ for LU. Use CG when $A$ is large and sparse (e.g., FD stiffness matrix, sparse covariance). Use LU when $A$ is small ( $n < 10^3$ ) or dense.
"Why does CG converge faster than steepest descent?" Answer: Convergence of steepest descent depends on $(\kappa-1)/(\kappa+1)$ ; CG depends on $(\sqrt\kappa-1)/(\sqrt\kappa+1)$ . CG eliminates redundant search directions by enforcing A-conjugacy, so it spans a new dimension of the Krylov subspace at each step. For $\kappa = 10^4$ : steepest descent needs $\sim 10^4$ iterations; CG needs $\sim 100$ .
"What does the Thomas algorithm do, and why is it used in FD pricing?" Answer: The Thomas algorithm is a direct $O(n)$ solver for tridiagonal systems. Crank-Nicolson PDE pricing produces one tridiagonal system per time step. Thomas is the standard solver — far faster than general LU ( $O(n^3)$ ) and exploits the exact sparsity pattern.

PRACTICE

L2 (Senior Quant / Quant Developer).

Expected depth: derivation of CG update, convergence proof sketch, preconditioning design.

"Derive the CG step size $\alpha_k$ from the A-conjugacy requirements." Answer: $x_{k+1} = x_k + \alpha_k p_k$ minimises $f(x_k + \alpha p_k)$ over $\alpha$ . Differentiating: $\frac{d}{d\alpha}f = p_k^\top(A(x_k + \alpha p_k) - b) = -p_k^\top r_k + \alpha p_k^\top A p_k = 0$ . So $\alpha_k = p_k^\top r_k / p_k^\top A p_k$ . By construction $p_0 = r_0$ , and subsequently $p_k^\top r_k = r_k^\top r_k$ (provable by induction from the Gram-Schmidt A-conjugation step).
"A Heston PDE on a $100 \times 100$ (spot × vol) grid gives a $10^4 \times 10^4$ sparse system at each time step. What solver strategy would you recommend?" Answer: ADI (alternating direction implicit) splits the 2D problem into two 1D tridiagonal solves per time step — each solved with the Thomas algorithm in $O(n)$ . If ADI is not accurate enough (e.g., near expiry with large mixed derivative terms), use PCG with IC0 preconditioning. Target residual tolerance $\sim 10^{-8}$ (matches pricing accuracy requirement of $10^{-8}$ relative price error).

PRACTICE

L3 (Quant Researcher / Model Risk).

Expected depth: Krylov subspace theory, connection to polynomials and eigenvalues, early stopping as regularisation.

"Explain why CG convergence depends on the distribution of eigenvalues of $A$ , not just on $\kappa_2(A)$ ." Answer: CG minimises $\|e_k\|_A$ over the Krylov subspace $\mathcal{K}_k$ . The error satisfies $\|e_k\|_A \leq \min_{p_k \in \mathcal{P}_k, p_k(0)=1} \max_{\lambda \in \sigma(A)} |p_k(\lambda)| \|e_0\|_A$ , where $\mathcal{P}_k$ is the set of degree- $k$ polynomials. If eigenvalues cluster in a few groups (e.g., $k$ distinct eigenvalues), CG converges in exactly $k$ steps regardless of $\kappa_2$ . The worst case is uniformly spread eigenvalues (where Chebyshev polynomials give the $\sqrt\kappa$ bound). For FD matrices, eigenvalues cluster near two values at fine grids — explaining why PCG with a good preconditioner converges in far fewer steps than the $\sqrt\kappa$ bound predicts.
"In what sense is early stopping of an iterative solver equivalent to regularisation, and how would you set the tolerance for a calibration problem?" Answer: CG at iteration $k$ produces the minimiser of $f$ over $\mathcal{K}_k$ . This is equivalent to applying the degree- $k$ polynomial filter $p_k(\lambda)$ to the eigenvalues of $A^{-1}$ . Small eigenvalues of $A$ (corresponding to large eigenvalues of $A^{-1}$ , i.e., ill-conditioned directions) are filtered out in early iterations — CG only resolves them once $k$ is large enough that $\mathcal{K}_k$ spans those directions. Stopping early thus suppresses contributions from the null-space-proximate directions, exactly like TSVD regularisation (Module 4). The optimal tolerance is $\varepsilon_\text{rel} \approx \|\delta b\|/\|b\|$ where $\delta b$ is the data noise level — solving beyond this adds computational cost without improving the solution quality beyond what the data supports.