Exercises

$\mathrm{rank}(\mathbf{A}) = 2 = \dim(\mathbb{R}^2)$, so $\mathcal{R}(\mathbf{A}) = \mathbb{R}^2$. Existence holds for every $\mathbf{y} \in \mathbb{R}^2$.

By rank-nullity, $\dim(\mathcal{N}(\mathbf{A})) = 3 - 2 = 1$. The null space is $\mathrm{span}\{(6, -3, 1)^T\} \neq \{0\}$. Uniqueness fails: $\mathbf{x}$ and $\mathbf{x} + c(6,-3,1)^T$ give the same $\mathbf{y}$.

On $\mathcal{R}(\mathbf{A}) = \mathbb{R}^2$, the pseudoinverse $\mathbf{A}^\dagger$ is a bounded linear map (any linear map on a finite-dimensional space is bounded). Stability holds.

The moral: in finite dimensions, the Hadamard stability condition is automatically satisfied. Ill-conditioning is a continuous-limit phenomenon captured by the condition number, not by Hadamard's framework itself.

$\sum_k \frac{|1/k^3|^2}{(1/k)^2} = \sum_k \frac{k^2}{k^6} = \sum_k k^{-4} < \infty$. Holds.

$\sum_k \frac{|1/k|^2}{(1/k)^2} = \sum_k 1 = \infty$. Fails.

$\sum_k \frac{|1/k^2|^2}{(1/k)^2} = \sum_k k^{-2} < \infty$. Holds. (The alternating sign does not affect convergence.)

$\mathbf{A}^{-1}\mathbf{y}^\delta = \begin{pmatrix} 3.1/3 \\ 0.5/0.01 \end{pmatrix} = \begin{pmatrix} 1.033 \\ 50 \end{pmatrix}.$$The second component is amplified by$1/0.01 = 100$.

$\mathbf{A}^T\mathbf{A} + 0.01\mathbf{I} = \begin{pmatrix} 9.01 & 0 \\ 0 & 0.0002 \end{pmatrix}, \qquad \mathbf{A}^T\mathbf{y}^\delta = \begin{pmatrix} 9.3 \\ 0.005 \end{pmatrix}.KATEXPLACEHOLDER0END\mathbf{x}_\alpha = \begin{pmatrix} 9.3/9.01 \\ 0.005/0.0002 \end{pmatrix} = \begin{pmatrix} 1.032 \\ 25 \end{pmatrix}.$$The well-determined component ($\sigma_1 = 3$) is barely affected. The ill-determined component ($\sigma_2 = 0.01$) is reduced from 50 to 25 — Tikhonov damps the amplification by the Tikhonov filter$\sigma_2^2/(\sigma_2^2 + \alpha) = 0.0001/0.0101 \approx 0.01$.

$\langle x^\dagger, v_k\rangle = (1/k^2)/(1/k) = 1/k$.

$e_K^2 = \sum_{k=K+1}^{100} \frac{1}{k^2} + \sum_{k=1}^{K} \frac{(0.01)^2}{(1/k)^2} = \sum_{k=K+1}^{100} k^{-2} + 10^{-4} \sum_{k=1}^{K} k^2.$$The bias$\approx 1/K$and variance$\approx 10^{-4} K^3/3$. Setting bias = variance:$1/K^2 \approx 10^{-4} K^2$, giving$K^4 \approx 10^4$, so$K^* \approx 10$.

The optimal Tikhonov cutoff occurs at $\sigma_K \approx \sqrt{\alpha}$, corresponding to $K^* \approx 10$ and $\sigma_{10} = 0.1$, so $\alpha^* \approx 0.01$.

Numerically, the TSVD error is slightly smaller than Tikhonov's because TSVD has infinite qualification while Tikhonov's qualification saturates at $\mu_0 = 2$. For this problem with $\mu \approx 1$ (solution coefficients $\sim 1/k$), both methods are near-optimal and the difference is modest.

$\mathbf{A}\mathbf{x}_\alpha - \mathbf{y}^\delta = -\sum_k \frac{\alpha}{\sigma_k^2 + \alpha} \langle \mathbf{y}^\delta, \mathbf{u}_k\rangle\, \mathbf{u}_k,$$giving$\varphi(\alpha) = \sum_k (\alpha/(\sigma_k^2+\alpha))^2 |\langle \mathbf{y}^\delta, \mathbf{u}_k\rangle|^2$.

