Inverse Problems in Vision

An inverse problem asks for the causes of an observation: the sharp image behind a blurred one, the 3D shape behind a shaded picture, the motion behind two frames, the tissue behind a set of X-ray projections. The direction from causes to observations, the forward model, is the part that physics and optics describe; the way back is harder, because the forward model loses information and the observations carry noise. Computer vision reports that Szeliski (2022, §1.1) calls vision itself an inverse problem, one that recovers unknowns from information insufficient to determine them fully. This article collects what the many inverse problems of vision share: what makes a problem well or ill posed, how ill-posedness shows up after discretization as ill-conditioning, why sampling a problem more finely can make it worse, and the three remedies that recur across the field, regularization, Bayesian priors and learned priors. The task articles own their problems: image restoration owns ill-posedness through the null space and the linear degradation family, singular value decomposition the pseudoinverse and truncation, and image deblurring the regularized filters and the stopping rules.

Definition

An inverse problem is the problem of estimating an unknown $x$ from data $y$ produced by a known forward model $F$ and noise, $y = F(x) + n$. The forward problem computes $y$ from $x$.

In this article’s definition, a problem is well posed if three conditions hold: a solution exists for every admissible datum, the solution is unique, and it depends continuously on the data, so that small changes in $y$ produce small changes in $x$. A problem that violates any of the three is ill posed. The three conditions carry the name of Jacques Hadamard; no source that states them was read for this article, so they are given here as its own definition.

Intuition

The forward models of vision can lose information. A camera projects 3D onto 2D, so depth is lost; Murphy (2022, §2.3.3) calls inferring a 3D shape from a 2D image a classic problem in visual scene understanding and a fundamentally ill-posed one, because several hidden shapes are consistent with the same image. A blur averages neighbouring pixels, so fine detail is attenuated, sometimes to nothing. A sensor with a colour filter array records one colour per pixel, so two thirds of the values are missing. The forward direction is many-to-one, or nearly so, and the way back must choose.

Two different things can go wrong. Information can be lost outright, so that different causes give exactly the same data (uniqueness fails). Or information can be kept but so faintly that noise swamps it (stability fails). A blur whose transfer function has no zeros is invertible on paper, yet image restoration notes that it is ill-conditioned: the data constrain its weakly transmitted frequencies only up to amplified noise. The first failure needs extra information to choose among exact solutions; the second needs it to suppress the noise that inversion would release. Regularization, priors and learned models are all ways of supplying that information.

Formal Definition

Linear inverse problems and the singular value decomposition

Many of the problems below are linear or linearized: $y = Ax + n$ with the unknown image or field written as a vector. Image Restoration lists the operators $A$ that the restoration tasks use: the identity for denoising, convolution for blur, a 0/1 mask for missing pixels, blur followed by subsampling for super-resolution. Murphy (2023, §28.2.4) writes the same family as a linear-Gaussian forward model $p(y \mid x) = \mathcal{N}(y \mid Wx, \sigma^2 I)$ and notes that inverting it is ill-posed because many latent images map to the same observation.

With the singular value decomposition $A = \sum_i \sigma_i u_i v_i^{\mathsf T}$, the three conditions can be read off. A solution exists for every $y$ only if the $u_i$ with $\sigma_i > 0$ span the data space; it is unique only if $A$ has no null space (no zero $\sigma_i$); and the least-squares solution

$$\hat x = \sum_{\sigma_i > 0} \frac{u_i^{\mathsf T} y}{\sigma_i}\, v_i = x + \sum_{\sigma_i > 0} \frac{u_i^{\mathsf T} n}{\sigma_i}\, v_i$$

(the second form when $A$ has no null space) divides each noise coefficient by $\sigma_i$, so its stability is measured by the condition number $\kappa = \sigma_{\max}/\sigma_{\min}$, which linear algebra defines. The SVD article works through the pseudoinverse, its minimum-norm choice in the null space, and truncation of small singular values; Image Restoration works through the split of the error between the null space and its complement.

Ill-posed and ill-conditioned

In this article’s usage, which follows Image Deblurring’s report of Szeliski (2022, §4.2), “ill-posed” describes the problem and “ill-conditioned” its discretized matrix: deblurring is ill-posed because many sharp images explain the same data, and ill-conditioned because small changes in the data produce large changes in the solution. A finite matrix always depends continuously on its data, so the third condition fails only in degree after discretization, through a large $\kappa$. How large $\kappa$ becomes depends on how fast the singular values fall. In our reading, and as the example below shows for two operators, a problem whose singular values fall like a power of the index ($\sigma_i \sim i^{-p}$) has a condition number that grows like a power of the number of unknowns, while one whose singular values fall faster than any power, as a Gaussian blur’s do, becomes numerically singular at a modest resolution. Differentiation of sampled data, the inverse of integration, is of the first kind; scale-space reports Lindeberg’s (2008) view that defining derivatives through Gaussian smoothing makes the ill-posed problem of differentiating sampled, noisy data well-posed.

Regularization in the Tikhonov sense

The classical remedy replaces the inverse by a nearby well-posed problem. Tikhonov regularization, in the form that image deblurring gives after Sage et al. (2017), minimizes

$$\lVert y - Ax \rVert^2 + \lambda \lVert L x \rVert^2 , \qquad \hat x_\lambda = (A^{\mathsf T}A + \lambda L^{\mathsf T}L)^{-1} A^{\mathsf T} y ,$$

with $L = I$ or a difference operator. With $L = I$ it acts on each singular direction by a filter factor: Murphy (2022, Eq. 11.69) shows that ridge regression multiplies the $j$-th component of the least-squares fit by $\sigma_j^2/(\sigma_j^2 + \lambda)$, which is close to 1 where $\sigma_j^2 \gg \lambda$ and close to 0 where $\sigma_j^2 \ll \lambda$. Truncating the SVD replaces the smooth factor with a 0/1 step, and stopping an iteration early, as Image Deblurring derives for Landweber’s iteration, gives another set of filter factors; early stopping compares gradient descent with ridge on a learned deblurring filter. The weight $\lambda$ trades fidelity to the data against stability. Choosing it from the data is its own problem: Image Deblurring describes a discrepancy rule (stop or regularize until the residual is as small as the noise alone would leave it), and total variation denoising applies the discrepancy principle for its weight. The word “regularization” also has a machine-learning sense, any change to training that aims to lower error on new data, which regularization covers; the two meet in ridge regression, which Murphy (2022, §11.3) also calls $\ell_2$ regularization or weight decay, Tikhonov regularization with $L = I$ on a model’s weights.

The Bayesian view

The same estimate has a probabilistic reading. Murphy (2022, §2.3.3) describes inverse probability as inverting the mapping from the state of the world $h$ to observations $y$, and the Bayesian answer as the posterior $p(h \mid y)$, which needs the forward model $p(y \mid h)$ and a prior $p(h)$ that rules out or down-weights implausible states. Murphy (2023, §27.1) puts the vision case in the same terms: when the forward mapping from scene to image is many-to-one (“perceptual aliasing”), the inverse mapping is ill-posed, and a prior is needed to make the estimate well defined; in simple settings a point estimate such as the MAP estimate $\hat z = \arg\max_z \log p(z) + \log p(x \mid z)$ can be used (in Murphy’s notation, $z$ the scene and $x$ the image), and in computer vision this approach is known as vision as inverse graphics or analysis by synthesis. With Gaussian noise the negative log-posterior is a data term plus a regularizer, and a Gaussian prior gives the Tikhonov estimate, as Image Restoration derives; maximum a posteriori estimation covers MAP in general and Bayesian inference the full posterior.

The posterior says more than its peak. MacKay (2003, Exercise 22.14) remarks that in ill-posed problems a typical posterior is often a weighted superposition of Gaussians with different means and widths, so that its density maximum sits near the narrowest component rather than the one with the most weight; in our reading, this is a reason to report posterior means, samples or uncertainty rather than only the MAP estimate when the data are weak.

Properties

Remedies share one structure

In our reading, every method in the sections below adds information that the data lack, in one of three forms:

