The Solver That Knows What It Doesn't Know

Linear algebra is inference

The tutorial's opening move (Hennig, Pfoertner & Weiland, Probabilistic Numerics, ICML Seoul 2026, pp. 7–8) is to refuse the usual staging. Solving $Ax_\ast=b$ is not a mechanical act that either finishes or fails; it is inference about an unknown quantity from data you choose to collect. Put a Gaussian prior on the solution,

$p(x)=\mathcal N(x;\,\mu_0,\Sigma_0),\qquad \mu_0=0\ \text{w.l.o.g.},$

and note what a solver actually gets to see. Never $x$. It picks an action $s_i\in\mathbb R^D$, multiplies once by $A$, and observes one projection of the right-hand side — a noiseless linear functional of the unknown:

$p(\bar b_i\mid x,A,s_i)\;=\;\delta\!\left(s_i^\top b-s_i^\top\!A\,x\right)\;=\;\delta\!\left(\bar b_i-y_i^\top x\right),\qquad y_i:=A^\top\! s_i.$

That is a regression observation with a delta likelihood — Bayes for Gaussians, in the noiseless limit. Condition on $D$ independent projections and the posterior collapses to $\delta(x-A^{-1}b)$: exact solving is the boring endpoint of inference. Everything interesting happens before the collapse, when the solver has spent $n<D$ matrix–vector products and must say what it believes about the rest.

Loading b one projection at a time

Because each observation is one linear functional, the posterior update is a rank-one surgery on the Gaussian (slides pp. 9–12). With $v_i$ the "new direction learned" and $\xi_i$ "how much of $b$ was found along it":

$v_i=\dfrac{\Sigma_{i-1}y_i}{\sqrt{y_i^\top\Sigma_{i-1}y_i}},\qquad \xi_i=\dfrac{\bar b_i-y_i^\top\mu_{i-1}}{\sqrt{y_i^\top\Sigma_{i-1}y_i}},\qquad \mu_i=\mu_{i-1}+\xi_i v_i,\qquad \Sigma_i=\Sigma_{i-1}-v_iv_i^\top.$

Run the first step by hand on the arena this page borrows from Shewchuk — $A=\begin{pmatrix}3&2\\2&6\end{pmatrix}$, $b=(2,-8)$, $x_\ast=(2,-2)$, eigenvalues $\{2,7\}$ — with the prior $\Sigma_0=A^{-1}=\frac{1}{14}\begin{pmatrix}6&-2\\-2&3\end{pmatrix}$ and the residual action $s_1=r_0=b=(2,-8)$. Everything stays rational:

$y_1=As_1=(-10,-44),\quad \Sigma_0y_1=s_1=(2,-8),\quad y_1^\top\Sigma_0y_1=332,$
$\bar b_1=s_1^\top b=68,\quad \mu_1=\tfrac{68}{332}\,(2,-8)=\left(\tfrac{34}{83},\,-\tfrac{136}{83}\right),\quad \Sigma_1=\tfrac{1}{1162}\begin{pmatrix}484&-110\\-110&25\end{pmatrix}.$

(With this prior, $\Sigma_0y_1=A^{-1}\!As_1$ is just $s_1$ again — file that away for §4.) The update looks humble, but the slides' Proposition (pp. 13–14) says it is secretly a classical algorithm: the $v_i$ come out orthonormal in the Mahalanobis inner product $\langle u,v\rangle_{\Sigma_0}=u^\top\Sigma_0^{-1}v$, so iterative conditioning implements Gram–Schmidt on the observed directions $\{y_i\}$ in $(\mathbb R^D,\langle\cdot,\cdot\rangle_{\Sigma_0})$, and the posterior is just the running tally in that frame:

$\Sigma_i=\Sigma_0-\textstyle\sum_{j\le i}v_jv_j^\top,\qquad \mu_i=\textstyle\sum_{j\le i}\xi_j v_j.$

Stop after any $i$ and you hold a posterior measure, not a point. This is the same move report 09 §5 made from the other shore — CG as sequential regression on precision-uncorrelated directions — now stated as a theorem about the conditioning process itself, before any policy is chosen.