Each term $f_k(\alpha) = (\alpha/(\sigma_k^2+\alpha))^2 |\langle \mathbf{y}^\delta, \mathbf{u}_k\rangle|^2$ has derivative

$$f_k'(\alpha) = \frac{2\alpha\sigma_k^2}{(\sigma_k^2+\alpha)^3} |\langle \mathbf{y}^\delta, \mathbf{u}_k\rangle|^2 \geq 0,$$

with equality only if $\alpha = 0$ or $\sigma_k = 0$. Since $\varphi = \sum_k f_k$, we have $\varphi'(\alpha) > 0$.

$\varphi(0) = 0$ (exact data case) and $\varphi(\infty) = \|\mathbf{y}^\delta\|^2$. By the intermediate value theorem, the equation $\varphi(\alpha) = \tau^2\delta^2$ has exactly one solution whenever $\tau\delta < \|\mathbf{y}^\delta\|$. $\blacksquare$

In the SVD basis, the $k$-th component of $x_n$ satisfies the recurrence (starting from $c_k^{(0)} = 0$):

$$c_k^{(n+1)} = c_k^{(n)} + \omega\sigma_k(d_k - \sigma_k c_k^{(n)}),$$

where $d_k = \langle y^\delta, u_k\rangle$. This gives $c_k^{(n)} = (1 - (1-\omega\sigma_k^2)^n)d_k/\sigma_k$, confirming the filter $F_n(\sigma^2) = (1-(1-\omega\sigma^2)^n)/\sigma^2$.

$|\sigma F_n(\sigma^2)| = |1 - (1-\omega\sigma^2)^n|/\sigma$. Using $1 - (1-t)^n \leq nt$ for $t \in [0,1]$ (with $t = \omega\sigma^2$): $|\sigma F_n| \leq n\omega\sigma$. Maximising over $\sigma \in (0, 1/\sqrt{\omega}]$ gives maximum at $\sigma = 1/\sqrt{\omega}$, yielding $|\sigma F_n| \leq n\sqrt{\omega}$. A tighter analysis using calculus gives $\sup |\sigma F_n| \leq \sqrt{n\omega}$. $\blacksquare$

We need $\langle w, v_k\rangle = \langle x^\dagger, v_k\rangle / \sigma_k = \langle x^\dagger, v_k\rangle \cdot k^2$ to be in $\ell^2$, i.e., $\sum_k k^4 |\langle x^\dagger, v_k\rangle|^2 < \infty$.

(a) $\sum_k k^4 (k^{-3})^2 = \sum_k k^{-2} < \infty$. Satisfies $\mu = 1$. $\checkmark$

(b) $\sum_k k^4 (k^{-2})^2 = \sum_k 1 = \infty$. Does not satisfy $\mu = 1$.

(c) $\sum_k k^4 (k^{-5})^2 = \sum_k k^{-6} < \infty$. Satisfies $\mu = 1$. (In fact, case (c) satisfies a higher-order source condition $\mu = 1.5$.) $\checkmark$

The KKT condition at component $i$ is $[\mathcal{A}^*r]_i + \lambda s_i = 0$ where $s_i \in \partial|\hat{x}_i|$. If $\hat{x}_i = 0$, then $s_i \in [-1, 1]$ and $|[\mathcal{A}^*r]_i| = \lambda|s_i| \leq \lambda$.

Conversely, if $|[\mathcal{A}^*r]_i| < \lambda$, we cannot have $\hat{x}_i \neq 0$ (which would require $|s_i| = 1$ and $|[\mathcal{A}^*r]_i| = \lambda$). Hence $\hat{x}_i = 0$.

The noise contribution to $\mathcal{A}^*\eta$ has each component $[\mathcal{A}^*\eta]_i \sim \mathcal{N}(0, \|\mathcal{A}^* e_i\|^2 \sigma_n^2) \approx \mathcal{N}(0, \sigma_n^2)$ (for near-isometric $\mathcal{A}^*$).