• A penalty or constraint on the solution: Tikhonov’s $\lVert Lx \rVert^2$, total variation, sparsity in a wavelet basis, non-negativity, or the smoothness terms of variational flow and shape from shading.
• A prior distribution, of which the penalty is the negative logarithm: Gaussian and Markov random field priors (Markov random field), the patch prior below, or the posterior samples that diffusion and other generative models provide.
• Learned or implicit priors: a network trained on pairs of clean and degraded images, or the structure of an untrained network.

Each changes the problem being solved. A regularized solution is the exact answer to a different, well-posed problem, and its error has two parts that the weight $\lambda$ trades against each other: bias from the penalty and noise passed through the inverse.

Learned and implicit priors

Murphy (2023, §28.2.4) describes a patch prior learned from data: Zoran and Weiss’s Gaussian mixture model of image patches, fitted to two million clean patches and used inside the objective $\frac{1}{2\sigma^2}\lVert Wx - y\rVert^2 - \text{EPLL}(x)$, where the expected patch log-likelihood sums the log-prior over overlapping patches, which Murphy notes is not a valid likelihood because it overcounts pixels, and which is optimized by half-quadratic splitting. He then describes two ways to replace the optimization by a network: training an amortized inference network on pairs produced by the optimizer, or training a discriminative network directly on pairs of clean images and synthetic corruptions, which works well but needs a different network for every operator and sometimes for every noise level (§28.2.4.2). Image Restoration and Image Deblurring describe half-quadratic splitting with a learned denoiser in place of the prior step (Zhang et al., 2017), which keeps $A$ explicit and learns only the prior.

The Deep Image Prior of Ulyanov, Vedaldi and Lempitsky (CVPR 2018) shows that a prior need not be learned at all. They write denoising, super-resolution and inpainting as $x^* = \min_x E(x; x_0) + R(x)$ (their Eq. 1) and replace the regularizer by the parametrization of a randomly initialized convolutional generator, fitting its weights to the single degraded image (Eq. 2); the implied $R(x)$ is 0 for images the network can produce and $+\infty$ otherwise. They report that the network fits natural images much faster than noise, “high impedance to noise and low impedance to signal” in their words, and restrict the number of iterations accordingly; in our reading, early stopping again acts as the regularizer. No aspect of the network is learned from data.

Blind problems

When the forward model itself is unknown, the problem is blind. Murphy (2023, §28.2.4.3) describes blind inverse problems as much harder, and an EM-like alternation between estimating the image and the operator. Image Deblurring covers blind deconvolution, where the scale trade between kernel and image is one of the degeneracies.

Finer sampling and the continuous problem

A discretized problem inherits the singular values of its continuous operator, plus more of them at each finer sampling. For an operator whose singular values fall quickly, sampling more finely adds directions that the data cannot constrain, so the condition number of the matrix grows with resolution while the information in the data does not. In our reading, this is why an unregularized method can get worse as the sampling becomes finer, and why a regularized one, which suppresses those directions, need not. The example below measures both effects.

Examples

Inverse problems across vision

Problem	Forward model	What fails	Where it is treated
Denoising	identity plus noise	none of the three in our reading ($\kappa = 1$); a prior is needed only to remove the noise	image denoising, Wiener filter, BM3D
Deblurring	convolution with a PSF	stability; uniqueness where the transfer function has zeros	image deblurring
Super-resolution	blur, then subsampling	uniqueness: the super-resolution article states that many HR images produce the same LR image	super-resolution
Demosaicing	colour filter array sampling	uniqueness: Ehret and Facciolo (2019), as Demosaicing reports, note the mosaic has a third as many values as the image	demosaicing
Inpainting	a 0/1 mask	uniqueness on the missing pixels	image inpainting
Dehazing	the haze model	uniqueness: $3N$ equations for $4N + 3$ unknowns, as Image Dehazing reports from He et al.	image dehazing
Matting	the compositing equation	uniqueness: three measurements and seven unknowns per colour pixel, as Image Matting reports	image matting
Optical flow	brightness constancy	uniqueness: one equation for two unknowns per pixel	brightness constancy, aperture problem, Horn–Schunck
Depth from one image	perspective projection	uniqueness: a scene scaled about the camera centre gives the same image	monocular depth estimation
Shape from shading	the reflectance map	contested: see below	shape from shading
CT and MRI	projections; Fourier samples	stability, and uniqueness when undersampled	medical image reconstruction, computed tomography, MRI

Some entries need care. For shape from shading, the article on it reports that Horn and Brooks (1989) argue the problem is ill-posed for an image patch without singular points or limbs but not necessarily when they are included, and that Prados and Faugeras (2005) showed that a perspective model with the light at the centre of projection becomes well-posed when the $1/r^2$ fall-off of the light is kept; it is not simply “ill-posed”. Photometric stereo removes the ambiguity by adding images under other lights until the per-pixel estimate is well posed, and the matting article reports that shooting a still object against two known backings gives a closed-form solution (Smith and Blinn, 1996). Adding measurements is the remedy that changes the problem itself rather than its solution. For medical imaging, Medical Image Reconstruction reports Adler and Öktem’s observation that for typical forward operators small changes in the data produce large changes in the reconstruction and that the maximum-likelihood solution over-fits the noise.

Finer sampling against a fixed operator

The script below fixes two continuous operators on the interval $[0, 1]$ and samples them ever more finely, N = 64, 128, 256 and 512 samples. It measures how conditioning and error change with N for an unregularized solver and for Tikhonov regularization at the best weight. Linear algebra shows the condition number of one $64 \times 64$ blur and what quantization does to its exact solve; here the blur’s width is fixed in units of $x$, so its standard deviation is 1.28 samples at N = 64 and 10.24 at N = 512, and the question is how the problem changes with the sampling.

The set-up, in the code’s conventions:

• Ground truth. A scanline $f$ equal to 0.2, except 0.8 on $[0.15, 0.35)$, a ramp from 0.2 to 0.8 on $[0.5, 0.8)$ and 0.5 on $[0.85, 0.95)$, sampled at the cell centres $x_j = (j + 0.5)/N$.
• Operators. A Gaussian blur of standard deviation 0.02 (units of $x$), written by the midpoint rule as $A_{ij} = h\, g(x_i - x_j)$ with $h = 1/N$, $f$ taken as zero outside $[0, 1]$ and rows not renormalized; integration, $A = h \cdot$ lower-triangular ones, the integral of $f$ from 0 to the right edge of each cell; and a control, $0.5 I + 0.5 \times$ the blur, which is well conditioned.
• Data. $y = Ax + n$ with the same matrix, so that the noise, Gaussian with standard deviation 0.01 per sample, is the only error in the data; seeds 0–9, with generator default_rng(1000 * seed + N).
• Solvers. NumPy’s lstsq with its default cutoff and with rcond=1e-15. Tikhonov in continuous units, $\lVert Ax - y \rVert^2 + \lambda \lVert Lx \rVert^2$ with $L = I$ (by SVD filter factors) or $L = D/h$, the forward difference as a derivative (by QR of the stacked matrix $[A; \sqrt{\lambda} L]$); $\lambda$ on a grid from $10^{-14}$ to $10^{2}$ in steps of 0.1 decade, with the oracle $\lambda$, the one with the smallest error against the known truth, chosen for each seed. The oracle is available only because the truth is known; it shows what the best fixed weight can do, not how to choose one.
• Error. The relative error $\lVert \hat x - x \rVert / \lVert x \rVert$.

NumPy 2.3.5’s source states that lstsq treats singular values smaller than rcond times the largest as zero, with a default rcond of machine precision times $\max(M, N)$, while pinv uses $10^{-15}$; the default cutoff therefore rises with N. The code prints each operator’s singular values, the ranks lstsq reports under both cutoffs, octave slopes $\log_2(\sigma_{2i}/\sigma_i)$ as a measure of decay (a power law $\sigma_i \sim 1/i$ gives slopes near −1; a Gaussian gives ever steeper ones), the errors, and a table of the data’s coefficients $|u_i^{\mathsf T} y|$ against $\sigma_i$. It also prints data and component checks and a verdict line for each planned hypothesis. Seeds, operators and every bar were fixed before the first run. It runs with NumPy 2.3.5 and SciPy 1.16.3 in about 15 seconds on one CPU thread.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"
import numpy as np
from scipy.linalg import qr, solve_triangular