Computation is active learning

The recursion accepts any actions. Choosing them is experimental design, and the slides name the design space's two ends — both of which this suite has already built, without knowing it was standing in the same room.

Corner case 1 — stability. Feed the solver the coordinate actions $s_n=e_n$ in order, keep the prior $\Sigma_0=A^{-1}$, and the bookkeeping matrix $R_{ij}=v_i^\top y_j$ comes out upper triangular — measured deviation 7.9×10⁻¹⁵ on the chain — with $R^\top R=A$ to 4.1×10⁻¹⁵ and $R^\top=\mathrm{chol}(A)$ including signs. The frame matrix $V$ is upper triangular too (1.3×10⁻¹⁶) with $VV^\top=A^{-1}$ (2.45×10⁻¹⁴): conditioning on coordinates is Cholesky, one pivot per observation [Hennig, SIOPT 2015]. The suite met this exact object as Shawe–Taylor's Algorithm 5.12 in the update-views explainer (§9–10), where partial Gram–Schmidt = incomplete Cholesky; here it returns as the stable corner of a probabilistic design space.

Corner case 2 — greedy information. Let the solver look at its own posterior and demand the observation that reduces residual error fastest: $s_n=r_{n-1}=b-A\mu_{n-1}$. For general $A$ this yields conjugate-direction methods on the normal equations [Dennis Jr & Turner 1987]. For SPD $A$ with the prior $\Sigma_0=A^{-1}$ it yields — exactly, iterate for iterate — conjugate gradients [Cockayne, Oates, Ipsen & Girolami 2019, BayesCG]. That is the headline, and the next section executes it.

Greedy actions + prior A⁻¹ = conjugate gradients, executed

Two arenas, one theorem. On the suite's Dirichlet chain — $A=\mathrm{tridiag}(-1,2,-1)/h^2$ at $n=32$, $h=\tfrac1{33}$, with the green-tents right-hand side (+2 W heater at node 7, −1.5 W chiller at node 25, entries $f/h$, so $b_7=66$, $b_{25}=-49.5$) — we ran the house CG (python/pcg.py, no preconditioner) and, separately, the §2 conditioning loop with greedy actions and prior $A^{-1}$. The posterior means match the CG iterates for every $n=0\dots32$ to a relative 2.1×10⁻¹⁴, and the learned frame is the CG geometry itself:

$v_n=\dfrac{p_n}{\sqrt{p_n^\top A\,p_n}}$ — the $A$-normalized search directions, all $n$, to 5.1×10⁻¹⁵.

And Shewchuk's central identity becomes a change of vocabulary, not of substance: the Gram matrix $p_i^\top A\,p_j$ (his conjugacy) and the Gram matrix $\langle v_i,v_j\rangle_{\Sigma_0}$ (the slides' Mahalanobis orthonormality) are the same matrix, entry for entry, to 2.9×10⁻¹⁴ — because with $\Sigma_0=A^{-1}$ the Mahalanobis inner product is the $A$-inner product.

On the 2×2 the entire double-run is exact rational arithmetic. CG's numbers (left) and the probabilistic solver's (right), same rows:

Both columns land on $x_\ast$ in exactly two steps, and $\Sigma_{n-1}y_n$ is parallel to $p_n$ at both steps — the conditioning rule manufactures Shewchuk's conjugate directions without ever being told about conjugacy. As a bonus, Parseval in the Mahalanobis frame: $\xi_1^2+\xi_2^2=\tfrac{1156}{83}+\tfrac{504}{83}=20=\|x_\ast\|_A^2$ — the energy of the solution is recovered coefficient by coefficient.

Why $\Sigma_0=A^{-1}$ is not a trick. Two independent justifications. Statistically, it is the unique prior on $x$ consistent with a Gaussian-process prior on the underlying regression function (slides pp. 21–23; Hennig, SIOPT 2015), and it carries the signature property $S_n^\top A\mu_n=S_n^\top b$: explored directions are residual-free (slides p. 24; Wenger et al., NeurIPS 2022), verified here to 3.6×10⁻¹³ across all three policies of §5. Physically, this suite has already measured what $A^{-1}$ is:

Uncertainty accounting: one unit per observation, any policy

Here is the part classical CG cannot say. Weight the posterior covariance by the operator itself and the bookkeeping becomes exact integer arithmetic: since each update subtracts $v_nv_n^\top$ with $v_n^\top Av_n=1$ (that is just the Mahalanobis normalization),

$\operatorname{tr}(A\Sigma_n)\;=\;32-n\qquad\text{exactly, for every policy.}$

Measured: greedy to 1.6×10⁻¹³; greedy, coordinate and random simultaneously to 1.2×10⁻¹¹. Every observation — clever or stupid — removes exactly one unit of $A$-weighted uncertainty. And the units mean something: $\mathbb E\,\|x-\mu_n\|_A^2=\operatorname{tr}(A\Sigma_n)$ over the prior, confirmed by Monte Carlo with $K=2000$ draws and fixed actions to 0.49% relative:

So if every policy spends the same budget, what separates them? Where the budget lands. The trace averages over all right-hand sides the prior can produce; a policy is good when it spends its 32 units in the directions this particular $b$ actually excites. Three policies, same prior, same conditioning rule:

Watch it happen node by node. The widget below runs the §2 conditioning loop live in your browser — dense 32×32 Gaussian conditioning, all three policies — and cross-checks itself against the JSON at load:

The honest footnote. The dashed band is a prior-average statement, not a per-$b$ guarantee: for this particular right-hand side the greedy error sits far below $\sqrt{32-k}$ once the policy locks onto the two sources, and no contradiction is involved — one draw need not equal the mean. Is the accounting honest as probability, though? Calibration spot-check: draw $x\sim\mathcal N(0,A^{-1})$ 200 times, condition every draw on the same 5 fixed actions, and test the truth against its own error bars. The Mahalanobis statistic $(x-\mu_5)^\top\Sigma_5^{+}(x-\mu_5)$ should be $\chi^2_{27}$ (rank of $\Sigma_5$ = 27, exactly): measured mean 26.82 vs 27, variance 53.42 vs 54, and 90.0% of draws fall under the $\chi^2_{27}$ 90% quantile 36.74. The error bars mean what they claim.

Shewchuk's pictures, given error bars

Every figure in the 1994 tutorial has a posterior reading. His quadratic-form contours (§3) are level sets of $\tfrac12x^\top Ax-b^\top x$ around $x_\ast$; the credible ellipses of $\mathcal N(\mu_n,\Sigma_n)$ have the same axes — eigenvectors $(2,-1)$ and $(1,2)$ of $A$, axis ratio $\sqrt{7/2}\approx1.87$ — because with $\Sigma_0=A^{-1}$ belief-geometry and energy-geometry are inverses of one another. His steepest-descent zigzag (§6) and his conjugate two-step solve (§7–8) become two information policies:

The "amnesiac" line is not a metaphor — it is an identity. Feed the §2 rule the residual action but reset the covariance to the prior first (condition, keep $\mu$, forget $\Sigma$). With $\Sigma_0=A^{-1}$ and $y=Ar$:

$v=\dfrac{\Sigma_0y}{\sqrt{y^\top\Sigma_0y}}=\dfrac{r}{\sqrt{r^\top\!Ar}},\qquad \xi=\dfrac{r^\top r}{\sqrt{r^\top\!Ar}},\qquad \mu'=\mu+\dfrac{r^\top r}{r^\top\!Ar}\,r$

— which is exactly Shewchuk's steepest-descent step. Steepest descent is the probabilistic solver with amnesia: same greedy action, no memory of which directions were already paid for, so it keeps re-buying overlapping information at the full price (verified live in the widget's self-test below). CG's only extra ingredient is the Gram–Schmidt bookkeeping that §2 gets for free from conditioning. Watch both on his own arena:

Where the suite has already been standing