By the maximum of $n$ i.i.d. Gaussians: $\mathbb{P}(\max_i |[\mathcal{A}^*\eta]_i| > \sigma_n\sqrt{2\log n + t}) \leq e^{-t/2}$.

Choosing $\lambda = \sigma_n\sqrt{2\log n}$ ensures all noise-only components are suppressed with high probability.

$\frac{dR^2}{d\alpha} = \frac{d}{d\alpha}\sum_k \left(\frac{\alpha}{\sigma_k^2+\alpha}\right)^2 |\langle y^\delta, u_k\rangle|^2 = 2\sum_k \frac{\alpha\sigma_k^2}{(\sigma_k^2+\alpha)^3} |\langle y^\delta, u_k\rangle|^2 > 0.$$

$\frac{dS^2}{d\alpha} = -2\sum_k \frac{\sigma_k^2}{(\sigma_k^2+\alpha)^3} |\langle y^\delta, u_k\rangle|^2 < 0.$$

$\rho'(\alpha) = \frac{1}{2R^2}\frac{dR^2}{d\alpha}, \qquad \eta'(\alpha) = \frac{1}{2S^2}\frac{dS^2}{d\alpha}.KATEXPLACEHOLDER0END\kappa(\alpha) = \frac{|\rho'\eta'' - \rho''\eta'|} {(\rho'^2 + \eta'^2)^{3/2}},$$where second derivatives follow by the quotient rule applied to the expressions above. This curvature is maximised at the L-curve corner.$\blacksquare$

The bias for the $k$-th mode is $\alpha/(\sigma_k^2+\alpha) \cdot \sigma_k^\mu|\langle w, v_k\rangle|$. We need $g(\sigma) = \alpha\sigma^\mu/(\sigma^2+\alpha)$.

For $\mu \leq 2$: the maximum occurs at $\sigma^* = \sqrt{\alpha\mu/(2-\mu)}$ (or $\sigma \to 0$ for $\mu = 2$), giving $\sup g = C_\mu \alpha^{\mu/2}$.

For $\mu > 2$: $g$ is monotone increasing and $\sup_{\sigma \leq \|\mathcal{A}\|} g = \alpha\|\mathcal{A}\|^\mu/(\|\mathcal{A}\|^2+\alpha) \leq \alpha\|\mathcal{A}\|^{\mu-2}$, so bias $= O(\alpha)$.

With bias $\leq C\alpha^{\mu/2}E$ and variance $\leq \delta^2/(4\alpha)$ (from $\sigma/(\sigma^2+\alpha) \leq 1/(2\sqrt{\alpha})$), optimising $\alpha^{\mu/2}E + \delta/\sqrt{\alpha}$ gives $\alpha^* \sim (\delta/E)^{2/(2\mu+1)}$ and rate $O(\delta^{2\mu/(2\mu+1)})$.

With bias $= O(\alpha E)$ and variance $= O(\delta/\sqrt{\alpha})$, optimising $\alpha E + \delta/\sqrt{\alpha}$ gives $\alpha^* \sim (\delta/E)^{2/3}$ and rate $O(\alpha^* E) = O((\delta E^2)^{1/3})$. More carefully (working with squared errors): bias $= O(\alpha^2 E^2)$, variance $= O(\delta^2/\alpha)$. Optimising: $2\alpha E^2 = \delta^2/\alpha^2$ gives $\alpha^3 = \delta^2/(2E^2)$, $\alpha^* \sim (\delta/E)^{2/3}$, and error $\sim (\delta/E)^{4/3}\cdot E^2 \cdot E^{-4/3}$, which in the $\|\cdot\|$ norm gives $O(\delta^{4/5})$ (the exponent $2\cdot2/(2\cdot2+1) = 4/5$). $\blacksquare$

Write the estimator as $\hat{y} = H_\alpha y^\delta$ where $H_\alpha = \mathcal{A}(\mathcal{A}^*\mathcal{A}+\alpha I)^{-1}\mathcal{A}^*$. The residual is $(I - H_\alpha)y^\delta$.