NS = [64, 128, 256, 512]                 # samples on [0, 1]
SB, SN = 0.02, 0.01                      # blur standard deviation (units of x), noise SD per sample
SEEDS = range(10)
LAMS = 10.0 ** (np.arange(-140, 21) / 10)  # 1e-14 ... 1e2, 161 values
EPS = np.finfo(float).eps
BREAKS = [0.15, 0.35, 0.5, 0.8, 0.85, 0.95]

def f(x):                                # steps and a ramp
    y = np.full_like(x, 0.2)
    y[(x >= 0.15) & (x < 0.35)] = 0.8
    m = (x >= 0.5) & (x < 0.8); y[m] = 0.2 + 0.6 * (x[m] - 0.5) / 0.3
    y[(x >= 0.85) & (x < 0.95)] = 0.5
    return y

def operators(N):                        # the same continuous operators at every N
    h = 1.0 / N; x = (np.arange(N) + 0.5) * h
    d = x[:, None] - x[None, :]
    B = h * np.exp(-0.5 * (d / SB) ** 2) / (np.sqrt(2 * np.pi) * SB)   # Gaussian blur, midpoint rule
    T = h * np.tril(np.ones((N, N)))                                   # integral from 0 to each cell's right edge
    return x, {"blur": B, "integration": T, "control": 0.5 * np.eye(N) + 0.5 * B}

def tik_I(U, s, Vt, Y, lam):             # min ||Ax - y||^2 + lam ||x||^2 by SVD filter factors
    return Vt.T @ ((s / (s ** 2 + lam))[:, None] * (U.T @ Y))
def tik_stack(A, L, Y, mu):              # min ||Ax - y||^2 + mu ||Lx||^2 by QR of [A; sqrt(mu) L]
    Q, R = qr(np.vstack([A, np.sqrt(mu) * L]), mode="economic")
    return solve_triangular(R, Q[:A.shape[0]].T @ Y)
rel = lambda X, x: np.linalg.norm(X - x[:, None], axis=0) / np.linalg.norm(x)

checks, out, S, coef_rms = [], {}, {}, {}
def check(name, ok, val): checks.append((name, bool(ok), val))
pooled0 = []
for N in NS:
    x, A = operators(N); xt = f(x); h = 1.0 / N
    D = (np.eye(N, k=1) - np.eye(N))[:-1] / h            # first difference in units of x
    dist = min(abs(b * N - 0.5 - round(b * N - 0.5)) for b in BREAKS)
    check(f"D1 N={N} breakpoint to nearest cell centre (cells) > 0.05", dist > 0.05, f"{dist:.2f}")
    NOISE = np.column_stack([SN * np.random.default_rng(1000 * s + N).standard_normal(N) for s in SEEDS])
    pooled0.append(NOISE[:, 0])
    B = A["blur"]
    check(f"C2 N={N} blur symmetric, row sums <= 1+1e-12, min eigenvalue >= -1e-12",
          np.array_equal(B, B.T) and B.sum(1).max() <= 1 + 1e-12 and np.linalg.eigvalsh(B).min() >= -1e-12,
          f"{B.sum(1).max() - 1:+.1e}, {np.linalg.eigvalsh(B).min():.1e}")
    for name, M in A.items():
        U, s, Vt = np.linalg.svd(M); S[name, N] = s
        if name == "integration":
            k = np.arange(1, N + 1)
            d1 = np.max(np.abs(s - h / (2 * np.sin((2 * k - 1) * np.pi / (4 * N + 2)))) / s)
            check(f"C1 N={N} integration singular values vs closed form < 1e-10", d1 < 1e-10, f"{d1:.1e}")
        b = M @ xt                                         # noise-free data
        x0 = np.linalg.lstsq(M, b)[0]
        res0 = np.linalg.norm(M @ x0 - b) / np.linalg.norm(b); e0 = rel(x0[:, None], xt)[0]
        check(f"C3 N={N} {name}: noise-free LS residual < 1e-10", res0 < 1e-10, f"{res0:.1e}")
        if name != "blur" or N == 64:
            check(f"C3 N={N} {name}: noise-free LS error < 1e-9", e0 < 1e-9, f"{e0:.1e}")
        else:
            out[name, N, "noise-free LS error"] = e0
        if name != "blur":
            eI = rel(tik_I(U, s, Vt, b[:, None], LAMS[0]), xt)[0]
            eD = rel(tik_stack(M, D, b[:, None], LAMS[0]), xt)[0]
            check(f"C4 N={N} {name}: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D",
                  max(eI, eD) < 1e-6, f"{eI:.1e}, {eD:.1e}")
        Y = b[:, None] + NOISE
        xd, _, rk_d, _ = np.linalg.lstsq(M, Y)                 # default rcond = eps * max(M, N)
        xf, _, rk_f, _ = np.linalg.lstsq(M, Y, rcond=1e-15)
        out[name, N, "rank"] = (rk_d, rk_f)
        out[name, N, "LS default"] = rel(xd, xt); out[name, N, "LS 1e-15"] = rel(xf, xt)
        c5 = np.max(np.abs(tik_I(U, s, Vt, Y[:, :1], 1e-4) - tik_stack(M, np.eye(N), Y[:, :1], 1e-4))) \
             / np.max(np.abs(tik_I(U, s, Vt, Y[:, :1], 1e-4)))
        check(f"C5 N={N} {name}: Tikhonov L=I by SVD vs by QR < 1e-8", c5 < 1e-8, f"{c5:.1e}")
        EI = np.array([rel(tik_I(U, s, Vt, Y, l), xt) for l in LAMS])         # (lam, seed)
        ED = np.array([rel(tik_stack(M, D, Y, l), xt) for l in LAMS])
        for lab, E in (("Tik L=I", EI), ("Tik L=D", ED)):
            k = E.argmin(0)
            out[name, N, lab] = E[k, np.arange(len(SEEDS))]; out[name, N, lab + " lam"] = LAMS[k]
            out[name, N, lab + " edges"] = int(np.sum((k == 0) | (k == len(LAMS) - 1)))
        C = U.T @ Y
        coef_rms[name, N] = np.sqrt(np.mean(C[N // 2:] ** 2, axis=0)) / SN
        out[name, N, "signal"] = np.sqrt(np.mean((U.T @ b)[N // 2:] ** 2)) / SN
        if N == 256 and name != "control":
            out[name, "table"] = [(i, s[i - 1], abs(C[i - 1, 0])) for i in [2 ** j for j in range(9)]]
sd0 = np.concatenate(pooled0).std()
check("D2 seed 0 noise SD (960 samples) within 0.0007 of 0.01", abs(sd0 - SN) < 0.0007, f"{sd0:.5f}")
for name in ("blur", "integration"):
    for N in (256, 512):
        v = out[name, N, "signal"]
        check(f"D3 N={N} {name}: signal RMS (i > N/2) / sigma_n < 0.05", v < 0.05, f"{v:.1e}")
edges = sum(out[k] for k in out if len(k) == 3 and isinstance(k[2], str) and k[2].endswith("edges"))
check("C6 oracle lam strictly inside the grid in every case", edges == 0, f"{edges} at an edge")

print("Checks")
for name, ok, val in checks:
    print(f"  {'held  ' if ok else 'FAILED'} {name:72s} {val}")

print("\nSingular values (rank = values lstsq keeps under each cutoff)")
print(f"  {'operator':11s} {'N':>4s} {'sigma_1':>10s} {'sigma_N':>10s} {'cond':>10s} {'rank':>9s}   sigma_i at i = 1, 4, 16, 32")
for name in ("blur", "integration", "control"):
    for N in NS:
        s = S[name, N]; rk = out[name, N, "rank"]
        print(f"  {name:11s} {N:4d} {s[0]:10.4e} {s[-1]:10.4e} {s[0] / s[-1]:10.4e} {rk[0]:4d}/{rk[1]:<4d}  "
              + " ".join(f"{s[i - 1]:.6e}" for i in (1, 4, 16, 32)))
c = [S["integration", N][0] / S["integration", N][-1] for N in NS]
print("  integration cond/N:", " ".join(f"{ci / N:.5f}" for ci, N in zip(c, NS)),
      "  cond(N)/cond(N/2):", " ".join(f"{c[j] / c[j - 1]:.5f}" for j in range(1, 4)))
for name in ("blur", "integration"):
    print(f"  {name}: change of sigma_i from N=256 to 512 (%):",
          " ".join(f"{100 * abs(S[name, 512][i - 1] / S[name, 256][i - 1] - 1):.4f}" for i in (1, 4, 16, 32)))

print("\nOctave slopes log2(sigma_2i / sigma_i), i = 1, 2, 4, ..., N/2 (* = sigma_2i < 1e-12 sigma_1)")
for name in ("blur", "integration"):
    for N in NS:
        s = S[name, N]
        print(f"  {name:11s} {N:4d} " + " ".join(
            f"{np.log2(s[2 * i - 1] / s[i - 1]):.4f}" + ("*" if s[2 * i - 1] < 1e-12 * s[0] else "")
            for i in [2 ** j for j in range(int(np.log2(N)))]))

print(f"\nRelative error ||x - x_true|| / ||x_true||, noise SD {SN}: min / mean / max over seeds 0-9")
for name in ("blur", "integration", "control"):
    for N in NS:
        r = [f"{lab} {out[name, N, lab].min():.4e} / {out[name, N, lab].mean():.4e} / {out[name, N, lab].max():.4e}"
             for lab in ("LS default", "LS 1e-15")]
        t = [f"{lab} {out[name, N, lab].min():.5f} / {out[name, N, lab].mean():.5f} / {out[name, N, lab].max():.5f}"
             f" (lam 1e{np.log10(out[name, N, lab + ' lam'].min()):.1f} to 1e{np.log10(out[name, N, lab + ' lam'].max()):.1f})"
             for lab in ("Tik L=I", "Tik L=D")]
        print(f"  {name:11s} {N:4d}  " + "; ".join(r))
        print(f"  {'':11s} {'':4s}  " + "; ".join(t))
for N in NS[1:]:
    print(f"  blur {N}: noise-free LS error (default cutoff) {out['blur', N, 'noise-free LS error']:.4e}")

print("\nCoefficients of the data, seed 0, N = 256: i, sigma_i, |u_i^T y|, |u_i^T y| / sigma_i")
for name in ("blur", "integration"):
    for i, si, ci in out[name, "table"]:
        print(f"  {name:11s} {i:4d} {si:10.3e} {ci:10.3e} {ci / si:10.3e}")
print("RMS of |u_i^T y| over i > N/2, divided by sigma_n: min / max over seeds")
for name in ("blur", "integration", "control"):
    print(f"  {name:11s} " + "  ".join(f"N={N}: {coef_rms[name, N].min():.4f} / {coef_rms[name, N].max():.4f}" for N in NS))

print("\nHypotheses")
def verdict(lab, ok, txt): print(f"  {'held  ' if ok else 'FAILED'} {lab:4s} {txt}")
cb = [S["blur", N][0] / S["blur", N][-1] for N in NS]
verdict("H1a", cb[0] < 1e8 and min(cb[1:]) > 1e12 and cb[1] > cb[0], "blur: cond(64) < 1e8, cond > 1e12 for N >= 128, cond(128) > cond(64)")
rb = {N: out["blur", N, "rank"][0] for N in NS}
verdict("H1b", rb[512] <= 1.1 * rb[256] and rb[256] < 128 and rb[512] < 256, f"blur rank (default cutoff) 256: {rb[256]}, 512: {rb[512]}")
cc = [S["control", N][0] / S["control", N][-1] for N in NS]
verdict("H1c", all(1.25 <= ci / N <= 1.30 for ci, N in zip(c, NS)) and all(abs(c[j] / c[j - 1] - 2) < 0.02 for j in (1, 2, 3))
        and max(cc) <= 2 + 1e-9, f"integration cond/N and doubling; control max cond {max(cc):.6f}")
ch = max(abs(S[nm, 512][i - 1] / S[nm, 256][i - 1] - 1) for nm in ("blur", "integration") for i in (1, 4, 16, 32))
verdict("H1d", ch < 0.01, f"max change of sigma_i (i = 1, 4, 16, 32) from 256 to 512: {100 * ch:.4f}%")
for name in ("blur", "integration"):
    for lab in ("LS default", "LS 1e-15"):
        E = np.array([out[name, N, lab] for N in NS])
        n_up = int(np.sum(np.all(np.diff(E, axis=0) > 0, axis=0)))
        verdict("H2a", n_up == 10, f"{name}, {lab}: error rises at every step of N in {n_up} of 10 seeds")
m = np.array([out["control", N, "LS default"].mean() for N in NS])
verdict("H2b", np.all(np.abs(m / m[0] - 1) < 0.10), "control: mean LS error / its N=64 value: " + " ".join(f"{v:.4f}" for v in m / m[0]))
bound = [2 * SN / np.sqrt(np.mean(f((np.arange(N) + 0.5) / N) ** 2)) for N in NS]
verdict("H2c", np.all(m <= 1.05 * np.array(bound)), "control: mean LS error vs 2 sigma_n / rms(x): " + " ".join(f"{a:.5f}/{b:.5f}" for a, b in zip(m, bound)))
for name in ("blur", "integration", "control"):
    for lab in ("Tik L=I", "Tik L=D"):
        E = np.array([out[name, N, lab] for N in NS]); mm = E.mean(1)
        verdict("H3a", np.all(mm[1:] <= 1.05 * mm[:-1]), f"{name}, {lab}: mean error ratio to previous N: " + " ".join(f"{v:.4f}" for v in mm[1:] / mm[:-1]))
        verdict("H3b", int(np.sum(E[3] <= E[0])) == 10, f"{name}, {lab}: error(512) <= error(64) in {int(np.sum(E[3] <= E[0]))} of 10 seeds")
        rmax = (E / E[0]).max(0)
        verdict("H3c", np.all(rmax <= 1.25), f"{name}, {lab}: largest error(N) / error(64) over seeds {rmax.max():.4f}")
ok4a = True
for N in NS:
    s = S["integration", N]
    sl = [np.log2(s[2 * i - 1] / s[i - 1]) for i in [2 ** j for j in range(2, int(np.log2(N)) - 2)]]
    ok4a &= all(-1.15 <= v <= -0.95 for v in sl)
verdict("H4a", ok4a, "integration: octave slopes from i = 4 to N/8 in [-1.15, -0.95] at every N")
ok4b, steep = True, {}
for N in NS:
    s = S["blur", N]
    sl = [np.log2(s[2 * i - 1] / s[i - 1]) for i in [2 ** j for j in range(int(np.log2(N)))] if s[2 * i - 1] > 1e-12 * s[0]]
    ok4b &= bool(np.all(np.diff(sl) < 0)); steep[N] = min(sl)
ok4b &= steep[256] < -5 and steep[512] < -5
verdict("H4b", ok4b, f"blur: slopes steepen octave by octave; steepest qualifying at 256 and 512: {steep[256]:.4f}, {steep[512]:.4f}")
for name in ("blur", "integration"):
    lo = min(coef_rms[name, N].min() for N in (256, 512)); hi = max(coef_rms[name, N].max() for N in (256, 512))
    verdict("H5", 0.8 <= lo and hi <= 1.2, f"{name}: RMS |u_i^T y| (i > N/2) / sigma_n at N = 256, 512 in [{lo:.4f}, {hi:.4f}]")

Output:

Checks
  held   D1 N=64 breakpoint to nearest cell centre (cells) > 0.05                 0.10
  held   C2 N=64 blur symmetric, row sums <= 1+1e-12, min eigenvalue >= -1e-12    +1.8e-14, 6.4e-04
  held   C3 N=64 blur: noise-free LS residual < 1e-10                             8.9e-16
  held   C3 N=64 blur: noise-free LS error < 1e-9                                 1.1e-13
  held   C5 N=64 blur: Tikhonov L=I by SVD vs by QR < 1e-8                        1.4e-14
  held   C1 N=64 integration singular values vs closed form < 1e-10               1.4e-15
  held   C3 N=64 integration: noise-free LS residual < 1e-10                      1.4e-15
  held   C3 N=64 integration: noise-free LS error < 1e-9                          3.2e-14
  held   C4 N=64 integration: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.6e-11, 2.1e-07
  held   C5 N=64 integration: Tikhonov L=I by SVD vs by QR < 1e-8                 1.6e-14
  held   C3 N=64 control: noise-free LS residual < 1e-10                          1.0e-15
  held   C3 N=64 control: noise-free LS error < 1e-9                              1.3e-15
  held   C4 N=64 control: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.1e-14, 5.9e-11
  held   C5 N=64 control: Tikhonov L=I by SVD vs by QR < 1e-8                     3.3e-15
  held   D1 N=128 breakpoint to nearest cell centre (cells) > 0.05                0.10
  held   C2 N=128 blur symmetric, row sums <= 1+1e-12, min eigenvalue >= -1e-12   +2.2e-16, 2.1e-14
  held   C3 N=128 blur: noise-free LS residual < 1e-10                            1.2e-15
  held   C5 N=128 blur: Tikhonov L=I by SVD vs by QR < 1e-8                       3.3e-14
  held   C1 N=128 integration singular values vs closed form < 1e-10              1.4e-15
  held   C3 N=128 integration: noise-free LS residual < 1e-10                     2.3e-15
  held   C3 N=128 integration: noise-free LS error < 1e-9                         6.1e-14
  FAILED C4 N=128 integration: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 4.6e-11, 2.4e-06
  held   C5 N=128 integration: Tikhonov L=I by SVD vs by QR < 1e-8                1.9e-14
  held   C3 N=128 control: noise-free LS residual < 1e-10                         1.4e-15
  held   C3 N=128 control: noise-free LS error < 1e-9                             1.9e-15
  held   C4 N=128 control: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.1e-14, 1.8e-10
  held   C5 N=128 control: Tikhonov L=I by SVD vs by QR < 1e-8                    3.6e-15
  held   D1 N=256 breakpoint to nearest cell centre (cells) > 0.05                0.10
  held   C2 N=256 blur symmetric, row sums <= 1+1e-12, min eigenvalue >= -1e-12   +2.2e-16, -1.9e-16
  held   C3 N=256 blur: noise-free LS residual < 1e-10                            1.3e-15
  held   C5 N=256 blur: Tikhonov L=I by SVD vs by QR < 1e-8                       6.4e-14
  held   C1 N=256 integration singular values vs closed form < 1e-10              3.7e-15
  held   C3 N=256 integration: noise-free LS residual < 1e-10                     1.7e-15
  held   C3 N=256 integration: noise-free LS error < 1e-9                         2.5e-13
  FAILED C4 N=256 integration: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.3e-10, 2.7e-05
  held   C5 N=256 integration: Tikhonov L=I by SVD vs by QR < 1e-8                5.9e-13
  held   C3 N=256 control: noise-free LS residual < 1e-10                         1.4e-15
  held   C3 N=256 control: noise-free LS error < 1e-9                             2.2e-15
  held   C4 N=256 control: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.2e-14, 5.1e-10
  held   C5 N=256 control: Tikhonov L=I by SVD vs by QR < 1e-8                    6.2e-15
  held   D1 N=512 breakpoint to nearest cell centre (cells) > 0.05                0.10
  held   C2 N=512 blur symmetric, row sums <= 1+1e-12, min eigenvalue >= -1e-12   +2.2e-16, -2.4e-16
  held   C3 N=512 blur: noise-free LS residual < 1e-10                            1.5e-15
  held   C5 N=512 blur: Tikhonov L=I by SVD vs by QR < 1e-8                       1.8e-13
  held   C1 N=512 integration singular values vs closed form < 1e-10              6.1e-15
  held   C3 N=512 integration: noise-free LS residual < 1e-10                     3.7e-15
  held   C3 N=512 integration: noise-free LS error < 1e-9                         8.9e-13
  FAILED C4 N=512 integration: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 3.7e-10, 3.0e-04
  held   C5 N=512 integration: Tikhonov L=I by SVD vs by QR < 1e-8                6.8e-13
  held   C3 N=512 control: noise-free LS residual < 1e-10                         1.6e-15
  held   C3 N=512 control: noise-free LS error < 1e-9                             2.4e-15
  held   C4 N=512 control: noise-free Tikhonov (lam 1e-14) error < 1e-6, L=I and L=D 1.1e-14, 1.5e-09
  held   C5 N=512 control: Tikhonov L=I by SVD vs by QR < 1e-8                    1.0e-14
  held   D2 seed 0 noise SD (960 samples) within 0.0007 of 0.01                   0.00988
  held   D3 N=256 blur: signal RMS (i > N/2) / sigma_n < 0.05                     7.9e-14
  held   D3 N=512 blur: signal RMS (i > N/2) / sigma_n < 0.05                     1.2e-14
  held   D3 N=256 integration: signal RMS (i > N/2) / sigma_n < 0.05              9.2e-03
  held   D3 N=512 integration: signal RMS (i > N/2) / sigma_n < 0.05              3.2e-03
  FAILED C6 oracle lam strictly inside the grid in every case                     5 at an edge

Singular values (rank = values lstsq keeps under each cutoff)
  operator       N    sigma_1    sigma_N       cond      rank   sigma_i at i = 1, 4, 16, 32
  blur          64 9.9812e-01 6.3543e-04 1.5708e+03   64/64    9.981246e-01 9.704126e-01 6.184462e-01 1.461284e-01
  blur         128 9.9812e-01 2.1440e-14 4.6555e+13  127/128   9.981187e-01 9.703215e-01 6.177353e-01 1.462755e-01
  blur         256 9.9812e-01 8.8623e-18 1.1263e+17  130/138   9.981172e-01 9.702990e-01 6.175662e-01 1.463233e-01
  blur         512 9.9812e-01 6.5574e-18 1.5221e+17  129/139   9.981169e-01 9.702934e-01 6.175244e-01 1.463357e-01
  integration   64 6.4161e-01 7.8148e-03 8.2102e+01   64/64    6.416092e-01 9.176727e-02 2.119636e-02 1.125600e-02
  integration  128 6.3911e-01 3.9065e-03 1.6360e+02  128/128   6.391105e-01 9.132880e-02 2.074021e-02 1.039964e-02
  integration  256 6.3786e-01 1.9532e-03 3.2658e+02  256/256   6.378642e-01 9.113029e-02 2.060716e-02 1.018788e-02
  integration  512 6.3724e-01 9.7657e-04 6.5253e+02  512/512   6.372417e-01 9.103624e-02 2.056391e-02 1.013068e-02
  control       64 9.9906e-01 5.0032e-01 1.9969e+00   64/64    9.990623e-01 9.852063e-01 8.092231e-01 5.730642e-01
  control      128 9.9906e-01 5.0000e-01 1.9981e+00  128/128   9.990593e-01 9.851607e-01 8.088676e-01 5.731377e-01
  control      256 9.9906e-01 5.0000e-01 1.9981e+00  256/256   9.990586e-01 9.851495e-01 8.087831e-01 5.731616e-01
  control      512 9.9906e-01 5.0000e-01 1.9981e+00  512/512   9.990584e-01 9.851467e-01 8.087622e-01 5.731678e-01
  integration cond/N: 1.28284 1.27813 1.27570 1.27448   cond(N)/cond(N/2): 1.99265 1.99621 1.99808
  blur: change of sigma_i from N=256 to 512 (%): 0.0000 0.0006 0.0068 0.0085
  integration: change of sigma_i from N=256 to 512 (%): 0.0976 0.1032 0.2099 0.5615

Octave slopes log2(sigma_2i / sigma_i), i = 1, 2, 4, ..., N/2 (* = sigma_2i < 1e-12 sigma_1)
  blur          64 -0.0081 -0.0325 -0.1300 -0.5200 -2.0814 -7.8453
  blur         128 -0.0082 -0.0326 -0.1304 -0.5211 -2.0783 -8.2371 -34.3963*
  blur         256 -0.0082 -0.0326 -0.1305 -0.5214 -2.0774 -8.1851 -31.6500* -14.0392*
  blur         512 -0.0082 -0.0326 -0.1305 -0.5214 -2.0772 -8.1730 -31.2750* -11.6929* -3.1680*
  integration   64 -1.5847 -1.2210 -1.0933 -1.0209 -0.9131 -0.5264
  integration  128 -1.5849 -1.2220 -1.0980 -1.0407 -0.9959 -0.8993 -0.5132
  integration  256 -1.5849 -1.2223 -1.0991 -1.0456 -1.0163 -0.9838 -0.8925 -0.5066
  integration  512 -1.5850 -1.2224 -1.0994 -1.0469 -1.0214 -1.0045 -0.9779 -0.8891 -0.5033

Relative error ||x - x_true|| / ||x_true||, noise SD 0.01: min / mean / max over seeds 0-9
  blur          64  LS default 4.4490e+00 / 6.9866e+00 / 8.6714e+00; LS 1e-15 4.4490e+00 / 6.9866e+00 / 8.6714e+00
                    Tik L=I 0.10619 / 0.11181 / 0.11923 (lam 1e-2.3 to 1e-1.8); Tik L=D 0.10586 / 0.10993 / 0.11472 (lam 1e-6.2 to 1e-5.5)
  blur         128  LS default 2.9349e+10 / 5.7770e+10 / 9.5751e+10; LS 1e-15 3.6284e+10 / 1.0053e+11 / 1.5142e+11
                    Tik L=I 0.11706 / 0.12050 / 0.12628 (lam 1e-2.3 to 1e-1.8); Tik L=D 0.11606 / 0.11889 / 0.12297 (lam 1e-6.2 to 1e-5.7)
  blur         256  LS default 8.2239e+09 / 2.2078e+10 / 3.6344e+10; LS 1e-15 6.2802e+11 / 1.3257e+12 / 2.4534e+12
                    Tik L=I 0.11285 / 0.11503 / 0.11928 (lam 1e-2.6 to 1e-2.0); Tik L=D 0.11290 / 0.11479 / 0.11811 (lam 1e-6.9 to 1e-6.0)
  blur         512  LS default 3.7855e+09 / 8.9979e+09 / 1.4636e+10; LS 1e-15 2.1934e+11 / 6.5345e+11 / 2.0173e+12
                    Tik L=I 0.11214 / 0.11484 / 0.11877 (lam 1e-2.9 to 1e-2.3); Tik L=D 0.11216 / 0.11488 / 0.11843 (lam 1e-7.2 to 1e-6.5)
  integration   64  LS default 1.2687e+00 / 1.9122e+00 / 2.1431e+00; LS 1e-15 1.2687e+00 / 1.9122e+00 / 2.1431e+00
                    Tik L=I 0.26430 / 0.28934 / 0.31990 (lam 1e-2.8 to 1e-2.6); Tik L=D 0.23570 / 0.26083 / 0.29828 (lam 1e-6.1 to 1e-5.4)
  integration  128  LS default 3.1351e+00 / 3.5522e+00 / 4.0216e+00; LS 1e-15 3.1351e+00 / 3.5522e+00 / 4.0216e+00
                    Tik L=I 0.22633 / 0.24959 / 0.27013 (lam 1e-3.0 to 1e-2.8); Tik L=D 0.19804 / 0.21921 / 0.24929 (lam 1e-6.7 to 1e-6.2)
  integration  256  LS default 6.3748e+00 / 7.0884e+00 / 7.5416e+00; LS 1e-15 6.3748e+00 / 7.0884e+00 / 7.5416e+00
                    Tik L=I 0.20422 / 0.22920 / 0.24137 (lam 1e-3.1 to 1e-3.0); Tik L=D 0.18990 / 0.21664 / 0.24066 (lam 1e-6.8 to 1e-6.4)
  integration  512  LS default 1.3644e+01 / 1.4237e+01 / 1.4865e+01; LS 1e-15 1.3644e+01 / 1.4237e+01 / 1.4865e+01
                    Tik L=I 0.20088 / 0.22013 / 0.23953 (lam 1e-3.2 to 1e-3.1); Tik L=D 0.17701 / 0.19724 / 0.23076 (lam 1e-7.1 to 1e-6.6)
  control       64  LS default 2.6115e-02 / 3.5269e-02 / 4.0066e-02; LS 1e-15 2.6115e-02 / 3.5269e-02 / 4.0066e-02
                    Tik L=I 0.02611 / 0.03516 / 0.04007 (lam 1e-14.0 to 1e-2.2); Tik L=D 0.02536 / 0.03479 / 0.03999 (lam 1e-14.0 to 1e-5.7)
  control      128  LS default 3.3771e-02 / 3.6815e-02 / 4.0464e-02; LS 1e-15 3.3771e-02 / 3.6815e-02 / 4.0464e-02
                    Tik L=I 0.03363 / 0.03659 / 0.04038 (lam 1e-14.0 to 1e-2.3); Tik L=D 0.03241 / 0.03547 / 0.03774 (lam 1e-6.8 to 1e-6.1)
  control      256  LS default 3.5970e-02 / 3.7810e-02 / 3.8967e-02; LS 1e-15 3.5970e-02 / 3.7810e-02 / 3.8967e-02
                    Tik L=I 0.03591 / 0.03746 / 0.03848 (lam 1e-2.8 to 1e-2.2); Tik L=D 0.03416 / 0.03531 / 0.03659 (lam 1e-6.7 to 1e-6.5)
  control      512  LS default 3.7121e-02 / 3.8555e-02 / 3.9637e-02; LS 1e-15 3.7121e-02 / 3.8555e-02 / 3.9637e-02
                    Tik L=I 0.03669 / 0.03829 / 0.03921 (lam 1e-2.6 to 1e-2.3); Tik L=D 0.03255 / 0.03387 / 0.03478 (lam 1e-6.9 to 1e-6.8)
  blur 128: noise-free LS error (default cutoff) 1.4453e-03
  blur 256: noise-free LS error (default cutoff) 5.6749e-02
  blur 512: noise-free LS error (default cutoff) 6.2606e-02

Coefficients of the data, seed 0, N = 256: i, sigma_i, |u_i^T y|, |u_i^T y| / sigma_i
  blur           1  9.981e-01  6.664e+00  6.676e+00
  blur           2  9.925e-01  4.133e-01  4.164e-01
  blur           4  9.703e-01  1.012e-02  1.043e-02
  blur           8  8.864e-01  2.497e-01  2.817e-01
  blur          16  6.176e-01  1.152e-01  1.865e-01
  blur          32  1.463e-01  3.876e-02  2.649e-01
  blur          64  5.028e-04  1.319e-02  2.623e+01
  blur         128  1.492e-13  2.584e-03  1.732e+10
  blur         256  8.862e-18  6.521e-03  7.358e+14
  integration    1  6.379e-01  4.115e+00  6.451e+00
  integration    2  2.126e-01  4.524e-01  2.128e+00
  integration    4  9.113e-02  2.705e-01  2.968e+00
  integration    8  4.254e-02  7.349e-03  1.728e-01
  integration   16  2.061e-02  1.281e-02  6.214e-01
  integration   32  1.019e-02  7.810e-03  7.666e-01
  integration   64  5.151e-03  8.016e-03  1.556e+00
  integration  128  2.775e-03  1.583e-02  5.703e+00
  integration  256  1.953e-03  1.889e-02  9.673e+00
RMS of |u_i^T y| over i > N/2, divided by sigma_n: min / max over seeds
  blur        N=64: 0.6665 / 1.2541  N=128: 0.7633 / 1.1333  N=256: 0.8737 / 1.0589  N=512: 0.9372 / 1.0175
  integration N=64: 0.6100 / 1.1846  N=128: 0.7926 / 1.1444  N=256: 0.8805 / 1.0697  N=512: 0.9276 / 1.0475
  control     N=64: 3.6288 / 4.1036  N=128: 2.8250 / 3.2151  N=256: 2.1501 / 2.3543  N=512: 2.0272 / 2.1649

Hypotheses
  held   H1a  blur: cond(64) < 1e8, cond > 1e12 for N >= 128, cond(128) > cond(64)
  FAILED H1b  blur rank (default cutoff) 256: 130, 512: 129
  held   H1c  integration cond/N and doubling; control max cond 1.998119
  held   H1d  max change of sigma_i (i = 1, 4, 16, 32) from 256 to 512: 0.5615%
  FAILED H2a  blur, LS default: error rises at every step of N in 0 of 10 seeds
  FAILED H2a  blur, LS 1e-15: error rises at every step of N in 1 of 10 seeds
  held   H2a  integration, LS default: error rises at every step of N in 10 of 10 seeds
  held   H2a  integration, LS 1e-15: error rises at every step of N in 10 of 10 seeds
  held   H2b  control: mean LS error / its N=64 value: 1.0000 1.0438 1.0720 1.0932
  held   H2c  control: mean LS error vs 2 sigma_n / rms(x): 0.03527/0.04036 0.03681/0.03974 0.03781/0.03962 0.03856/0.03977
  FAILED H3a  blur, Tik L=I: mean error ratio to previous N: 1.0777 0.9546 0.9983
  FAILED H3b  blur, Tik L=I: error(512) <= error(64) in 2 of 10 seeds
  held   H3c  blur, Tik L=I: largest error(N) / error(64) over seeds 1.1892
  FAILED H3a  blur, Tik L=D: mean error ratio to previous N: 1.0816 0.9655 1.0008
  FAILED H3b  blur, Tik L=D: error(512) <= error(64) in 1 of 10 seeds
  held   H3c  blur, Tik L=D: largest error(N) / error(64) over seeds 1.1616
  held   H3a  integration, Tik L=I: mean error ratio to previous N: 0.8626 0.9183 0.9604
  held   H3b  integration, Tik L=I: error(512) <= error(64) in 10 of 10 seeds
  held   H3c  integration, Tik L=I: largest error(N) / error(64) over seeds 1.0000
  held   H3a  integration, Tik L=D: mean error ratio to previous N: 0.8405 0.9883 0.9104
  held   H3b  integration, Tik L=D: error(512) <= error(64) in 10 of 10 seeds
  held   H3c  integration, Tik L=D: largest error(N) / error(64) over seeds 1.0210
  held   H3a  control, Tik L=I: mean error ratio to previous N: 1.0407 1.0238 1.0220
  FAILED H3b  control, Tik L=I: error(512) <= error(64) in 1 of 10 seeds
  FAILED H3c  control, Tik L=I: largest error(N) / error(64) over seeds 1.5020
  held   H3a  control, Tik L=D: mean error ratio to previous N: 1.0196 0.9955 0.9591
  FAILED H3b  control, Tik L=D: error(512) <= error(64) in 7 of 10 seeds
  FAILED H3c  control, Tik L=D: largest error(N) / error(64) over seeds 1.4838
  held   H4a  integration: octave slopes from i = 4 to N/8 in [-1.15, -0.95] at every N
  held   H4b  blur: slopes steepen octave by octave; steepest qualifying at 256 and 512: -8.1851, -8.1730
  held   H5   blur: RMS |u_i^T y| (i > N/2) / sigma_n at N = 256, 512 in [0.8737, 1.0589]
  held   H5   integration: RMS |u_i^T y| (i > N/2) / sigma_n at N = 256, 512 in [0.8805, 1.0697]

Checks. The data checks held: no breakpoint of $f$ lies closer than 0.10 cells to a cell centre, the seed-0 noise has standard deviation 0.00988, and the signal’s part in the upper half of the singular directions is negligible for the blur and integration at N = 256 and 512 (at most $9.2 \times 10^{-3}\,\sigma_n$), so that those coefficients are noise. The integration matrix’s singular values agree with the closed form $h / (2 \sin((2k-1)\pi/(4N+2)))$ (our derivation: for the integration matrix $T$, $h^2 (TT^{\mathsf T})^{-1}$ is the tridiagonal matrix with $-1, 2, -1$ and a last diagonal entry of 1, whose eigenvalues are $4\sin^2((2k-1)\pi/(4N+2))$) to $6.1 \times 10^{-15}$ (C1 held), the blur matrix is symmetric with row sums at most $1 + 1.8 \times 10^{-14}$ and no eigenvalue below $-2.4 \times 10^{-16}$ (C2 held), noise-free data are fitted to a relative residual of at most $3.7 \times 10^{-15}$ and recovered to at most $8.9 \times 10^{-13}$ for the control and integration at every N and for the blur at N = 64 (C3 held), and the two Tikhonov solvers agree to $6.8 \times 10^{-13}$ (C5 held), which is the only evidence the example gives that the stacked-QR solver is right. Two checks failed:

• C4 failed for integration with $L = D/h$: the noise-free error at $\lambda = 10^{-14}$ was $2.4 \times 10^{-6}$, $2.7 \times 10^{-5}$ and $3.0 \times 10^{-4}$ at N = 128, 256 and 512, against a bar of $10^{-6}$. It held for integration with $L = I$, for integration with $L = D/h$ at N = 64 ($2.1 \times 10^{-7}$) and for the control under both penalties. The bar was set from a bias estimate of about $\lambda/\sigma_{\min}^2$, which covers only $L = I$. In our derivation, $(D/h)^{\mathsf T}(D/h)$ has eigenvalues up to about $4N^2$, so the bias in the highest integration mode is about $\lambda \cdot 4N^2/\sigma_{\min}^2 \approx 16 \lambda N^4$, about $10^{-2}$ at N = 512; the bar was mis-set for $L = D/h$. This concerns the check’s design and says nothing either way about the solver.
• C6 failed: five oracle values of $\lambda$ fell on the grid’s lower edge, $10^{-14}$, all for the control: with $L = I$ at N = 64 and 128 and with $L = D/h$ at N = 64, where the printed $\lambda$ ranges start at $10^{-14.0}$. For those three cases the printed Tikhonov errors include seeds whose best weight is at the edge of the grid, essentially no regularization; at N = 64 the control’s smallest and largest Tikhonov errors with $L = I$, 0.02611 and 0.04007, equal its least-squares ones, 0.026115 and 0.040066, to the printed digits.

Conditioning against N. The blur’s singular values at a fixed index barely change with N (at most 0.0085% between N = 256 and 512, and integration’s at most 0.5615%), so the finer matrices approximate the same continuous operators (H1d held). What changes is how many of them there are. The blur’s condition number is $1.57 \times 10^{3}$ at N = 64 and $4.66 \times 10^{13}$ at N = 128; at N = 256 and 512 its smallest singular values, $8.9 \times 10^{-18}$ and $6.6 \times 10^{-18}$, are below double-precision resolution relative to $\sigma_1 \approx 1$ (H1a held). The ranks lstsq keeps under its default cutoff are 64, 127, 130 and 129: from N = 128 on, about 130 directions survive the cutoff however many samples are taken. Even without noise, the default cutoff leaves an error of $1.4 \times 10^{-3}$, $5.7 \times 10^{-2}$ and $6.3 \times 10^{-2}$ at N = 128, 256 and 512: from N = 128 on, lstsq with its default cutoff is a truncated SVD, which by construction drops the part of the scanline that lies in the discarded directions. H1b’s bar asked for a rank at N = 512 at most 1.1 times that at N = 256, which holds, and a rank below N/2 at N = 256 and 512, which does not at N = 256 (130 against 128): H1b failed. Integration’s condition number grows in proportion to N, with $\kappa/N$ between 1.27448 and 1.28284 and doubling ratios 1.99265 to 1.99808, and the control’s stays below 2 (1.998119 at most; H1c held). The octave slopes show the two kinds of decay: integration’s lie between −1.0994 and −0.9779 from $i = 4$ to $N/8$ at every N (H4a held), a $1/i$ law; the blur’s steepen from −0.0082 to −8.1730 over the octaves that stay above $10^{-12}\sigma_1$ (H4b held), the pattern of a decay faster than any power. The starred slopes beyond that point are at the rounding floor.

Unregularized least squares. For integration the error roughly doubles with each doubling of N, means 1.9122, 3.5522, 7.0884 and 14.237, rising at every step in all ten seeds under both cutoffs (H2a held for integration). For the blur the error jumps from about 7 at N = 64 to about $10^{10}$–$10^{11}$ at N = 128, and then does not keep rising: under the default cutoff the means are $5.78 \times 10^{10}$, $2.21 \times 10^{10}$ and $9.00 \times 10^{9}$ at N = 128, 256 and 512, rising at every step in 0 of 10 seeds, and under rcond=1e-15 they are $1.01 \times 10^{11}$, $1.33 \times 10^{12}$ and $6.53 \times 10^{11}$, rising at every step in 1 of 10: H2a failed for the blur under both cutoffs. The control’s error stays near the noise level, with means rising by 4.38%, 7.20% and 9.32% relative to N = 64, inside H2b’s 10% (H2b held), and below the bound $2\sigma_n/\mathrm{rms}(x)$ that its smallest singular value of 0.5 gives (H2c held; our derivation: the expected squared error is $\sigma_n^2 \sum_i 1/\sigma_i^2 \le 4N\sigma_n^2$).

Tikhonov at the oracle weight. Regularized errors stay small at every N: the largest in any row is 0.31990 (integration, N = 64, $L = I$), against least-squares errors up to $10^{12}$ (the control’s rows with $L = I$ at N = 64 and 128 and with $L = D/h$ at N = 64 include oracle weights at the grid edge, as noted under C6). For integration the error falls as N grows, mean ratios 0.8626, 0.9183, 0.9604 ($L = I$) and 0.8405, 0.9883, 0.9104 ($L = D/h$), and in all ten seeds the error at N = 512 is at most that at N = 64 (H3a and H3b held for integration). For the blur, H3a’s bar of at most 5% growth from one N to the next failed: the mean rises by factors of 1.0777 ($L = I$) and 1.0816 ($L = D/h$) from 64 to 128; the later ratios are 0.9546 and 0.9983 for $L = I$ and 0.9655 and 1.0008 for $L = D/h$. H3b, that every seed’s error at N = 512 be at most its error at N = 64, failed for the blur (2 and 1 of 10 seeds) and for the control (1 and 7 of 10), while H3a held for the control (ratios at most 1.0407 for $L = I$ and 1.0196 for $L = D/h$). H3c, a bound of 1.25 on the ratio to N = 64, held for the blur (1.1892, 1.1616) and integration, and failed for the control (1.5020 and 1.4838). The batch plan’s hypothesis, “the oracle-λ Tikhonov error does not” grow with N, thus held only for integration; what the rows show for the blur is a regularized error that in no seed exceeds 1.19 times its N = 64 value (H3c’s 1.1892 and 1.1616) while the unregularized one moves by ten orders of magnitude.

Coefficients. For the blur at N = 256 the data’s coefficients $|u_i^{\mathsf T} y|$ range from 6.7 down to 0.010 for $i \le 32$ and are 0.013, 0.0026 and 0.0065 at $i = 64$, 128 and 256, near the noise level of 0.01, while $\sigma_i$ falls to $1.5 \times 10^{-13}$ and $8.9 \times 10^{-18}$; their ratios, which the unregularized solution multiplies onto the singular vectors, reach $7 \times 10^{14}$. Over the upper half of the indices the coefficients’ root mean square lies between 0.8737 and 1.0589 times $\sigma_n$ for the blur and 0.8805 and 1.0697 for integration at N = 256 and 512 (H5 held): they are noise, which the data checks confirm. Every regularizer above is a way of not dividing those coefficients by small $\sigma_i$.

The example is a model in one dimension, with data generated by the same matrix that is inverted, so it shows the effect of noise alone; a real camera adds model error. Its process is below.

Process. The batch plan’s example entry was turned into this set-up with these departures, each fixed before the first run: (1) the scanline was fixed and “width 0.02” read as the Gaussian’s standard deviation; (2) both operators are discretized by the midpoint rule with weight $h$, integration to each cell’s right edge; (3) “identity plus a small blur” was read as $0.5 I + 0.5 \times$ the blur; (4) noise standard deviation 0.01 and seeds 0–9 were chosen, and the data use the same matrix that is inverted; (5) the error is the relative Euclidean error; (6) Tikhonov is written in continuous units with $L = D/h$, solved by SVD filter factors ($L = I$) and stacked QR ($L = D/h$) on the stated grid with the oracle per seed; (7) decay is measured by octave slopes, for the blur only above $10^{-12}\sigma_1$; (8) the ranks, $\sigma_i$ at fixed indices, integration’s $\kappa/N$, the noise-free blur error and the coefficient RMS were added as measures; (9) the checks D1–D3 and C1–C6 were added; (10) the coefficient table is at N = 256 for seed 0; (11) the plan’s hypotheses were made operational as H1a–H1d, H2a–H2b, H3a–H3b and H4a–H4b, H5 was added, and H2c and H3c were added after a planning script; (12) no rule for choosing $\lambda$ is included, as the entry required. The hypotheses and bars were written before a planning script ran; that script computed expected (noise-averaged) errors on the same matrices and continuous-operator predictions. Its output already showed five bars failing, all kept unchanged: H2a for the blur (expected default-cutoff errors $6.57 \times 10^{10}$, $2.67 \times 10^{10}$ and $9.16 \times 10^{9}$ at N = 128, 256 and 512; with rcond=1e-15 $1.14 \times 10^{12}$ at 256 falling to $8.88 \times 10^{11}$ at 512), H2b (expected ratio 1.162 at N = 512, against which the first run’s seed means held the bar), and H3a and H3b for the blur and the control (expected oracle errors rising from 0.11041 to 0.11814 and from 0.10897 to 0.11721 for the blur between N = 64 and 128, and from 0.03340 to 0.03629 ($L = I$) and 0.03279 to 0.03500 ($L = D/h$) for the control between N = 64 and 128; and, against H3b, expected errors at N = 512 of 0.11440 and 0.11430 for the blur and 0.03866 and 0.03423 for the control, above their N = 64 values). The fifth is H1b: the script printed a rank of 130 at N = 256, which already failed H1b’s bar, although the plan recorded the bar as met. H2c and H3c were written after that script and are illustrations rather than tests. The pilot ran this script unchanged, and the first run, a copy of it, printed the same output; nothing was changed after the first run.

Common Misconceptions

• “Ill-posed means underdetermined.” A square blur with no zeros in its transfer function has a unique solution and is still ill-conditioned, as Image Restoration notes; the example’s blur has full rank at N = 64 and a condition number of $1.6 \times 10^3$.
• “More pixels give a better reconstruction.” Without regularization, the example’s error grows with N for integration and jumps by ten orders of magnitude for the blur between 64 and 128 samples. Finer sampling adds unknowns faster than it adds information when the singular values fall quickly.
• “lstsq gives the unregularized solution.” Its default cutoff, machine precision times $\max(M, N)$ relative to the largest singular value, is a truncated SVD, and the example’s printed ranks show it discarding 126 of the blur’s 256 directions at N = 256 and 383 of 512 at N = 512. State the cutoff used.
• “Regularization is a numerical trick.” In the Bayesian reading it is a prior, and the weight $\lambda$ expresses how much the prior is trusted relative to the noise: for ridge regression Murphy (2022, §11.3) sets $\lambda = \sigma^2/\tau^2$, the noise variance over the prior variance. Image Restoration and maximum likelihood estimation show the correspondence between priors and penalties.
• “A learned prior needs a trained network.” The Deep Image Prior uses only the structure of an untrained network and a limited number of iterations.
• “Every vision problem called ill-posed is ill-posed in the same way.” Some lose information exactly (super-resolution, matting, depth from one image), some only attenuate it (deblurring, integration), and for shape from shading whether the problem is ill-posed at all depends on the model and the boundary information, as the shape-from-shading article reports.

Where It Is Used

• Restoration: image restoration, image deblurring, image denoising, super-resolution, image inpainting, image dehazing, demosaicing, with priors from the Wiener filter, total variation and BM3D.
• Reconstruction: medical image reconstruction, computed tomography and MRI.
• Geometry and motion: shape from shading and variational shape from shading, photometric stereo, stereo matching, monocular depth estimation, and optical flow through brightness constancy and Horn–Schunck smoothness.
• Segmentation and compositing: image matting.
• Numerics: singular value decomposition, least squares, optimization and early stopping supply the tools; scale-space supplies the stable derivative.