By Stein's lemma, $\mathbb{E}[\|(I - H_\alpha)y^\delta\|^2] = \|(\mathcal{I} - H_\alpha)\mathcal{A}x^\dagger\|^2 + \sigma_n^2 \mathrm{tr}((I - H_\alpha)^2) + 2\sigma_n^2 \mathrm{tr}(H_\alpha - H_\alpha^2) - m\sigma_n^2$.

Rearranging recovers the SURE formula. The expectation of SURE equals the prediction MSE.

For $m = n = 1$: $\widehat{\mathrm{MSE}}(\alpha) = -\sigma_n^2 + (\alpha d/(\sigma^2+\alpha))^2 + 2\sigma_n^2\sigma^2/(\sigma^2+\alpha)$.

Setting the derivative to zero and solving gives

$$\alpha^* = \frac{\sigma^2(\sigma_n^2/d^2 - 1)}{1 - \sigma_n^2/d^2} = \sigma_n^2 \sigma^2 / (d^2 - \sigma_n^2)$$

when $|d| > \sigma_n$ (the data exceeds the noise). For $|d| \leq \sigma_n$ (pure noise), $\alpha^* \to \infty$ and the solution is set to zero.

For $w \in L^2(\text{receiver domain})$:

$$[(\mathcal{F}'(x_n))^* w](\mathbf{r}) = k_0^2 \overline{E_n(\mathbf{r})} \int_{\text{rx}} \overline{G_0(\mathbf{r}_{\text{rx}}, \mathbf{r})} w(\mathbf{r}_{\text{rx}})\,d\mathbf{r}_{\text{rx}}.$$

This is a backpropagation: field $w$ recorded at receivers is backpropagated to domain $\Omega$ via the conjugate Green's function, then weighted by $\overline{E_n}$.

The IRGNM normal equation is $(\mathcal{F}'(x_n)^*\mathcal{F}'(x_n) + \alpha_n I)\delta\chi = \mathcal{F}'(x_n)^*(y^\delta - \mathcal{F}(x_n))$.

The $(i,j)$ entry of the system matrix involves a volume integral $\int_\Omega G_0^*(\mathbf{r}_{\text{rx}},\mathbf{r}_i)E_n(\mathbf{r}_i) G_0(\mathbf{r}_{\text{rx}},\mathbf{r}_j)E_n(\mathbf{r}_j)d\mathbf{r}_{\text{rx}}$, a sum over receivers of outer products of the forward Green's functions.

The forward computation $\mathcal{F}'(x_n)\delta\chi$ evaluates the field at $N_{\text{rx}}$ receiver positions, each requiring $O(N_\Omega)$ work (one Green's function integral per receiver per domain point): total $O(N_{\text{rx}} N_\Omega)$.

The adjoint $(\mathcal{F}'(x_n))^* w$ maps a vector of length $N_{\text{rx}}$ to a function on $\Omega$ of size $N_\Omega$: total $O(N_{\text{rx}} N_\Omega)$ as well — the same cost.

The key is that one forward/adjoint pair costs $2 \times O(N_{\text{rx}} N_\Omega)$ rather than $O(N_{\text{rx}}^2 N_\Omega)$ for the full Jacobian. This is the basis of the efficient IRGNM implementation.

The log-posterior is

$$\log p(x|y^\delta) = -\frac{1}{2\sigma_n^2}\|y^\delta - \mathcal{A}x\|^2 - \frac{1}{2}\langle x-x_0, \mathcal{C}_{\text{prior}}^{-1}(x-x_0)\rangle + \mathrm{const}.$$

In the SVD basis, the $k$-th component satisfies:

$$\langle x_{\text{MAP}}, v_k\rangle = \frac{\sigma_k^2/\sigma_n^2}{\sigma_k^2/\sigma_n^2 + 1/\lambda_k^{\text{pr}}} \cdot \frac{\langle y^\delta, u_k\rangle}{\sigma_k} + \frac{1/\lambda_k^{\text{pr}}}{\sigma_k^2/\sigma_n^2 + 1/\lambda_k^{\text{pr}}} \cdot \langle x_0, v_k\rangle.$$

This is a weighted average of the data-driven and prior-mean estimates.

The $k$-th eigenvalue of $\mathcal{C}_{\text{post}}$ is

$$\lambda_k^{\text{post}} = \frac{1}{\sigma_k^2/\sigma_n^2 + 1/\lambda_k^{\text{pr}}}.$$

When $\sigma_k^2 \gg \sigma_n^2/\lambda_k^{\text{pr}}$ (strong signal): $\lambda_k^{\text{post}} \approx \sigma_n^2/\sigma_k^2$ — data-determined, uncertainty set by noise.

When $\sigma_k^2 \ll \sigma_n^2/\lambda_k^{\text{pr}}$ (weak signal): $\lambda_k^{\text{post}} \approx \lambda_k^{\text{pr}}$ — prior-determined, data provides no information.

As $\lambda_k^{\text{pr}} \to \infty$, $1/\lambda_k^{\text{pr}} \to 0$ and the regularisation parameter $\alpha_k = \sigma_n^2/\lambda_k^{\text{pr}} \to 0$. The MAP estimate becomes $\langle x_{\text{MAP}}, v_k\rangle \to \langle y^\delta, u_k\rangle/\sigma_k = \langle \mathcal{A}^\dagger y^\delta, v_k\rangle$, recovering the pseudoinverse. $\blacksquare$

For a single step of height $h$ at position $k$, any piecewise-constant solution $\hat{x}$ taking value $a$ on $\{i < k\}$ and $b$ on $\{i \geq k\}$ has $\mathrm{TV}(\hat{x}) = |b - a|$.

The ROF functional becomes $J(a,b) = \frac{k}{2}a^2 + \frac{n-k}{2}(b-h)^2 + \lambda|b-a|$.

Since we can independently choose $a$ and $b$ (the model is decoupled at the edge), the optimal $a = 0$ (minimum $\ell_2$ from 0) and $b$ satisfies $\min_b \frac{n-k}{2}(b-h)^2 + \lambda|b|$.

$\min_b \frac{1}{2}(b-h)^2 + \frac{\lambda}{n-k}|b|$ is a soft thresholding problem with threshold $\lambda/(n-k)$. For a single step ($n-k = 1$): solution is $b^* = \mathrm{sign}(h)(|h| - \lambda)^+$.

Therefore $\hat{x}_i = 0$ for $i < k$ and $\hat{x}_i = (h-\lambda)^+ \cdot \mathrm{sign}(h)$.

The minimum surviving edge height is $h_{\min} = \lambda$. Edges with $|h| > \lambda$ are preserved (with height reduced by $\lambda$); weaker features are suppressed.

For Tikhonov regularisation, the step would be blurred over $O(1/\sqrt{\alpha})$ pixels — the edge is not preserved, but spread out. For radar scenes with discrete targets or building edges, this blurring destroys the spatial precision needed for target localisation. TV preserves the edge location exactly. $\blacksquare$

Using the same grid and solver means $y = \mathcal{A}_h\chi^\dagger$ exactly lies in the range of the discrete forward operator $\mathcal{A}_h$. The reconstruction algorithm can therefore find $\hat{\chi}$ such that $\mathcal{A}_h\hat{\chi} = y$ exactly (to machine precision), regardless of whether this corresponds to a physically correct solution.

The inverse crime inflates performance: discretization errors (which would appear in reality when the forward model is imperfect) are absent, giving a reconstruction quality that far exceeds what any real system could achieve. Published error metrics from such simulations are misleading.

1. Different grids: Generate data on a fine grid ($200 \times 200$) and reconstruct on a coarser grid ($100 \times 100$). The grid mismatch introduces a model error comparable to measurement noise.

2. Different solvers: Use FEM for data generation and a simpler Born series for reconstruction. The model mismatch reflects the real-world imperfection of the forward model used in practice.

3. Measured calibration data: Use a calibrated measurement setup rather than a pure simulation.

With a fine forward model (high-accuracy data generation) and a coarser reconstruction grid, the mismatch error is $O(h_{\text{coarse}}^p)$ for a $p$-th order solver. If this dominates the measurement noise $\delta$, the effective noise level seen by the reconstruction algorithm is $\delta_{\text{eff}} \approx h_{\text{coarse}}^p$, which requires a larger $\alpha$ than the measurement noise alone would suggest.

In practice, the regularization parameter must be tuned to balance both measurement noise and model error — not just noise.