Part II · Module 13

Inverse problems and regularization

Sound-speed maps, deconvolved RF lines and sparse-recovery beamformers are all solutions of \(\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{n}\) with an operator that throws information away. This module shows where the information goes, how Tikhonov, iteration count, L1 and total variation put prior knowledge back, how to set their strength without the truth, and what the answer can and cannot resolve. It ends with a limited-angle echo-mode sound-speed tomography you can break.

You will be able to

Diagnose ill-posedness from the singular spectrum and name the null space of an ultrasound operator.
Write Tikhonov in filter-factor form, convert between the λ and λ² conventions, and predict its noise gain and bias.
Choose λ by the discrepancy principle, GCV or the L-curve, and say when each fails.
Use iteration count, L1 (ISTA/FISTA) and total variation as regularizers and state the bias each introduces.
Read a model-resolution matrix and explain the axial smearing of limited-angle echo-mode sound-speed images.

Assumes

Module 05 · Matrices as operators: least squares, κ, null space
Module 06 · Eigendecomposition & SVD: SVD, truncated SVD and the Picard plot (Demo 06.5)
Module 11 · Correlation & time-delay estimation: delay estimation between angles
Module 12 · Estimation theory & the Cramér–Rao bound: likelihood and MAP, briefly

Exact formulaclosed form plotted as isSimulationinversion or Monte Carlo run here (JS in the page, or the named Octave script)Approx. modelsimplified physics (Gaussian pulse spectrum, vertical-path average, straight rays linear in slowness); assumptions stated under the badgeSchematicdrawing for intuition, not computedNo curve on this page is hand-drawn: every plotted line comes from the formula or the solver named next to it.

Convention used in this module

The Tikhonov penalty is written \(\lambda^2\lVert\mathbf{L}\mathbf{x}\rVert^2\), so the filter factor is \(f_i = \sigma_i^2/(\sigma_i^2 + \lambda^2)\) and \(\lambda\) has the units of a singular value (Hansen 2010; Module 06 uses the same form with α). Module 05 wrote the penalty as \(\lambda\lVert\mathbf{x}\rVert^2\), whose λ is our \(\lambda^2\). Ridge regression, many ML texts and most "diagonal loading" values also quote the squared quantity. Before you copy a regularization value from a paper, find out which one it is: the two differ by a square, which is a factor 100 for λ = 0.01.

13.1Forward models and why inversion fails

Every model-based method in the target areas starts from a linear (or linearized) forward model and asks for \(\mathbf{x}\) back:

\[ \begin{gathered} \mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{n}, \qquad \mathbf{A} \in \mathbb{R}^{m\times n} \\ \mathbf{n} \sim \mathcal{N}(0,\ \sigma_n^2\mathbf{I}) \end{gathered} \]

(13.1)

\(\mathbf{y}\): data: m samples (RF line, delays between angles, average sound speeds)
\(\mathbf{x}\): unknowns: n values (reflectivity, slowness map, bubble amplitudes)
\(\mathbf{A}\): forward operator (the physics), m × n
\(\sigma_n\): noise standard deviation, units of \(\mathbf{y}\)
\(\mathbf{n}\): noise vector, m values, same units as \(\mathbf{y}\)
\(m,\ n\): number of data samples and number of unknowns
\(\mathbf{I}\): identity matrix
\(\mathcal{N}(0,\ \sigma_n^2\mathbf{I})\): zero-mean white Gaussian distribution, variance \(\sigma_n^2\)
\(\mathbb{R}^{m\times n}\): real matrices with m rows and n columns

\(\mathbf{A}\) is the physics (convolution with the pulse, ray integrals of slowness, a PSF dictionary); \(\mathbf{n}\) is electronic noise plus everything the model leaves out (speckle decorrelation, refraction, model error).

Problem	A	\(\mathbf{x}\)	What A loses
RF deconvolution	convolution with the pulse (Toeplitz; circulant if padded)	reflectivity along depth	DC and everything outside the transducer band: \(\sigma = \lvert H(f)\rvert \approx 0\) there
Average-to-local sound speed	running mean of slowness with depth	local slowness \(s(z)\)	fine detail at depth: \(\sigma_i\) decay, differentiation amplifies noise \(\propto z\)
Echo-mode SoS tomography	differences of straight-ray integrals between angle pairs	slowness map \(\Delta s(x, z)\)	axial (\(k_z\)) structure and laterally invariant layers
Sparse beamforming, dense ULM	PSF dictionary on a fine grid	scatterer or bubble amplitudes	uniqueness: more unknowns than data (\(m < n\)), neighbouring columns nearly parallel

Hadamard's three conditions. A problem is well posed if a solution exists, is unique, and depends continuously on the data. Discretized imaging problems usually satisfy the first two formally (least squares exists; add a minimum-norm condition for uniqueness) and fail the third in practice: the inverse exists, but its amplification is so large that the answer is noise. The SVD shows it in one line (Module 06):

\[ \begin{aligned} \mathbf{x}_{\mathrm{naive}} = \mathbf{A}^{+}\mathbf{y} &= \sum_i \frac{\mathbf{u}_i^{T}\mathbf{y}}{\sigma_i}\,\mathbf{v}_i \\ &= \mathbf{x} + \sum_i \frac{\mathbf{u}_i^{T}\mathbf{n}}{\sigma_i}\,\mathbf{v}_i \\ \operatorname{std}(\text{noise term } i) &= \sigma_n/\sigma_i \end{aligned} \]

(13.2)

\(\sigma_i,\ \mathbf{u}_i,\ \mathbf{v}_i\): singular values and left/right singular vectors of \(\mathbf{A}\) (Module 06)
\(\mathbf{A}^{+}\): pseudo-inverse
\(\mathbf{x}_{\mathrm{naive}}\): unregularized estimate of \(\mathbf{x}\)
\(\mathbf{y},\ \mathbf{x},\ \mathbf{n}\): as in (13.1)
\(\mathbf{u}_i^{T}\): transpose of \(\mathbf{u}_i\); \(\mathbf{u}_i^{T}\mathbf{y}\) is the projection of the data on \(\mathbf{u}_i\)
\(\operatorname{std}(\cdot)\): standard deviation (same units as \(\mathbf{x}\))
\(\sigma_n\): as in (13.1)

White noise puts the same \(\sigma_n\) into every coefficient \(\mathbf{u}_i^{T}\mathbf{n}\); dividing by small \(\sigma_i\) blows it up. The continuous analogue: smoothing (compact) operators have \(\sigma_i \to 0\), so their inverse is unbounded. The discrete Picard condition of Module 06 says where the signal coefficients \(\lvert\mathbf{u}_i^{T}\mathbf{y}\rvert\) sink into the noise floor; we do not redo that plot here.

Ultrasound link

The three operators above fail differently. The pulse has an exact null space (frequencies it never transmits), the running mean has a smooth decay, and the limited-angle tomography operator loses a whole cone of spatial frequencies. Averaging more frames lowers \(\sigma_n\) but cannot recover a null space: only a new measurement (more bandwidth, more angles) or a prior can.

Pitfall

A small residual is not evidence of a good solution. The naive inverse fits the data, noise included, exactly. Judge an inversion by its error against a known truth in simulation, by its resolution matrix, and by its stability to a new noise draw.

13.2Tikhonov regularization: spectral filtering, bias and resolution

Add a penalty on the size of the solution and minimize both:

\[ \begin{aligned} \mathbf{x}_\lambda &= \arg\min_{\mathbf{x}}\ \lVert \mathbf{A}\mathbf{x} - \mathbf{y} \rVert^2 + \lambda^2\lVert\mathbf{x}\rVert^2 \\ &= (\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{I})^{-1}\mathbf{A}^{T}\mathbf{y} \\ &= \sum_i f_i\,\frac{\mathbf{u}_i^{T}\mathbf{y}}{\sigma_i}\,\mathbf{v}_i \\ f_i &= \frac{\sigma_i^2}{\sigma_i^2 + \lambda^2} \end{aligned} \]

(13.3)

\(\lambda\): regularization parameter, units of a singular value (Hansen convention)
\(f_i\): filter factor of component i, 0 to 1
\(\sigma_i,\ \mathbf{u}_i,\ \mathbf{v}_i\): SVD of \(\mathbf{A}\)
\(\mathbf{x}_\lambda\): regularized (Tikhonov) estimate
\(\arg\min\): the \(\mathbf{x}\) that minimizes the expression
\(\lVert\cdot\rVert\): Euclidean (2-) norm
\(\mathbf{A},\ \mathbf{y},\ \mathbf{x},\ \mathbf{I}\): as in (13.1)
\(^{T},\ ^{-1}\): matrix transpose, matrix inverse

Equivalent stacked least squares: \(\min \lVert[\mathbf{A};\ \lambda\mathbf{I}]\,\mathbf{x} - [\mathbf{y};\ 0]\rVert\), which you solve by QR without forming \(\mathbf{A}^{T}\mathbf{A}\) (Module 05). \(f_i \approx 1\) for \(\sigma_i \gg \lambda\), \(f_i \approx \sigma_i^2/\lambda^2\) for \(\sigma_i \ll \lambda\), and exactly ½ at \(\sigma_i = \lambda\): λ is a soft singular-value threshold. Truncated SVD is the hard version.

Split the error into what the filter removes from the truth and what it lets through from the noise. With \(\mathbf{A}^{\#} = (\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{I})^{-1}\mathbf{A}^{T}\) and the model-resolution matrix \(\mathbf{R} = \mathbf{A}^{\#}\mathbf{A} = \mathbf{V}\operatorname{diag}(f)\,\mathbf{V}^{T}\):

\[ \begin{gathered} \mathbf{x}_\lambda - \mathbf{x} = (\mathbf{R} - \mathbf{I})\,\mathbf{x} + \mathbf{A}^{\#}\mathbf{n} \\ \text{noise gain per component:} \\ \frac{f_i}{\sigma_i} = \frac{\sigma_i}{\sigma_i^2 + \lambda^2} \le \frac{1}{2\lambda} \end{gathered} \]

(13.4)

\(\mathbf{R}\): model-resolution matrix \(\mathbf{A}^{\#}\mathbf{A}\); its rows are the reconstruction's point-spread functions
\(\mathbf{A}^{\#}\): the regularized inverse, \((\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{I})^{-1}\mathbf{A}^{T}\)
\(\mathbf{x}_\lambda,\ \mathbf{x},\ \mathbf{n}\): regularized estimate, true unknowns and noise, as in (13.1) and (13.3)
\(\mathbf{I}\): identity matrix
\(f_i,\ \sigma_i,\ \lambda\): as in (13.3)
\(\lvert\cdot\rvert,\ \le\): noise gain per component is bounded by \(1/(2\lambda)\), peak at \(\sigma_i = \lambda\)

The first term is the regularization error (a bias: \(\mathrm{E}[\mathbf{x}_\lambda] = \mathbf{R}\mathbf{x}\), not \(\mathbf{x}\)); the second is the perturbation error. Raising λ trades the second for the first. The cap \(1/(2\lambda)\), reached at \(\sigma_i = \lambda\), is the most any component of the noise can be amplified.

Deconvolution is the circulant case. For a pulse \(h\) and a periodic (or padded) line, the SVD is the DFT and \(\sigma = \lvert H(f)\rvert\). Tikhonov becomes a frequency-domain filter, and \(\mathbf{R}\) becomes a convolution: the deconvolved line is the true reflectivity convolved with a new pulse \(r(t)\) whose spectrum is the filter factor.

\[ \begin{gathered} X_\lambda(f) = \frac{H^*(f)\,Y(f)}{\lvert H(f)\rvert^2 + \lambda^2} \\ \mathbf{x}_\lambda = r * \mathbf{x} + \text{noise} \\ R(f) = \frac{\lvert H(f)\rvert^2}{\lvert H(f)\rvert^2 + \lambda^2} \end{gathered} \]

(13.5)

\(H(f)\): pulse spectrum (two-way transfer function), normalized to peak 1
\(Y(f),\ X_\lambda(f)\): spectra of the RF line and of the estimate
\(r(t),\ R(f)\): resolution kernel and its spectrum
\(H^*(f)\): complex conjugate of \(H(f)\)
\(\lvert H(f)\rvert\): magnitude of the pulse spectrum
\(\lambda\): regularization parameter, in units of \(\lvert H\rvert\) (peak 1)
\(\mathbf{x}_\lambda,\ \mathbf{x}\): estimate and true reflectivity series
\(*\): convolution in time
\(f\): frequency, MHz

This is the Wiener deconvolution filter for a white reflectivity of variance \(\sigma_x^2\) in white noise, with \(\lambda^2 = \sigma_n^2/\sigma_x^2\) (13.6 gives the Bayesian reason). \(\operatorname{trace}\mathbf{R} = \sum f_i\) counts the effective number of parameters the data determine.

Demo 13.1

Tikhonov deconvolution of an RF line: solution, residual, filter factors, resolution kernel

SimulationInversion run in the page (JS, FFT): Gaussian-spectrum 5 MHz pulse, random reflectivity, white noise; filter factors from (13.3).

true reflectivityx_λ / filter factors / kernelRF data yresidual y − Ax_λ

λ (× σ_max)

SNR of the RF line

Try: 5 MHz pulse, 60 % fractional bandwidth (−6 dB), \(f_s\) = 40 MHz, 512 samples (9.9 mm). Push λ to \(10^{-4}\): noise in the band edges where \(\lvert H\rvert \approx 10^{-4}\) is amplified up to \(1/(2\lambda)\) = 5000 times and the estimate fills with high-frequency garbage while the residual vanishes. Push λ to 1: only the band centre survives, the kernel is as long as the pulse, and the residual reaches the signal level. The error is measured against the band-limited truth (the reflectivity filtered to where \(\lvert H\rvert \ge 0.01\), about 1.1–8.9 MHz), because nothing outside that band is recoverable without a prior.

Intuition

Read the filter-factor panel as "which frequencies of the reflectivity do I believe?" Where the pulse is stronger than λ the estimate takes the data at face value; where it is weaker the estimate returns roughly zero. The passband edges where \(\lvert H\rvert = \lambda\) set the new axial resolution, so deconvolution buys resolution only up to the frequency where the pulse sinks to the noise.

Pitfall

"Super-resolved" deconvolution with a quadratic penalty is band widening up to the \(\lvert H\rvert = \lambda\) crossing, never beyond. Recovering detail outside the band needs a prior that extrapolates (sparsity, positivity, piecewise constancy; 13.6–13.7), and then the result is only as right as the prior.

13.3General-form Tikhonov: priors on shape

Penalize a derivative instead of the size, and the regularizer prefers a shape instead of zero:

\[ \begin{aligned} \mathbf{x}_\lambda &= \arg\min_{\mathbf{x}}\ \lVert \mathbf{A}\mathbf{x} - \mathbf{y} \rVert^2 + \lambda^2\lVert\mathbf{L}\mathbf{x}\rVert^2 \\ &= (\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{L}^{T}\mathbf{L})^{-1}\mathbf{A}^{T}\mathbf{y} \\ &\text{unique iff } \mathcal{N}(\mathbf{A}) \cap \mathcal{N}(\mathbf{L}) = \{0\} \end{aligned} \]

(13.6)

\(\mathbf{L}\): penalty operator: \(\mathbf{I}\), first difference \(\mathbf{D}_1\), second difference \(\mathbf{D}_2\), or a 2D gradient
\(\mathcal{N}(\cdot)\): null space
\(\gamma_i\): generalized singular values of \((\mathbf{A}, \mathbf{L})\)
\(\mathbf{x}_\lambda,\ \mathbf{A},\ \mathbf{y},\ \lambda\): as in (13.3)
\(\{0\}\): the zero vector alone; the null spaces of \(\mathbf{A}\) and \(\mathbf{L}\) share no other vector
\(\cap\): intersection of the two null spaces

\(\mathbf{L} = \mathbf{I}\) pulls toward \(\mathbf{x} = 0\) (here: toward the reference \(c_0\)); \(\mathbf{L} = \mathbf{D}_1\) (first difference) leaves constants free and pulls toward flat profiles; \(\mathbf{L} = \mathbf{D}_2\) leaves straight lines free. The generalized SVD of \((\mathbf{A}, \mathbf{L})\) gives filter factors \(\gamma_i^2/(\gamma_i^2 + \lambda^2)\), so everything in 13.2 carries over.

A one-dimensional sound-speed example. Depth-resolved autofocus gives an average sound speed \(\bar c(z)\). In the vertical-path approximation the average is a travel-time average, which is linear in slowness \(s = 1/c\):

\[ \begin{aligned} \frac{1}{\bar c(z)} &= \frac{1}{z}\int_0^z s(z')\,dz' \\ \Rightarrow\ s(z_i) &= \frac{z_i/\bar c(z_i) - z_{i-1}/\bar c(z_{i-1})}{\Delta z} \\ \operatorname{std}(c) &\approx \left(\frac{c}{\bar c}\right)^{\!2} \frac{\sqrt{z_i^2 + z_{i-1}^2}}{\Delta z}\,\sigma_{\bar c} \\ &\approx \sqrt{2}\,\frac{z}{\Delta z}\,\sigma_{\bar c} \end{aligned} \]

(13.7)

\(\bar c(z)\): average (travel-time) sound speed from the surface to depth z, m/s
\(s(z)\): local slowness \(1/c\), s/m
\(\Delta z\): depth sample spacing, mm
\(\sigma_{\bar c}\): std of each average-sound-speed estimate, m/s
\(z,\ z_i\): depth and depth of sample i, mm
\(z\prime\): integration variable (depth), mm
\(c\): local sound speed, m/s
\(\operatorname{std}\): standard deviation, m/s

The naive inverse differentiates, and the noise gain grows linearly with depth. For Δz = 0.4 mm and \(\sigma_{\bar c}\) = 1 m/s: \(\sqrt{2}\,z/\Delta z\) gives 35, 71 and 106 m/s at 10, 20 and 30 mm, and 40, 75 and 109 m/s with the factor \((c/\bar c)^2\) (1.03–1.15 for the layered profile of m13_tv_vs_tikhonov.m, which measures 39, 77, 111). The discretized operator (100 samples) has κ = 310. Real average-sound-speed estimators average over the aperture's angular range, so model-based methods use a more complete forward model (for example Jakovljevic et al., JASA, 2018); the vertical-path version keeps the algebra visible.

Demo 13.2

Average-to-local sound speed: naive inverse vs Tikhonov with L = I, D₁, D₂

SimulationApprox. modelVertical-path travel-time average 1/c̄ = (1/z)∫s dz, independent Gaussian errors on c̄.

truthnoisy c̄ data / naive inverseTikhonov estimate

Profile

Penalty L

Noise on c̄

show naive inverse

Try: 100 depth samples, Δz = 0.4 mm. The naive inverse (amber) is usable only in the first few millimetres: its noise grows linearly with depth, about 35 m/s at 5 mm and 210 m/s at 30 mm for 2 m/s on c̄. Press Oracle λ for each L: D₁ keeps the fat/muscle step sharpest among the quadratic penalties, L = I biases deep values toward 1540 m/s, and on the ramp profile D₂ follows the slope that D₁ flattens. None of them keeps a sharp step and low noise at the same time; 13.7 shows the penalty that can.

Ultrasound link

The same algebra sits under any "average to local" step: local attenuation from cumulative attenuation (Module 27), local sound speed from depth-dependent autofocus, local delay from cumulative phase. Each differentiates a running integral, so each needs regularization whose strength should grow with depth.

13.4Choosing λ: discrepancy, GCV, L-curve

In simulation you can scan λ against the truth (the oracle). On data you need a rule. The three standard ones use different information and fail differently:

\[ \begin{gathered} \text{Discrepancy:}\quad \lVert\mathbf{y} - \mathbf{A}\mathbf{x}_\lambda\rVert = \tau\lVert\mathbf{n}\rVert \\ \approx \tau\sqrt{m}\,\sigma_n, \qquad \tau \ge 1 \end{gathered} \]

(13.8)

\(\mathbf{x}_\lambda,\ \mathbf{A},\ \mathbf{y}\): as in (13.3)
\(\mathbf{n}\): noise vector, as in (13.1)
\(\sigma_n\): noise std per sample, units of \(\mathbf{y}\)
\(m\): number of data samples
\(\tau\): safety factor, 1 or slightly above
\(\lVert\cdot\rVert\): Euclidean norm

Needs the noise level. Do not fit the data better than the noise allows. It tends to over-smooth slightly (τ = 1 here; theory uses τ a little above 1).

\[ \begin{gathered} \text{GCV:}\quad G(\lambda) = \frac{m\,\lVert\mathbf{y} - \mathbf{A}\mathbf{x}_\lambda\rVert^2}{\big[m - \operatorname{trace}(\mathbf{A}\mathbf{A}^{\#})\big]^2} \\ \operatorname{trace}(\mathbf{A}\mathbf{A}^{\#}) = \sum_i f_i \end{gathered} \]

(13.9)

\(G(\lambda)\): GCV score, minimized over \(\lambda\)
\(m\): number of data samples
\(\mathbf{x}_\lambda,\ \mathbf{A},\ \mathbf{y}\): as in (13.3)
\(\mathbf{A}^{\#}\): regularized inverse, as in (13.4)
\(\operatorname{trace}\): sum of the diagonal; here the effective number of parameters fitted
\(f_i\): filter factors, as in (13.3)

Leave-one-out prediction error made rotation invariant (Golub, Heath & Wahba, Technometrics, 1979). No \(\sigma_n\) needed, but it assumes white noise; correlated noise looks like signal and GCV then under-smooths badly. Its minimum can also be very flat.

\[ \begin{gathered} \text{L-curve:}\quad \big(\log\lVert\mathbf{y} - \mathbf{A}\mathbf{x}_\lambda\rVert,\ \log\lVert\mathbf{L}\mathbf{x}_\lambda\rVert\big) \\ \text{corner = point of maximum curvature} \end{gathered} \]

(13.10)

\(\mathbf{x}_\lambda,\ \mathbf{A},\ \mathbf{y}\): as in (13.3)
\(\mathbf{L}\): penalty operator, as in (13.6)
\(\log\): natural logarithm; \(\lVert\mathbf{y}-\mathbf{A}\mathbf{x}_\lambda\rVert\) is the residual norm, \(\lVert\mathbf{L}\mathbf{x}_\lambda\rVert\) the solution (penalty) norm
\(\lVert\cdot\rVert\): Euclidean norm

Vertical branch: noise-dominated (λ too small); horizontal branch: over-smoothed (Hansen & O'Leary, SIAM J. Sci. Comput., 1993). Heuristic: it can over-smooth for smooth solutions and large problems (Vogel, Inverse Problems, 1996) and has no clear corner at very low or very correlated noise.

Demo 13.3

L-curve, GCV and the discrepancy principle against the oracle

SimulationSame JS inversion as Demo 13.1, scanned over λ; white or AR(1) noise.

oraclediscrepancy (τ = 1)GCVL-curve cornercurves

SNR

Noise

Try: same RF deconvolution as Demo 13.1. With white noise GCV and the corner land near the oracle and the discrepancy principle sits a little to the right (smoother). Switch to AR(1) noise: GCV jumps to the smallest λ and the error explodes; the L-curve corner often follows. The discrepancy principle, told the true noise norm, still works. At 25 dB, m13_tikhonov_deconv_lcurve.m prints oracle 0.020, GCV 0.020, corner 0.025, discrepancy 0.071 for white noise, and GCV 0.0003 (error 31 times the oracle's) for AR(1).

Ultrasound link

Delay estimates between steered acquisitions are correlated over the correlation-kernel size and the speckle size (Module 11), so the white-noise assumption behind GCV is violated in echo-mode sound-speed tomography. Either whiten the data with an estimated covariance, or use the discrepancy principle with a noise level measured from repeated acquisitions, and report the oracle λ alongside in simulation studies.

13.5Iterative regularization: the iteration count is λ

Large problems (2D or 3D tomography, model-based beamforming) never form \(\mathbf{A}^{T}\mathbf{A}\); they apply \(\mathbf{A}\) and \(\mathbf{A}^{T}\) as functions (ray tracing and back-projection; \(\mathbf{A}^{T}\) is DAS-like, Module 05). Started at \(\mathbf{x}_0 = 0\), the simplest iteration is a gradient step on \(\tfrac12\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2\):

\[ \begin{aligned} &\text{Landweber:} \\ &\mathbf{x}_{k+1} = \mathbf{x}_k + \omega\,\mathbf{A}^{T}(\mathbf{y} - \mathbf{A}\mathbf{x}_k) \\ &0 < \omega < 2/\sigma_1^2 \\ &\Rightarrow\ f_i^{(k)} = 1 - (1 - \omega\sigma_i^2)^k \\ &\phantom{\Rightarrow\ f_i^{(k)}} \approx 1 - e^{-k\omega\sigma_i^2} \end{aligned} \]

(13.11)

\(\omega\): step size; \(\sigma_1\) is the largest singular value
\(k\): iteration count
\(f_i^{(k)}\): filter factor after k steps
\(\mathbf{x}_k\): iterate after k steps, starting from \(\mathbf{x}_0 = 0\)
\(\mathbf{A},\ \mathbf{y}\): as in (13.1)
\(^{T}\): matrix transpose

Half power where \(k\omega\sigma_i^2 = \ln 2\), so k Landweber steps act like Tikhonov with \(\lambda^2 \approx \ln 2/(k\omega)\). With \(\omega = 1/\sigma_1^2\): \(k \approx \ln 2\,(\sigma_1/\lambda)^2\). At 25 dB the oracle \(\lambda = 0.020\,\sigma_1\) predicts 1720 steps; the error minimum is at 1823 (m13_cgls_semiconvergence.m).

\[ \begin{aligned} &\text{CGLS:}\quad \mathbf{x}_k = \arg\min_{\mathbf{x} \in \mathcal{K}_k} \lVert\mathbf{y} - \mathbf{A}\mathbf{x}\rVert \\ &\mathcal{K}_k = \operatorname{span}\{\mathbf{A}^{T}\mathbf{y},\ (\mathbf{A}^{T}\mathbf{A})\mathbf{A}^{T}\mathbf{y}, \\ &\qquad\quad \dots,\ (\mathbf{A}^{T}\mathbf{A})^{k-1}\mathbf{A}^{T}\mathbf{y}\} \end{aligned} \]

(13.12)

\(\mathbf{x}_k\): CGLS iterate after k steps
\(\mathcal{K}_k\): Krylov subspace of dimension k
\(\operatorname{span}\{\cdots\}\): all linear combinations of the listed vectors
\(\mathbf{A},\ \mathbf{y}\): as in (13.1)
\(\arg\min\): the \(\mathbf{x}\) that minimizes the residual norm
\(^{T}\): matrix transpose

Conjugate gradients on the normal equations, implemented with \(\mathbf{A}\) and \(\mathbf{A}^{T}\) only. Its filter is a data-adapted polynomial, so it reaches the same error in far fewer steps: 35 instead of 1823 here, roughly the square root (polynomial acceleration; Engl, Hanke & Neubauer 1996).

Both show semiconvergence: the error falls while the large-σ (signal) components converge, then rises as the small-σ (noise) components follow. Stop with the discrepancy principle on the residual \(\lVert\mathbf{r}_k\rVert \le \tau\sqrt{m}\,\sigma_n\).

Demo 13.4

Semiconvergence of CGLS and Landweber

SimulationCGLS iterated in the page (JS) up to k = 400; Landweber from its closed-form filter factors (13.11).

CGLSLandweber (ω = 1/σ₁²)Tikhonov oracle error / band-limited truthdiscrepancy stop

Show iterate of

Iteration k

SNR

Try: slide k past the minimum of the CGLS curve and watch the iterate fill with noise; that is "more converged" and worse. Lower the SNR: both minima move to fewer iterations, as a larger λ would. CGLS iterates are computed up to k = 400; the Landweber curve uses its closed-form filter factors.

Pitfall

Papers that report "we ran 50 iterations of CG" without a stopping rule have chosen a regularization parameter silently. When you review or reproduce one, ask for the stopping rule, and when you write one, report k with the rule that chose it.

13.6Sparsity, L1, and the Bayesian reading

When the object is a few point reflectors (or bubbles), penalize the number of nonzeros. The convex stand-in is the L1 norm. To keep symbols apart, μ weights L1 and TV penalties and λ stays the Tikhonov parameter:

\[ \begin{aligned} &\min_{\mathbf{x}}\ F(\mathbf{x}) = \tfrac12\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2 + \mu\lVert\mathbf{x}\rVert_1 \\ &\operatorname{prox}_{\tau\lvert\cdot\rvert}(v) = S_\tau(v) \\ &\phantom{\operatorname{prox}_{\tau\lvert\cdot\rvert}(v)} = \operatorname{sign}(v)\cdot\max(\lvert v\rvert - \tau,\ 0) \end{aligned} \]

(13.13)

\(\mu\): L1 weight, units of \(\mathbf{A}^{T}\mathbf{y}\)
\(\lVert\mathbf{x}\rVert_1\): \(\sum_j \lvert x_j\rvert\)
\(S_\tau\): soft-thresholding operator with threshold \(\tau\)
\(F(\mathbf{x})\): objective (data misfit plus L1 penalty)
\(\mathbf{A},\ \mathbf{x},\ \mathbf{y}\): as in (13.1)
\(\operatorname{prox}_{\tau\lvert\cdot\rvert}(v)\): proximal operator of \(\tau\lvert\cdot\rvert\), applied to each element
\(v\): input scalar (element of a vector)
\(\tau\): threshold, same units as \(v\)
\(\operatorname{sign}\): +1, 0 or \(-1\) for positive, zero or negative argument
\(\max(\cdot,\cdot)\): larger of the two arguments

Soft thresholding is the proximal map of the L1 norm: the exact minimizer of \(\tfrac12(x - v)^2 + \tau\lvert x\rvert\). It sets \(\lvert v\rvert \le \tau\) to zero and shrinks the rest by \(\tau\).

\[ \begin{aligned} &\text{ISTA:} \\ &\mathbf{x}_{k+1} = S_{\mu/L}\!\left(\mathbf{x}_k - \mathbf{A}^{T}(\mathbf{A}\mathbf{x}_k - \mathbf{y})/L\right) \\ &L = \lVert\mathbf{A}\rVert^2 = \sigma_1^2 \\[4pt] &\text{FISTA: same step at} \\ &\mathbf{z}_k = \mathbf{x}_k + \frac{t_{k-1} - 1}{t_k}\,(\mathbf{x}_k - \mathbf{x}_{k-1}) \\ &t_k = \frac{1 + \sqrt{1 + 4t_{k-1}^2}}{2} \end{aligned} \]

(13.14)

\(\mathbf{x}_k\): ISTA iterate after k steps
\(S_{\mu/L}\): soft threshold at \(\mu/L\), as in (13.13)
\(\mu\): L1 weight, as in (13.13)
\(\mathbf{A},\ \mathbf{y}\): as in (13.1)
\(L\): Lipschitz constant of the gradient, \(\lVert\mathbf{A}\rVert^2 = \sigma_1^2\) (step size is \(1/L\))
\(\lVert\mathbf{A}\rVert\): spectral norm, equal to the largest singular value \(\sigma_1\)
\(\mathbf{z}_k\): FISTA extrapolated point where the step is taken
\(t_k\): FISTA momentum sequence, \(t_0 = 1\)

Worst-case objective gaps from \(\mathbf{x}_0\): \(F - F^* \le L\lVert\mathbf{x}_0 - \mathbf{x}^*\rVert^2/(2k)\) for ISTA and \(\le 2L\lVert\mathbf{x}_0 - \mathbf{x}^*\rVert^2/(k + 1)^2\) for FISTA (Beck & Teboulle 2009; ISTA's convergence itself was proved by Daubechies, Defrise & De Mol 2004). A starting μ for noise \(\sigma_n\) and equal column norms \(\lVert\mathbf{a}\rVert\) is the universal threshold \(\mu \approx \sigma_n\lVert\mathbf{a}\rVert\sqrt{2\ln N}\) (exact rationale for orthonormal A, Donoho & Johnstone 1994).

L1 is biased. At an isolated spike the optimality condition gives \(x_j = (\mathbf{a}_j^{T}\mathbf{y} - \mu\,\operatorname{sign} x_j)/\lVert\mathbf{a}_j\rVert^2\): every surviving amplitude is pulled toward zero by \(\mu/\lVert\mathbf{a}_j\rVert^2\). Use L1 to find the support, then refit amplitudes by least squares on that support (debiasing).

The Bayesian reading. With Gaussian noise, the maximum a posteriori estimate minimizes the negative log-likelihood plus the negative log-prior (Module 12):

\[ \begin{gathered} \mathbf{x}_{\mathrm{MAP}} = \arg\min_{\mathbf{x}}\ \frac{\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2}{2\sigma_n^2} - \log p(\mathbf{x}) \\ \mathcal{N}(0, \sigma_x^2\mathbf{I}):\ \ \lambda^2 = \frac{\sigma_n^2}{\sigma_x^2} \\ \text{Laplace } p \propto e^{-\lvert x\rvert/\beta}:\ \ \mu = \frac{\sigma_n^2}{\beta} \end{gathered} \]

(13.15)

\(p(\mathbf{x})\): prior density of the unknowns
\(\sigma_x\): prior std per unknown (Gaussian prior)
\(\beta\): scale of the Laplace prior, same units as \(x\)
\(\sigma_n\): noise std

Tikhonov = Gaussian prior; Tikhonov with \(\mathbf{D}_1\) = Gaussian prior on increments (a random walk); L1 = Laplace prior; TV = Laplace prior on increments. Sparsity is a property of the MAP estimator, not of the prior: a Laplace draw is never exactly zero, and only \(1 - e^{-0.1}\) = 9.5 % of its samples have \(\lvert x\rvert < 0.1\beta\).

Demo 13.5

Sparse-spike deconvolution: L2 vs L1, ISTA vs FISTA

SimulationISTA and FISTA iterated in the page (JS) on a 256-sample convolution problem; Tikhonov by FFT.

true spikesTikhonov (oracle λ)L1 estimate / FISTAISTA

μ (× universal threshold)

Algorithm

Iterations

SNR

debias (LS on support)

Try: 256 samples, same 5 MHz pulse (−6 dB envelope 294 ns = 11.8 samples). Seven spikes include an opposite-sign pair 6 samples apart (at 70 and 76) and a same-sign pair 4 samples apart (at 170 and 174: 0.1 µs, half a period, λ/4 in depth). L1 resolves the first pair, which Tikhonov smears. The second pair is lost by both: its two echoes nearly cancel in band (spectrum 1 − 0.9 at f₀), so the data barely contain it. Switch to ISTA at 100 iterations, then FISTA: same fixed point, very different speed.

Ultrasound link

ULM at high bubble density, compressed-sensing beamforming and sparse-spike deconvolution all solve this L1 problem with a PSF dictionary. The demo's lesson transfers: L1 resolves sources closer than the PSF only when their echoes remain distinguishable in band, it biases amplitudes (debias before quantifying intensity), and its answer depends on μ, which should be set from the noise.

13.7Total variation: piecewise-constant priors

Tissue sound speed is closer to piecewise constant (fat, muscle, liver, a lesion) than to smooth. Total variation penalizes the L1 norm of the gradient, so it allows a few large jumps and suppresses many small ones (Rudin, Osher & Fatemi, Physica D, 1992):

\[ \begin{aligned} &\min_{\mathbf{x}}\ \tfrac12\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2 + \mu\,\mathrm{TV}(\mathbf{x}) \\ &\mathrm{TV}_{\mathrm{1D}} = \sum_i \lvert x_{i+1} - x_i\rvert \\ &\mathrm{TV}_{\mathrm{2D,aniso}} = \sum \big(\lvert\partial_x x\rvert + \lvert\partial_z x\rvert\big) \\ &\mathrm{TV}_{\mathrm{2D,iso}} = \sum \sqrt{(\partial_x x)^2 + (\partial_z x)^2} \end{aligned} \]

(13.16)

\(\mathbf{A},\ \mathbf{x},\ \mathbf{y}\): as in (13.1)
\(\mu\): TV weight, in units of \(\mathbf{A}^{T}\mathbf{y}\)
\(\mathrm{TV}(\mathbf{x})\): total variation: sum of the magnitudes of the differences (gradient) of \(\mathbf{x}\)
\(x_i\): unknown at sample i
\(\partial_x,\ \partial_z\): finite differences along lateral (x) and axial (z) pixel index
\(\lvert\cdot\rvert\): absolute value

Anisotropic TV favours axis-aligned edges; isotropic TV is rotation invariant. Smoothed TV, \(\sum\sqrt{(\mathbf{D}\mathbf{x})^2 + \varepsilon^2}\), is differentiable and can be minimized with gradient or lagged-diffusivity iterations.

\[ \begin{aligned} &\text{ADMM, split } \mathbf{z} = \mathbf{D}\mathbf{x}: \\ &\mathbf{x} \leftarrow (\mathbf{A}^{T}\mathbf{A} + \rho\,\mathbf{D}^{T}\mathbf{D})^{-1} \\ &\qquad\times\big(\mathbf{A}^{T}\mathbf{y} + \rho\,\mathbf{D}^{T}(\mathbf{z} - \mathbf{u})\big) \\ &\mathbf{z} \leftarrow S_{\mu/\rho}(\mathbf{D}\mathbf{x} + \mathbf{u}) \\ &\mathbf{u} \leftarrow \mathbf{u} + \mathbf{D}\mathbf{x} - \mathbf{z} \end{aligned} \]

(13.17)

\(\mathbf{D}\): difference (gradient) operator
\(\rho\): ADMM penalty parameter
\(\mathbf{u}\): scaled dual variable
\(\mathbf{x}\): unknown image/vector, as in (13.1)
\(\mathbf{A},\ \mathbf{y},\ \mu\): as in (13.1) and (13.13)
\(\mathbf{z}\): auxiliary copy of \(\mathbf{D}\mathbf{x}\), the variable that carries the L1 penalty
\(S_{\mu/\rho}\): soft threshold at \(\mu/\rho\), as in (13.13)
\(\leftarrow\): assignment: update the left variable with the right side
\(^{T},\ ^{-1}\): matrix transpose, matrix inverse

One factorization of \(\mathbf{A}^{T}\mathbf{A} + \rho\mathbf{D}^{T}\mathbf{D}\) serves every iteration (Boyd et al., Found. Trends Mach. Learn., 2011). ρ changes the speed, not the answer: on the 1D sound-speed problem ρ = 0.1 converges in about 1000 iterations, ρ = 2 is still 1 % off in objective after 3000.

Two signatures to recognise. TV keeps edges but loses contrast: denoising an isolated plateau of width w samples lowers it by exactly \(2\mu/w\) (m13_tv_vs_tikhonov.m: μ = 1, w = 20 gives 0.900). And it turns ramps into staircases, because a staircase has the same total variation as the ramp and fits the data better locally.

Demo 13.6

TV vs Tikhonov-D₁ on the average-to-local sound-speed problem

SimulationApprox. modelVertical-path travel-time average, independent Gaussian errors on c̄; Tikhonov and ADMM-TV solved in the page (JS).

truthTikhonov, L = D₁TV (ADMM, ρ = 0.1, 1500 iterations)

Profile

Tikhonov λ

TV μ

Noise on c̄

Try: Oracle on the layers: TV cuts the RMSE by a factor of 2 to 3 and keeps the 10–90 % width of the 8 mm step within about 1 mm (often one 0.4 mm sample), against about 3 mm for Tikhonov. On the ramp profile compare the readout for z > 24 mm: TV staircases there and most of its advantage disappears (m13_tv_vs_tikhonov.m finds TV and D₂ tied at 7.2 m/s on that segment). On the thin layer raise μ and watch the plateau sink before it disappears, the contrast-loss signature.

Pitfall

A TV reconstruction that shows crisp, flat regions looks convincing whether or not the tissue is piecewise constant: the prior draws the edges. Check the staircase test (does a known gradient come back as steps?) and report values averaged over regions well inside each plateau, where contrast loss is smallest.

13.8Echo-mode sound-speed tomography with limited angles

The measurement. Beamform the same region with an assumed \(c_0\) from several transmit angles \(\theta_t\) and receive angles \(\theta_r\). If the true slowness differs by \(\Delta s\), the echo from pixel \(\mathbf{r}\) is misaligned by the travel-time error along the transmit path plus the receive path. In the straight-ray, small-perturbation limit each path contributes a line integral from \(\mathbf{r}\) back to the array:

\[ \begin{gathered} \tau(\mathbf{r};\, \theta_t, \theta_r) = L_{\theta_t}(\mathbf{r}) + L_{\theta_r}(\mathbf{r}) \\ L_\theta(\mathbf{r}) = \int_{\text{path from } \mathbf{r} \text{ at } \theta} \Delta s\ dl \\ \Delta s = \frac{1}{c} - \frac{1}{c_0} \end{gathered} \]

(13.18)

\(\mathbf{r}\): pixel position \((x, z)\)
\(\theta_t,\ \theta_r\): transmit and receive steering angles from vertical
\(\Delta s\): slowness deviation, ns/mm (1 ns/mm = 10⁻⁶ s/m)
\(c_0\): beamforming sound speed, m/s
\(\tau\): delay error of the echo, ns
\(L_\theta(\mathbf{r})\): delay accumulated along the one-way path through \(\mathbf{r}\) at angle \(\theta\), ns
\(\theta\): ray angle from vertical, deg
\(L_{\theta_t},\ L_{\theta_r}\): transmit-path and receive-path delay errors, ns
\(dl\): path length element along the ray, mm
\(c\): true local sound speed, m/s

Linear in slowness, not in c: \(\Delta c \approx -c_0^2\,\Delta s\), so 1 ns/mm of slowness is 2.37 m/s at 1540 m/s, and a +50 m/s inclusion is −20.4 ns/mm (10 mm of it advances an echo by 204 ns, a full period at 5 MHz). Refraction and the position shift of the misfocused echo are neglected.

The absolute τ is not observable, because the echo phase of a speckle pixel is random. What is measured, as a phase shift or delay between beamformed images (Module 11), is the difference between angle pairs whose speckle stays correlated. Echoes from pairs with the same mid-angle (\(\theta_t + \theta_r = 2\theta_m\)) sample the same region of k-space and stay correlated, so a common-midangle difference is a natural datum. With K angles in \([-\theta_{\max}, \theta_{\max}]\):

\[ \begin{gathered} d_{ij}(\mathbf{r}) = L_i(\mathbf{r}) + L_j(\mathbf{r}) - 2L_m(\mathbf{r}) \\ m = (i + j)/2 \\ \Rightarrow\quad \mathbf{d} = \mathbf{A}\,\Delta\mathbf{s} + \mathbf{n} \end{gathered} \]

(13.19)

\(d_{ij}(\mathbf{r})\): common-midangle delay difference at pixel \(\mathbf{r}\) for angle pair (i, j), ns
\(L_i,\ L_j,\ L_m\): one-way path delay errors along the rays of angles i, j and m, as in (13.18), ns
\(i,\ j,\ m\): angle indices; m is the mid-angle of the pair
\(\mathbf{d}\): stacked delay differences of all pairs and pixels
\(\mathbf{A}\): ray-weight matrix: each row is a difference of three ray-length vectors, mm
\(\Delta\mathbf{s}\): slowness deviation map, ns/mm, as in (13.18)
\(\mathbf{n}\): noise vector, ns

Each row of \(\mathbf{A}\) is a difference of three ray-weight vectors. Computed ultrasound tomography in echo mode (CUTE; Jaeger et al., Ultrasound Med. Biol., 2015) measured echo phase shifts between transmit-steered images with a fixed receive aperture; the common-midangle pairing of transmit and receive angles used here was introduced with the improved forward model of Stähli et al. (Ultrasonics, 2020). Spatial-domain formulations were developed in parallel (Sanabria et al., Phys. Med. Biol., 2018). The toy here keeps their linear-algebra skeleton, not their exact acquisition schemes.

Why depth is smeared. By the Fourier-slice theorem, line integrals along a direction measure the spatial frequencies perpendicular to it. Rays within \(\pm\theta_{\max}\) of vertical therefore sample only a wedge of k-space of half-angle \(\theta_{\max}\) around the \(k_x\) axis. The missing cone around \(k_z\) means horizontal edges and axial detail are poorly determined: inclusions stretch in depth, and the axial width of the resolution kernel grows roughly as \(1/\sin\theta_{\max}\) times the lateral width. Differencing makes layered media nearly invisible as well. For \(\Delta s(z)\) only, every ray from depth z collects \(L_\theta = S(z)/\cos\theta\), where \(S(z)\) is the vertical one-way delay, and the symmetric pair gives

\[ \begin{aligned} d &= L_{-\theta} + L_{\theta} - 2L_0 \\ &= 2(\sec\theta - 1)\,S(z) \approx \theta^2 S(z) \end{aligned} \]

(13.20)

\(d\): delay difference of the symmetric pair, ns
\(L_{-\theta},\ L_{\theta},\ L_0\): path delay errors at angles \(-\theta\), \(+\theta\) and vertical, ns
\(\theta\): ray angle from vertical, rad
\(S(z)\): vertical one-way delay error accumulated to depth z, ns
\(\sec\theta\): \(1/\cos\theta\), path-length stretch of a slanted ray

Only 0.031, 0.128 and 0.309 of the vertical delay \(S(z)\) at 10°, 20° and 30°: second order in θ.

SchematicDrawing for intuition. Array, grid and depth range to scale (8 px/mm); straight rays; the k-space wedge is drawn for θmax = 20°.

Common-midangle datum for one pixel: the transmit ray (−20°) and receive ray (+20°) collect slowness along their paths, the 0° reference with the same mid-angle is subtracted twice. All rays from a linear array lie within ±θmax of vertical, so their line integrals sample only the shaded k-space wedge around kₓ.

What the inversion can see. With a gradient penalty \(\mathbf{G}\) (Hansen convention), \(\mathbf{R} = (\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{G}^{T}\mathbf{G})^{-1}\mathbf{A}^{T}\mathbf{A}\). Row k of \(\mathbf{R}\) is the point-spread function of the reconstruction at pixel k (how its estimate mixes the true values around it), and \(\operatorname{trace}\mathbf{R}\) counts the parameters actually determined. m13_sos_tomography_limited_angle.m (two inclusions 12 mm apart in depth, 5 ns delay noise, oracle λ) gives:

\(\theta_{\max}\)	data rows	R row FWHM axial (mm)	lateral (mm)	trace R / 480	RMSE Tikhonov (m/s)	RMSE TV (m/s)
10°	3908	9.4	1.2	66	7.7	7.6
20°	3462	4.6	1.3	132	6.7	5.6
30°	2982	3.3	1.3	185	4.8	2.4

The lateral width is pinned near the 1.6 mm pixel. Wider angles give fewer rows (more rays leave the aperture) but much better conditioning. TV helps a lot at 30° and not at all at 10°: a prior can only sharpen what the data partly constrain.

Demo 13.7

Limited-angle echo-mode sound-speed tomography: reconstruction, resolution row, missing cone

SimulationApprox. modelStraight rays, linearized in slowness, no refraction; common-midangle differences of ray integrals.

transmit ray −θ_maxreceive ray +θ_maxcommon-midangle reference 0°other angles / truth profile

Truth

Angle span ±θ_max (7 angles)

Delay noise

Regularizer

Tikhonov λ (× σ_max)

TV μ

Try: 24 × 20 pixels of 1.6 mm under a 38.4 mm array (128 × 0.3 mm), z = 3–35 mm, K = 7 angles, every pixel a measurement point. Data come from 0.1 mm ray steps through a smooth-edged truth, the inversion uses bilinear 0.4 mm steps on the grid (no inverse crime). Start with the stacked pair at 10°: one vertical streak, and a resolution row many pixels tall. Step to 30° and they separate. Click any map to move the resolution row; its 2D spectrum (the |FFT₂| panel) fills the wedge between the dashed \(\pm\theta_{\max}\) lines around \(k_x\). Choose the fat layer: it is sensed only weakly and leaks into streaks. Switch to TV and press Oracle parameter at 30°, then at 10°.

Ultrasound link

When you read an echo-mode sound-speed paper, look for four things: the angle span (it sets the missing cone), what is differenced (it sets the blind spots, such as layers), the regularizer and how its weight was chosen, and whether accuracy is reported for inclusions large compared with the axial resolution. Mean values inside a 10 mm inclusion can be accurate while a 3 mm inclusion is smeared to half its contrast; the resolution row tells you which regime you are in.

Pitfall

The inverse crime: simulating data with the same discretization you invert with. The data then lie exactly in the range of A, and reconstructions look far better than they will on phantoms. Generate data with a finer grid, a different ray model (or a full-wave simulator), and noise with realistic correlation.

13.9Visual summary

Twelve questions an inverse-problems reviewer answers in five seconds, each with the picture behind the answer. Every curve is computed from the formula in the card.

13.10How experts think about it

An expert looks at an inversion as three objects: the singular spectrum of the forward operator (what the data can determine), the prior encoded by the regularizer (what fills in the rest), and the resolution matrix that results (what you actually see). The regularization weight is a statement about noise level and prior strength, so it is chosen from the noise, reported with its convention, and checked against the oracle in simulation.

Ill-posed means small singular values, not a singular matrix: the inverse exists and is useless.

noise gain \(\sigma_n/\sigma_i\)\(\lvert H\rvert < 0.01\) outside ≈ 1.1–8.9 MHz\(1/\lvert H\rvert\) ≈ 130 dB at 12 MHz

The naive solution is \(\mathbf{x} + \sum(\mathbf{u}_i^{T}\mathbf{n}/\sigma_i)\,\mathbf{v}_i\) (13.2). For a 5 MHz, 60 % bandwidth Gaussian pulse \(\sigma = \lvert H(f)\rvert = \exp(-(f - f_0)^2/2\sigma_f^2)\) with \(\sigma_f\) = 1.27 MHz, so \(\lvert H\rvert = 0.01\) at \(f_0 \pm 3.87\) MHz and \(\lvert H(12\ \mathrm{MHz})\rvert = e^{-15.1}\), −131 dB. Any noise there is amplified by \(10^{6.5}\). A matrix can be numerically invertible (κ = 310 for the running mean) and still multiply noise by 100.

Approx. modelGaussian pulse spectrum, 5 MHz, 60 % −6 dB bandwidth.

κ alone hides where the small σ are; the Picard plot and the singular vectors tell you which features are lost.

λ is a soft threshold on the singular values: components with \(\sigma \gg \lambda\) pass, \(\sigma \ll \lambda\) are damped as \(\sigma^2/\lambda^2\).

\(f = \sigma^2/(\sigma^2 + \lambda^2)\)\(f = \tfrac12\) at \(\sigma = \lambda\)Module 05's λ = our \(\lambda^2\)

Tikhonov's normal equations diagonalize in the SVD basis: \(x_i = \sigma_i(\mathbf{u}_i^{T}\mathbf{y})/(\sigma_i^2 + \lambda^2)\) (13.3). With σ = (1, 0.1, 0.01) and λ = 0.03: f = 0.999, 0.917, 0.100. Writing the penalty as \(\lambda^2\lVert\mathbf{x}\rVert^2\) gives λ the units of σ, so it can be compared directly with \(\lvert H(f)\rvert\) or with the noise floor of a Picard plot.

Exact formula

Quoted values differ by a square between conventions; with L ≠ I the comparison is with generalized singular values, not σ.

Tikhonov never amplifies noise by more than \(1/(2\lambda)\).

\(\max\ \sigma/(\sigma^2 + \lambda^2) = 1/(2\lambda)\)λ = 0.02 → 25×

The gain from \(\mathbf{u}_i^{T}\mathbf{n}\) to the solution is \(f_i/\sigma_i = \sigma_i/(\sigma_i^2 + \lambda^2)\) (13.4). Its derivative vanishes at \(\sigma = \lambda\), where it equals \(1/(2\lambda)\). The worst components are those right at the threshold, not the smallest ones (which are suppressed). With \(\lambda = 0.02\,\sigma_1\), no noise component grows more than 25 times, versus \(1/\sigma_{\min}\) for the naive inverse.

Exact formula

Per component; the total noise is \(\sigma_n\big(\sum f_i^2/\sigma_i^2\big)^{1/2}\), which depends on how many \(\sigma_i\) sit near λ.

A regularized inversion returns the truth seen through \(\mathbf{R}\): for deconvolution, a new pulse whose spectrum is the filter factor.

\(\mathbf{R} = \mathbf{V}\operatorname{diag}(f)\mathbf{V}^{T}\)bias \(= (\mathbf{R} - \mathbf{I})\mathbf{x}\)\(\operatorname{trace}\mathbf{R} = \sum f_i\) = dof

\(\mathrm{E}[\mathbf{x}_\lambda] = \mathbf{A}^{\#}\mathbf{A}\mathbf{x} = \mathbf{R}\mathbf{x}\). For a circulant \(\mathbf{A}\), \(\mathbf{R}\) is a convolution with \(r(t) = \mathcal{F}^{-1}\{\lvert H\rvert^2/(\lvert H\rvert^2 + \lambda^2)\}\), a band-pass kernel whose band ends where \(\lvert H\rvert = \lambda\) (13.5). Its width is the axial resolution of the deconvolved line; \(\operatorname{trace}\mathbf{R}\) counts the effective number of parameters the data determine (Demo 13.1 shows it).

Approx. modelGaussian pulse spectrum, 5 MHz, 60 % −6 dB bandwidth.

Exact for linear estimators only; L1 and TV have data-dependent resolution, which you probe with point-target simulations instead.

Choose L to say what you believe: I pulls toward the reference, D₁ toward flat, D₂ toward straight lines.

\(\mathcal{N}(\mathbf{D}_1)\) = constants\(\mathcal{N}(\mathbf{D}_2)\) = linesneed \(\mathcal{N}(\mathbf{A}) \cap \mathcal{N}(\mathbf{L}) = \{0\}\)

The penalty is blind to its null space, so those components are set by the data alone (13.6). On the layered sound-speed profile with 2 m/s noise, oracle Tikhonov gives a 10–90 % fat/muscle edge of 6.4 mm with \(\mathbf{L} = \mathbf{I}\) and 2.8 mm with \(\mathbf{D}_1\); on the ramp segment \(\mathbf{D}_2\) gives 7.2 m/s RMSE against 11.4 for \(\mathbf{D}_1\) (m13_tv_vs_tikhonov.m).

Exact formulaPenalty spectra \(\lvert L(k)\rvert^2\) of I, D₁, D₂ on a uniform grid.

If A also cannot see the null space of L (for example a constant offset invisible to differential data), the solution is not unique and the solver returns an arbitrary member.

Differentiating a running average is the everyday worst case: the noise gain grows as \(\sqrt{2}\,z/\Delta z\).

1 m/s → 110 m/s at 30 mmΔz = 0.4 mmκ = 310

\(s_i = (z_i\bar s_i - z_{i-1}\bar s_{i-1})/\Delta z\), so independent errors in \(\bar s\) add with weights \(z/\Delta z\): \(\operatorname{std}(s) \approx \sqrt{2}\,(z/\Delta z)\operatorname{std}(\bar s)\), and in m/s the factor \((c/\bar c)^2\) (1.03–1.15 for the layered profile) multiplies it (13.7). Simulated: 39, 77, 111 m/s at 10, 20, 30 mm for 1 m/s on \(\bar c\), against 40, 75, 109 from the formula with that factor (the script prints both).

Approx. modelCurve \(\sqrt{2}\,z/\Delta z\) (independent errors per depth sample, \((c/\bar c)^2\) set to 1); dots from m13_tv_vs_tikhonov.m.

Assumes independent errors per depth sample; real autofocus estimates are correlated over the estimation window, which lowers the gain for differences over shorter spans and makes GCV unreliable.

The discrepancy principle needs σ and over-smooths, GCV needs white noise, the L-curve needs a corner.

\(\lVert\mathbf{r}\rVert = \tau\sqrt{m}\,\sigma_n\)\(G = m\lVert\mathbf{r}\rVert^2/(m - \sum f_i)^2\)AR(1): GCV error 31× oracle

On the RF deconvolution at 25 dB with white noise, the oracle is λ = 0.020, GCV 0.020, the corner 0.025, the discrepancy principle 0.071. With AR(1) noise (ρ = 0.9) GCV picks 0.0003 and the relative error is 11.8 against 0.38 for the oracle, because correlated noise is predictable and GCV mistakes it for signal.

Approx. modelExpected norms for the Gaussian-pulse spectrum, white reflectivity, \(\sigma_n\) = 0.02.

All three can fail on the same data; in simulation always report the oracle too, and on data check that the chosen λ is stable to a new noise draw.

Iteration count is a regularization parameter: Landweber needs \(k \approx \ln 2\,(\sigma_1/\lambda)^2\), CGLS about the square root of that.

\(f^{(k)} = 1 - (1 - \omega\sigma^2)^k\)1720 predicted, 1823 foundCGLS: 35

With \(\omega = 1/\sigma_1^2\), \(1 - (1 - \sigma^2)^k \approx 1 - e^{-k\sigma^2}\) reaches ½ at \(\sigma^2 = \ln 2/k\) (13.11). The 25 dB oracle λ = 0.020 predicts k = 1720; the Landweber error minimum is at 1823 with the same error (0.33). CGLS reaches it in 35 steps because its Krylov polynomial adapts to the spectrum.

Approx. modelSame model problem as the L-curve; Landweber filter in closed form.

The CGLS count depends on the singular-value distribution and on the data; without a stopping rule (discrepancy, or a validation set) the iteration count is an unreported tuning knob.

Soft thresholding is the proximal map of L1, and it shrinks every surviving amplitude by \(\mu/\lVert\mathbf{a}\rVert^2\).

\(S_\tau(v) = \operatorname{sign} v\cdot\max(\lvert v\rvert - \tau, 0)\)shrink 0.02 in Demo 13.5

At a nonzero \(x_j\), the optimality condition \(\mathbf{a}_j^{T}(\mathbf{y} - \mathbf{A}\mathbf{x}) = \mu\,\operatorname{sign} x_j\) gives \(x_j = (\mathbf{a}_j^{T}\mathbf{y} - \mu\,\operatorname{sign} x_j)/\lVert\mathbf{a}_j\rVert^2\) for an isolated spike. With \(\lVert h\rVert\) = 0.336 and μ = 0.0023 the shrink is 0.020; in m13_ista_sparse_spikes.m the isolated spikes come back 0.02–0.03 low (the unit spike at 0.97), because noise and pulse overlap with neighbouring estimates add to the single-spike formula. Least squares on the support removes it (debiasing).

Exact formula

Neighbouring spikes interact through \(\mathbf{A}^{T}\mathbf{A}\), so the shrink is only approximately \(\mu/\lVert\mathbf{a}\rVert^2\); debiasing amplifies noise if the support has nearly parallel columns.

ISTA converges as \(1/k\) and FISTA as \(1/k^2\), both with step \(1/\lVert\mathbf{A}\rVert^2\).

\(L\lVert\mathbf{x}_0 - \mathbf{x}^*\rVert^2/(2k)\)\(2L\lVert\mathbf{x}_0 - \mathbf{x}^*\rVert^2/(k + 1)^2\)FISTA \(10^{-6}\) gap in 274 it

The bounds are worst-case (Beck & Teboulle 2009); both hold in m13_ista_sparse_spikes.m. There the relative gap after 1000 iterations is \(6.5\times10^{-3}\) for ISTA and \(10^{-11}\) for FISTA, which crosses \(10^{-6}\) at iteration 274. Matching FISTA's bound at 100 iterations would take ISTA's bound 2550.

Exact formula

FISTA's objective is not monotone (it ripples); with strong convexity ISTA itself is linear, and restarted FISTA is the practical default.

Start μ at the universal threshold \(\sigma_n\lVert\mathbf{a}\rVert\sqrt{2\ln N}\), then adjust.

N = 256: \(3.33\,\sigma_n\lVert\mathbf{a}\rVert\)μ = 0.0023 in Demo 13.5

For orthonormal \(\mathbf{A}\) the noise coefficients are N independent \(\mathcal{N}(0, \sigma_n^2)\) values and their maximum is close to \(\sigma_n\sqrt{2\ln N}\), so thresholding there removes pure-noise coefficients with high probability (Donoho & Johnstone 1994). For a convolution dictionary the columns have norm \(\lVert h\rVert\), which scales the threshold.

Exact formula

Correlated columns (a band-limited pulse) and correlated noise break the derivation: expect a few false spikes at the universal value (Demo 13.5 shows them) and tune with the discrepancy principle.

No prior recovers what the data cancel: same-sign reflectors λ/4 apart nearly vanish in band.

1 − 0.9 = 0.1 at f₀77 µm at 5 MHz4 samples at 40 MHz

Two echoes half a period apart (0.1 µs, \(c\Delta t/2\) = 77 µm = λ/4 in depth) have spectrum \(1 + 0.9\,e^{-j\pi f/f_0}\), which is 0.1 at \(f_0\): the pair is ten times weaker in band than either spike. L1 cannot recover it; it returns spikes of the wrong sign in the wrong places (+0.40 and +0.36 at 166 and 178 in m13_ista_sparse_spikes.m).

Exact formula

Near-cancellation depends on spacing and amplitude ratio; opposite-sign pairs at the same spacing reinforce instead and are resolved.

Tikhonov is a Gaussian prior, L1 a Laplace prior, TV a Laplace prior on gradients; sparsity belongs to the MAP estimator, not to the prior.

\(\lambda^2 = \sigma_n^2/\sigma_x^2\)\(\mu = \sigma_n^2/\beta\)\(P(\lvert x\rvert < 0.1\beta)\) = 9.5 %

MAP minimizes \(\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2/(2\sigma_n^2) - \log p(\mathbf{x})\) (13.15). A Gaussian prior gives \(\lambda^2 = \sigma_n^2/\sigma_x^2\), the Wiener ratio; a Laplace prior of scale β gives \(\mu = \sigma_n^2/\beta\). Laplace samples are dense (only 9.5 % fall within 0.1β of zero), and the posterior mean under that prior is not sparse either: the zeros come from taking the mode.

Exact formula

The MAP reading justifies the weights only if the prior is right; for uncertainty, the posterior spread matters, which MAP does not provide.

TV keeps edges and pays in contrast: an isolated plateau of width w drops by \(2\mu/w\), and ramps turn into stairs.

\(1 - 2\mu/w\)μ = 1, w = 20: 0.900

For \(\tfrac12\lVert\mathbf{x} - \mathbf{y}\rVert^2 + \mu\sum\lvert x_{i+1} - x_i\rvert\), lowering a plateau by δ reduces the penalty by \(2\mu\delta\) (two edges) and costs about \(w\delta^2/2\) in fit; the optimum \(\delta = 2\mu/w\). Verified to three digits for μ = 0.5–4. A staircase has the same total variation as the ramp it replaces, so TV has no reason to prefer the ramp.

Exact formula

Exact for 1D denoising of an isolated plateau with long flat surroundings; with a blurring A or in 2D (perimeter-to-area ratio) the loss differs but keeps the same trend.

With rays within ±θ of vertical, k-space is a wedge: axial resolution is worse than lateral by roughly \(1/\sin\theta_{\max}\).

\(\mathrm{FWHM}_z\) 9.4 / 4.6 / 3.3 mm\(\mathrm{FWHM}_x\) ≈ 1.3 mm10° / 20° / 30°

Fourier slice: a projection along direction θ measures k perpendicular to it, at angle θ from \(k_x\). Wedge half-angle \(\theta_{\max}\) with radial extent K gives lateral detail ≈ 1/K and axial ≈ \(1/(K\sin\theta_{\max})\). In the simulation the ratio is 7.8, 3.5, 2.5 for \(1/\sin\theta_{\max}\) = 5.8, 2.9, 2.0, with the lateral width limited by the 1.6 mm pixel. trace R rises from 66 to 185 of 480.

Schematic

The ratio also grows with λ; curved rays, real angle weighting and correlated delay noise change the constant.

Differential angle data are nearly blind to layers: their sensitivity is \(2(\sec\theta - 1) \approx \theta^2\).

0.031 at 10°0.128 at 20°0.309 at 30°

For \(\Delta s(z)\) only, every ray from depth z collects \(S(z)/\cos\theta\). The common-midangle difference \(L_{-\theta} + L_\theta - 2L_0 = 2(\sec\theta - 1)\,S(z)\), second order in θ (13.20). A 6 mm fat layer (1450 m/s) delays each vertical path by 242 ns but changes the ±20° datum by only 31 ns, and a ±10° pair (sensitivity 0.031) by 7.5 ns, close to the 5 ns noise of a single pixel.

Exact formula

For non-symmetric pairs the coefficient is \(\sec\theta_i + \sec\theta_j - 2\sec\theta_m\), still second order; absolute-delay (autofocus) data see layers directly, which is why combined approaches are attractive.

Invert slowness, not speed, and remember 1 ns/mm ≈ 2.4 m/s.

\(\Delta c \approx -c_0^2\,\Delta s\)+50 m/s = −20.4 ns/mm10 mm → 204 ns

Travel time is \(\int s\,dl\), linear in s, so the forward model, the Gaussian noise model and quadratic penalties are all natural in slowness (13.18). Near 1540 m/s, \(d\Delta s/dc = -1/c_0^2\) = −0.42 ns/mm per m/s. A 10 mm path through a +50 m/s inclusion advances the echo by 204 ns, a full period at 5 MHz: absolute phases wrap, which is another reason to measure differences between nearby angles.

Exact formula

Linearization around c₀ holds for small contrast (a few %); straight rays ignore refraction, which bends rays at fat/muscle boundaries by several degrees at oblique incidence.

Never commit the inverse crime: simulate data with a finer, different model than the one you invert.

data: 0.05–0.1 mm rays, smooth edgesmodel: 0.4 mm bilinear

If the data are generated by the same \(\mathbf{A}\) you invert, they lie exactly in its range, model error is zero, and every regularizer looks better than it is. The tomography demo and script integrate a smooth-edged truth with fine steps and invert with coarse bilinear weights on a 1.6 mm grid, so discretization error behaves like correlated noise, as it does in practice.

Exact formulatanh-edged truth and its bilinear interpolation on a 1.6 mm grid, computed directly.

A finer straight-ray model still shares the straight-ray physics; validating against full-wave simulations and phantoms is the real test.

13.11Cheat sheet

Concept	Formula	Ultrasound use
Forward model	\(\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{n}\)	RF = pulse ∗ reflectivity; delays = ray integrals of slowness
Naive inverse error	\(\sum_i (\mathbf{u}_i^{T}\mathbf{n}/\sigma_i)\,\mathbf{v}_i\)	why out-of-band and deep components explode
Tikhonov (standard form)	\((\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{I})^{-1}\mathbf{A}^{T}\mathbf{y}\), \(f_i = \sigma_i^2/(\sigma_i^2 + \lambda^2)\)	deconvolution, diagonal loading (Module 05 writes λ for our \(\lambda^2\))
Tikhonov (general form)	\((\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{L}^{T}\mathbf{L})^{-1}\mathbf{A}^{T}\mathbf{y}\)	smooth or layered sound-speed maps (L = D₁, D₂, gradient)
Circulant / Wiener	\(X = H^*Y/(\lvert H\rvert^2 + \lambda^2)\), \(\lambda^2 = \sigma_n^2/\sigma_x^2\)	RF and PSF deconvolution in the frequency domain
Noise gain cap	\(\max\ \sigma/(\sigma^2 + \lambda^2) = 1/(2\lambda)\)	worst-case amplification of a regularized inverse
Model resolution	\(\mathbf{R} = \mathbf{A}^{\#}\mathbf{A}\), bias \((\mathbf{R} - \mathbf{I})\mathbf{x}\), \(\operatorname{trace}\mathbf{R} = \sum f_i\)	resolution and effective dof of SoS maps
Discrepancy principle	\(\lVert\mathbf{y} - \mathbf{A}\mathbf{x}_\lambda\rVert = \tau\sqrt{m}\,\sigma_n\)	λ or stopping when the noise level is known
GCV	\(m\lVert\mathbf{r}\rVert^2/(m - \sum f_i)^2\)	λ without \(\sigma_n\); white noise only
L-curve	max curvature of \((\log\lVert\mathbf{r}\rVert,\ \log\lVert\mathbf{L}\mathbf{x}\rVert)\)	visual sanity check of λ
Landweber	\(\mathbf{x} \leftarrow \mathbf{x} + \omega\mathbf{A}^{T}(\mathbf{y} - \mathbf{A}\mathbf{x})\), \(\lambda^2 \approx \ln 2/(k\omega)\)	matrix-free tomography, back-projection iterations
CGLS	Krylov least squares, stop at \(\lVert\mathbf{r}_k\rVert \le \tau\sqrt{m}\,\sigma_n\)	large SoS and model-based beamforming problems
Soft threshold	\(S_\tau(v) = \operatorname{sign} v\cdot\max(\lvert v\rvert - \tau, 0)\)	sparse spikes, ULM at high density
ISTA / FISTA	step \(1/\lVert\mathbf{A}\rVert^2\); gaps \(O(1/k)\) / \(O(1/k^2)\)	L1 deconvolution and sparse beamforming
Universal threshold	\(\mu \approx \sigma_n\lVert\mathbf{a}\rVert\sqrt{2\ln N}\)	first guess for the L1 weight
TV by ADMM	\(\tfrac12\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2 + \mu\lVert\mathbf{D}\mathbf{x}\rVert_1\); \(\mathbf{z} \leftarrow S_{\mu/\rho}(\mathbf{D}\mathbf{x} + \mathbf{u})\)	layered tissue, inclusions with sharp edges
TV contrast loss	plateau drops \(2\mu/w\) (1D denoising)	why small inclusions lose contrast under TV
Slowness	\(\Delta s = 1/c - 1/c_0\); 1 ns/mm ≈ 2.37 m/s	units of SoS tomography
Limited angle	k-space wedge \(\pm\theta_{\max}\); layer sensitivity \(2(\sec\theta - 1)\)	axial smearing and layer blindness of echo-mode SoS

13.12Interview and qualifying-exam questions

I1Why is deconvolution of an RF line ill-posed, and what is in its null space? · screen

Convolution with the pulse is diagonal in frequency with singular values \(\lvert H(f)\rvert\). A transducer passes a band around \(f_0\) and nothing at DC or far above, so those frequencies of the reflectivity are (numerically) in the null space: for a 5 MHz, 60 % pulse, \(\lvert H\rvert < 0.01\) outside about 1.1–8.9 MHz. Dividing by H amplifies noise there by \(1/\lvert H\rvert\), over 100 dB at 12 MHz. Tikhonov returns the minimum-norm, band-limited estimate; anything outside the band must come from a prior.

Follow-up: "Then how do sparse deconvolution methods claim higher resolution?" They assume few reflectors; the prior extrapolates the spectrum, and it fails when the assumption does (for example two same-sign reflectors λ/4 apart, which nearly cancel in band).

I2Write the Tikhonov filter factor. A colleague says "we used λ = 0.01": what do you ask? · screen

\(f_i = \sigma_i^2/(\sigma_i^2 + \lambda^2)\) for the penalty \(\lambda^2\lVert\mathbf{x}\rVert^2\). I ask whether the penalty was \(\lambda\lVert\mathbf{x}\rVert^2\) or \(\lambda^2\lVert\mathbf{x}\rVert^2\) (a factor 100 apart at 0.01), whether A was normalized (λ is meaningful only relative to \(\sigma_1\)), what L was, and how λ was chosen. With λ relative to \(\sigma_1\) in the squared convention, components below 1 % of \(\sigma_1\) are suppressed and no noise component is amplified more than \(1/(2\lambda)\) = 50 times.

Follow-up: "Which components are amplified most?" Those with σ = λ, not the smallest ones.

I3How would you choose λ on in-vivo data, and how can each method fail? · onsite

Discrepancy principle if I can measure the noise level (repeat acquisitions, or a noise-only region): \(\lVert\mathbf{r}\rVert = \tau\sqrt{m}\,\sigma_n\); it is robust but over-smooths a little and is wrong if \(\sigma_n\) is wrong. GCV needs no \(\sigma_n\) but assumes white noise: with correlated noise (delay estimates correlated over the kernel, speckle-induced errors) it picks a tiny λ; in the RF demo with AR(1) noise its error is 31 times the oracle's. The L-curve corner is a heuristic that can over-smooth for smooth solutions and has no corner at very low noise. I would cross-check two rules, test stability to a new noise draw, and in simulation report the oracle too.

Follow-up: "Can you fix GCV for correlated noise?" Whiten with an estimate of the noise covariance first, then apply GCV to the whitened problem.

I4Your CG reconstruction gets worse after 40 iterations. Is that a bug? · onsite

No, it is semiconvergence. CGLS first fits the large singular components (signal) and then the small ones (noise), so the error has a minimum. The iteration count plays the role of λ: for Landweber with \(\omega = 1/\sigma_1^2\), \(k \approx \ln 2\,(\sigma_1/\lambda)^2\); CGLS gets there in roughly the square root of that (35 versus 1823 in the RF example). Stop with the discrepancy principle on \(\lVert\mathbf{r}_k\rVert\), or add an explicit penalty and iterate to convergence.

Follow-up: "Which is better, early stopping or explicit Tikhonov?" Equivalent in error in this example; explicit penalties make the result independent of the solver and let you compute R, early stopping is cheaper for very large problems.

I5Compare L2 and L1 regularization for point reflectors. Why soft thresholding, why FISTA, and what about the amplitudes? · onsite

L2 spreads energy over all components and returns band-limited blobs; L1 prefers few nonzeros and can resolve reflectors closer than the pulse length if their echoes stay distinguishable in band. Soft thresholding is the proximal map of \(\mu\lvert x\rvert\), so ISTA alternates a gradient step on the data term with thresholding; FISTA adds Nesterov momentum and improves the worst-case gap from \(O(1/k)\) to \(O(1/k^2)\). L1 shrinks surviving amplitudes by about \(\mu/\lVert\mathbf{a}\rVert^2\), so I would debias by least squares on the recovered support before quantifying intensities, and start μ at \(\sigma_n\lVert\mathbf{a}\rVert\sqrt{2\ln N}\).

Follow-up: "What breaks L1 recovery?" Highly coherent dictionaries (fine grids, band-limited pulses), off-grid sources, and pairs that cancel in band.

I6When would you use TV for a sound-speed map, and what artifacts would you warn a reviewer about? · onsite

When the tissue is approximately piecewise constant (abdominal wall layers, a lesion in a homogeneous organ) and edges matter. TV preserves edges but lowers small plateaus by about \(2\mu/w\) (1D denoising), turns gradients into staircases, and with the anisotropic form favours axis-aligned edges. In the layered example it roughly halves the RMSE compared with the best quadratic penalty; on a ramp it ties \(\mathbf{D}_2\). I would show a ramp phantom and a small-inclusion phantom to expose both artifacts and quote values from plateau interiors.

Follow-up: "How do you solve it?" ADMM with \(\mathbf{z} = \mathbf{D}\mathbf{x}\): one factorization of \(\mathbf{A}^{T}\mathbf{A} + \rho\mathbf{D}^{T}\mathbf{D}\), then soft thresholding; ρ affects speed, not the minimizer.

I7Give the Bayesian interpretation of Tikhonov and of L1. What does it tell you about λ? · research

With Gaussian noise of variance \(\sigma_n^2\), MAP minimizes \(\lVert\mathbf{A}\mathbf{x} - \mathbf{y}\rVert^2/(2\sigma_n^2) - \log p(\mathbf{x})\). A zero-mean Gaussian prior with variance \(\sigma_x^2\) gives Tikhonov with \(\lambda^2 = \sigma_n^2/\sigma_x^2\), the noise-to-signal ratio, which is the Wiener filter in the circulant case; a Laplace prior of scale β gives L1 with \(\mu = \sigma_n^2/\beta\). So λ should rise when the noise rises and fall when the expected contrast rises. MAP gives a point estimate; the sparsity comes from the mode, not from the prior, whose samples are not sparse.

Follow-up: "How would you get error bars?" From the posterior: for the Gaussian case its covariance is \(\sigma_n^2(\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{L}^{T}\mathbf{L})^{-1}\); for non-Gaussian priors, sampling (MCMC) or approximations.

I8Why do echo-mode sound-speed images smear inclusions in depth, and what would you change in the acquisition? · research

Straight rays from each pixel to a linear array lie within \(\pm\theta_{\max}\) of vertical. By the Fourier-slice theorem, their integrals determine only spatial frequencies in a wedge of half-angle \(\theta_{\max}\) around \(k_x\); the cone around \(k_z\) is missing, so the resolution kernel is elongated axially, roughly by \(1/\sin\theta_{\max}\) (9.4 versus 1.2 mm at 10°, 3.3 versus 1.3 mm at 30° in the toy). Differences between angle pairs also make layered structures nearly invisible (sensitivity \(2(\sec\theta - 1)\)). I would widen the angle span as far as the element directivity allows, add angle pairs with large mid-angle differences, combine with absolute-delay information such as depth-resolved autofocus for layers, and use a prior that matches the anatomy.

Follow-up: "What limits \(\theta_{\max}\)?" Element directivity, grating lobes at large steering, aperture coverage of the oblique rays at depth (fewer valid rows), and speckle decorrelation between widely different angles.

I9A reviewer asks for the resolution of your sound-speed map. What do you report? · research

For a linear (Tikhonov-type) reconstruction, rows of the model-resolution matrix \(\mathbf{R} = (\mathbf{A}^{T}\mathbf{A} + \lambda^2\mathbf{L}^{T}\mathbf{L})^{-1}\mathbf{A}^{T}\mathbf{A}\) at representative pixels: their axial and lateral FWHM, how they vary with depth, and \(\operatorname{trace}\mathbf{R}\) as the effective number of resolved parameters. I would add the noise standard deviation map from the posterior covariance, and for nonlinear regularizers (TV, L1) point-target and edge-response simulations at the same λ. Inclusion-mean accuracy is not resolution.

Follow-up: "Can R be computed for a 3D problem?" Not densely; compute selected rows by solving with unit vectors (one solve each), or estimate the diagonal with randomized probing.

I10How do you validate a regularized reconstruction method without fooling yourself? · research

Avoid the inverse crime: generate data on a finer grid with a different (preferably full-wave) model and realistic, correlated noise. Report the oracle parameter and the parameter from a practical rule. Test priors against phantoms that violate them (gradients for TV, extended targets for L1). Show stability to new noise draws and to model error (wrong c₀, refraction). Then phantoms with independently known sound speed, then in vivo with a reference such as a transmission measurement or a known anatomical value.

Follow-up: "What is the signature of the inverse crime?" Near-perfect reconstructions at unrealistically small λ, and errors that fall without bound as noise goes to zero.

I11Why invert for slowness rather than sound speed, and what does the linearization ignore? · screen

Travel time is the path integral of slowness, so the forward model is linear in \(s = 1/c\) and the Gaussian delay-noise model maps to a quadratic data term. Near 1540 m/s, 1 ns/mm of slowness equals 2.37 m/s (\(\Delta c \approx -c_0^2\,\Delta s\)). The linear model assumes straight rays (no refraction), small contrast and that the misfocused echo stays at the same pixel; larger contrasts call for iterating with ray bending or a wave-based model.

Follow-up: "How big is the delay error from a 50 m/s inclusion?" 20.4 ns/mm, so 204 ns over 10 mm, a full period at 5 MHz.

13.13MATLAB lab

Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. Paste into the editor and press Run. Printed numbers below are from Octave; MATLAB's random generator gives slightly different noise draws.

Download the scripts: m13_tikhonov_deconv_lcurve.m · m13_cgls_semiconvergence.m · m13_ista_sparse_spikes.m · m13_tv_vs_tikhonov.m · m13_sos_tomography_limited_angle.m. Base MATLAB R2019b+ or GNU Octave 8, no toolboxes; all scripts in the repository.

m13_tikhonov_deconv_lcurve.mTikhonov RF deconvolution; oracle vs discrepancy vs GCV vs L-curve for white and AR(1) noise · change snr_dB, bw, rho_ar · runtime ≈ 1 s

%% m13_tikhonov_deconv_lcurve.m
% Tikhonov deconvolution of an RF line and four ways to choose lambda.
% Purpose : RF = pulse * reflectivity + noise (circulant model, so the SVD is
%           the DFT and sigma_i = |H(f_i)|). Tikhonov in Hansen's convention
%              x_lam = argmin ||A x - b||^2 + lam^2 ||x||^2,
%              filter factors f_i = sigma_i^2 / (sigma_i^2 + lam^2),
%           so lam is directly comparable to the singular values (|H| <= 1).
%           Compare the oracle lambda (needs the truth) with the discrepancy
%           principle, GCV and the L-curve corner, for white and for
%           correlated (AR(1), rho = 0.9) noise.
% Change  : snr_dB, the pulse bandwidth bw, the AR coefficient rho_ar.
% Expect  : white noise: GCV and L-curve land near the oracle, the
%           discrepancy principle (tau = 1) over-smooths a little.
%           AR(1) noise: GCV and the L-curve pick a tiny lambda and the error
%           explodes (they assume white noise); the discrepancy principle,
%           given the true noise norm, still works.
% Runtime : about 1 s.
clear; close all; rng(1);

%% Pulse and forward operator (circulant, applied by FFT)
N = 512; fs = 40e6; f0 = 5e6; bw = 0.6;        % samples, Hz, Hz, -6 dB fractional bandwidth
sf = bw*f0/2/sqrt(2*log(2));                   % spectral std of the Gaussian amplitude spectrum (Hz)
st = 1/(2*pi*sf);                              % envelope std (s)
t  = ((0:N-1)' - N/2)/fs;
h  = exp(-t.^2/(2*st^2)).*cos(2*pi*f0*t);
h  = circshift(h, -N/2);                       % pulse centred at sample 0 (zero-phase)
H  = fft(h); H = H/max(abs(H));                % normalise so sigma_max = 1
f  = (0:N-1)'*fs/N; f(f >= fs/2) = f(f >= fs/2) - fs;   % frequency axis (Hz)
Aop  = @(v) real(ifft(H.*fft(v)));

%% Reflectivity: a few interfaces plus weak diffuse scattering
x = 0.08*randn(N,1);
pos = [70 150 160 260 330 345 420]; amp = [1 -0.7 0.6 0.9 -0.5 0.5 0.8];
x(pos) = x(pos) + amp(:);
band = abs(H) >= 0.01;                         % band within 40 dB of the peak
X = fft(x); X(~band) = 0; xB = real(ifft(X));  % recoverable (band-limited) truth
fprintf('Recoverable band (|H| >= 0.01): %.2f to %.2f MHz, holds %.0f %% of the reflectivity energy\n', ...
    min(f(band & f > 0))/1e6, max(f(band & f > 0))/1e6, 100*norm(xB)^2/norm(x)^2);
b0 = Aop(x);

lams = logspace(-4, 0.5, 200);
snr_dB = 25; rho_ar = 0.9;
names = {'white', 'AR(1) rho=0.9'};
for nt = 1:2
    w = randn(N,1);
    if nt == 2, n = filter(1, [1 -rho_ar], w); else, n = w; end
    n = n*norm(b0)/norm(n)*10^(-snr_dB/20);    % set SNR
    sig = norm(n)/sqrt(N);                     % noise std (known for the discrepancy principle)
    b = b0 + n; B = fft(b);
    E = zeros(size(lams)); R = E; Xn = E; G = E;
    for k = 1:numel(lams)
        F  = abs(H).^2./(abs(H).^2 + lams(k)^2);           % filter factors
        xl = real(ifft(conj(H)./(abs(H).^2 + lams(k)^2).*B));
        r  = b - Aop(xl);
        E(k) = norm(xl - xB)/norm(xB);  R(k) = norm(r);  Xn(k) = norm(xl);
        G(k) = N*sum(r.^2)/(N - sum(F))^2;                 % GCV function
    end
    [~, io] = min(E); [~, ig] = min(G);
    [~, id] = min(abs(R - sqrt(N)*sig));                   % discrepancy, tau = 1
    lr = log(R); lx = log(Xn); tl = log(lams);             % L-curve curvature (log-log)
    d1r = gradient(lr, tl); d1x = gradient(lx, tl); d2r = gradient(d1r, tl); d2x = gradient(d1x, tl);
    kap = (d1r.*d2x - d2r.*d1x)./(d1r.^2 + d1x.^2).^1.5;
    [~, il] = max(kap);
    fprintf('\n%s noise, SNR %d dB (noise std %.4f)\n', names{nt}, snr_dB, sig);
    fprintf('  %-22s lambda = %.4f   rel. error = %.3f\n', 'oracle (needs truth)', lams(io), E(io));
    fprintf('  %-22s lambda = %.4f   rel. error = %.3f\n', 'discrepancy (tau=1)', lams(id), E(id));
    fprintf('  %-22s lambda = %.4f   rel. error = %.3f\n', 'GCV', lams(ig), E(ig));
    fprintf('  %-22s lambda = %.4f   rel. error = %.3f\n', 'L-curve corner', lams(il), E(il));
    if nt == 1
        Rw = R; Xw = Xn; Ew = E; Gw = G; iw = [io id ig il]; bw_ = b;
    end
end

%% Figures (white-noise case)
figure('Position', [50 50 1000 700]);
subplot(2,2,1);
loglog(Rw, Xw, 'k-'); hold on;
mk = {'o','s','d','^'}; cl = {'b','m','r',[0 0.6 0.4]};
for q = 1:4, loglog(Rw(iw(q)), Xw(iw(q)), mk{q}, 'Color', cl{q}, 'MarkerFaceColor', cl{q}); end
xlabel('residual norm ||b - A x_\lambda||'); ylabel('solution norm ||x_\lambda||');
title('L-curve (white noise)'); legend('L-curve','oracle','discrepancy','GCV','corner','Location','northeast');
subplot(2,2,2);
semilogx(lams, Ew, 'b-'); hold on;
semilogx(lams, Gw/max(Gw)*max(Ew), 'r--');
xlabel('\lambda (units of \sigma_{max})'); ylabel('relative error vs band-limited truth');
title('error (blue) and GCV (red, scaled)'); ylim([0 1.5]);
xl = real(ifft(conj(H)./(abs(H).^2 + lams(iw(1))^2).*fft(bw_)));
tt = (0:N-1)/fs*1e6;
subplot(2,2,3);
plot(tt, xB, 'Color', [0.6 0.6 0.6]); hold on; plot(tt, xl, 'b');
xlabel('time (\mus)'); ylabel('reflectivity'); title('band-limited truth (grey) and oracle Tikhonov (blue)');
subplot(2,2,4);
fp = f(1:N/2)/1e6; F = abs(H(1:N/2)).^2./(abs(H(1:N/2)).^2 + lams(iw(1))^2);
plot(fp, abs(H(1:N/2)), 'k--'); hold on; plot(fp, F, 'b');
xlabel('frequency (MHz)'); ylabel('gain'); title('|H(f)| (dashed) and filter factors at oracle \lambda');

What to look for: recoverable band 1.17–8.83 MHz. White noise at 25 dB: oracle λ = 0.0203 (error 0.334), GCV 0.0203 (0.334), corner 0.0249 (0.341), discrepancy 0.0707 (0.445). AR(1): GCV 0.0003 with error 11.8; the corner 0.028 happens to survive this draw.

m13_cgls_semiconvergence.mCGLS and Landweber error vs iteration, discrepancy stop, Landweber k vs ln 2/λ² · change snr_dB, tau · runtime ≈ 1 s

%% m13_cgls_semiconvergence.m
% Iteration count as the regularization parameter: CGLS and Landweber.
% Purpose : Same RF deconvolution problem as m13_tikhonov_deconv_lcurve.m.
%           Run CGLS (conjugate gradients on the normal equations, never
%           forming A'A) and Landweber x <- x + w A'(b - A x), w = 1/sigma_max^2.
%           The error against the band-limited truth first falls, then rises
%           (semiconvergence). Landweber after k steps has filter factors
%           1 - (1 - w sigma^2)^k, half-power at sigma^2 ~ ln2/k, so it acts
%           like Tikhonov with lambda^2 ~ ln2/k.
% Change  : snr_dB; kmax; the discrepancy factor tau.
% Expect  : at 25 dB the best CGLS iterate is reached in tens of iterations,
%           Landweber needs about ln2/lambda_opt^2 (over a thousand), and
%           both reach about the oracle-Tikhonov error.
% Runtime : about 1 s.
clear; close all; rng(1);

N = 512; fs = 40e6; f0 = 5e6; bw = 0.6;
sf = bw*f0/2/sqrt(2*log(2)); st = 1/(2*pi*sf);
t = ((0:N-1)' - N/2)/fs;
h = circshift(exp(-t.^2/(2*st^2)).*cos(2*pi*f0*t), -N/2);
H = fft(h); H = H/max(abs(H));
Aop  = @(v) real(ifft(H.*fft(v)));
Atop = @(v) real(ifft(conj(H).*fft(v)));        % adjoint (correlation with the pulse)

x = 0.08*randn(N,1);
pos = [70 150 160 260 330 345 420]; amp = [1 -0.7 0.6 0.9 -0.5 0.5 0.8];
x(pos) = x(pos) + amp(:);
band = abs(H) >= 0.01; X = fft(x); X(~band) = 0; xB = real(ifft(X));
b0 = Aop(x);
w = randn(N,1);   % one draw reused so the three SNRs differ only in level

tau = 1.0; kmax = 400;
for snr_dB = [40 25 10]
    n = w*norm(b0)/norm(w)*10^(-snr_dB/20); sig = norm(n)/sqrt(N); b = b0 + n;
    % ---- CGLS
    xk = zeros(N,1); r = b; s = Atop(r); p = s; g = s'*s;
    Ec = zeros(kmax,1); Rc = Ec;
    for k = 1:kmax
        q = Aop(p); a = g/(q'*q);
        xk = xk + a*p; r = r - a*q;
        s = Atop(r); gn = s'*s; p = s + (gn/g)*p; g = gn;
        Ec(k) = norm(xk - xB)/norm(xB); Rc(k) = norm(r);
    end
    [ec, kc] = min(Ec); kdp = find(Rc <= tau*sqrt(N)*sig, 1);
    % ---- Landweber, closed form in the Fourier domain (w = 1)
    ks = unique(round(logspace(0, 5, 300))); El = zeros(size(ks)); B = fft(b);
    Hs = H; Hs(abs(Hs) < 1e-300) = 1;
    for q = 1:numel(ks)
        F = 1 - (1 - abs(H).^2).^ks(q);
        El(q) = norm(real(ifft(F./Hs.*B)) - xB)/norm(xB);
    end
    [el, ql] = min(El);
    % ---- Tikhonov oracle for comparison
    lams = logspace(-4, 0.5, 300); Et = zeros(size(lams));
    for q = 1:numel(lams)
        Et(q) = norm(real(ifft(conj(H)./(abs(H).^2 + lams(q)^2).*B)) - xB)/norm(xB);
    end
    [et, qt] = min(Et);
    fprintf('SNR %2d dB: CGLS best k = %3d (err %.3f), discrepancy stop k = %3d (err %.3f), err at k=%d: %.2f\n', ...
        snr_dB, kc, ec, kdp, Ec(kdp), kmax, Ec(end));
    fprintf('           Landweber best k = %5d (err %.3f); Tikhonov oracle lambda = %.4f (err %.3f); ln2/lambda^2 = %.0f\n', ...
        ks(ql), el, lams(qt), et, log(2)/lams(qt)^2);
    if snr_dB == 25, Ec25 = Ec; El25 = El; et25 = et; kc25 = kc; kdp25 = kdp; end
end

figure('Position', [50 50 900 380]);
subplot(1,2,1);
semilogx(1:kmax, Ec25, 'b-'); hold on; plot(kc25, Ec25(kc25), 'bo'); plot(kdp25, Ec25(kdp25), 'ms');
plot([1 kmax], [et25 et25], 'k:'); ylim([0 1.5]);
xlabel('CGLS iteration k'); ylabel('relative error'); title('CGLS, 25 dB: semiconvergence');
legend('error', 'best', 'discrepancy stop', 'Tikhonov oracle', 'Location', 'northwest');
subplot(1,2,2);
semilogx(ks, El25, 'r-'); hold on; plot([1 1e5], [et25 et25], 'k:'); ylim([0 1.5]);
xlabel('Landweber iteration k'); ylabel('relative error'); title('Landweber (\omega = 1/\sigma_{max}^2), 25 dB');

What to look for: at 40 / 25 / 10 dB, CGLS best k = 71 / 35 / 6, Landweber best k = 8841 / 1823 / 33 against \(\ln 2/\lambda_{\mathrm{opt}}^2\) = 9076 / 1720 / 35, and the same minimum error as oracle Tikhonov. At k = 400 the CGLS error is 1.5 / 8.5 / 47: semiconvergence.

m13_ista_sparse_spikes.mL1 sparse-spike deconvolution with ISTA and FISTA, convergence bounds, debiasing, Tikhonov comparison · change spikes, snr_dB, lam · runtime ≈ 5 s

%% m13_ista_sparse_spikes.m
% Sparse-spike deconvolution: L2 (Tikhonov) vs L1 (ISTA and FISTA).
% Purpose : Recover point reflectors from a band-limited RF line by
%              min_x 0.5 ||A x - b||^2 + lam ||x||_1
%           with ISTA  x <- soft(x - A'(A x - b)/Lip, lam/Lip)
%           and FISTA (Beck & Teboulle 2009, Nesterov momentum), Lip = sigma_max^2 = 1.
%           Check the worst-case bounds  F(x_k) - F* <= Lip ||x0 - x*||^2 / (2k)    (ISTA)
%                                         F(x_k) - F* <= 2 Lip ||x0 - x*||^2 / (k+1)^2 (FISTA).
%           lam starts from the universal threshold sigma_n ||h|| sqrt(2 ln N).
%           Debias: least squares on the recovered support removes the L1 shrinkage.
% Change  : spike list, snr_dB, lam_factor, number of iterations.
% Expect  : FISTA reaches a 1e-6 relative gap in a few hundred iterations,
%           ISTA does not within 1000; L1 recovers the isolated spikes and the
%           opposite-sign pair 6 samples apart; the same-sign pair 4 samples apart
%           (0.1 us = half a period at 5 MHz, lambda/4 in depth) is lost: its two
%           echoes nearly cancel in band (1 - 0.9 = 0.1 at f0), so no prior
%           can recover it and L1 returns wrong-sign spikes at the wrong positions.
% Runtime : about 2 s.
clear; close all; rng(2);

N = 256; fs = 40e6; f0 = 5e6; bw = 0.6;
sf = bw*f0/2/sqrt(2*log(2)); st = 1/(2*pi*sf);
t = ((0:N-1)' - N/2)/fs;
h = circshift(exp(-t.^2/(2*st^2)).*cos(2*pi*f0*t), -N/2);
H = fft(h); H = H/max(abs(H)); h = real(ifft(H));
Aop  = @(v) real(ifft(H.*fft(v)));
Atop = @(v) real(ifft(conj(H).*fft(v)));
soft = @(v, th) sign(v).*max(abs(v) - th, 0);
fprintf('Pulse: -6 dB envelope duration %.0f ns = %.1f samples; ||h|| = %.3f\n', ...
    2*sqrt(2*log(2))*st*1e9, 2*sqrt(2*log(2))*st*fs, norm(h));

x = zeros(N,1);
pos = [30 70 76 120 170 174 215]; amp = [1 -0.8 0.7 0.6 -1 -0.9 0.5];
x(pos + 1) = amp(:);                             % 0-based sample positions
b0 = Aop(x); n = randn(N,1); snr_dB = 25;
n = n*norm(b0)/norm(n)*10^(-snr_dB/20); sig = norm(n)/sqrt(N); b = b0 + n;
lamU = sig*norm(h)*sqrt(2*log(N));
lam = 1.0*lamU;
Fobj = @(v) 0.5*sum((Aop(v) - b).^2) + lam*sum(abs(v));
fprintf('noise std %.4f, universal lambda %.4f, using lambda = %.4f\n', sig, lamU, lam);

%% Reference minimizer (long FISTA run)
xs = zeros(N,1); y = xs; tk = 1;
for k = 1:20000
    xn = soft(y - Atop(Aop(y) - b), lam); tn = (1 + sqrt(1 + 4*tk^2))/2;
    y = xn + ((tk - 1)/tn)*(xn - xs); xs = xn; tk = tn;
end
Fs = Fobj(xs); R2 = sum(xs.^2);

%% ISTA vs FISTA from x0 = 0
K = 1000; gi = zeros(K,1); gf = gi;
xi = zeros(N,1); xf = zeros(N,1); y = xf; tk = 1;
for k = 1:K
    xi = soft(xi - Atop(Aop(xi) - b), lam); gi(k) = Fobj(xi) - Fs;
    xn = soft(y - Atop(Aop(y) - b), lam); tn = (1 + sqrt(1 + 4*tk^2))/2;
    y = xn + ((tk - 1)/tn)*(xn - xf); xf = xn; tk = tn; gf(k) = Fobj(xf) - Fs;
end
kk = (1:K)'; bI = R2./(2*kk); bF = 2*R2./(kk + 1).^2;
fprintf('relative gap (F-F*)/F* at k = 10, 100, 1000:  ISTA %.1e %.1e %.1e   FISTA %.1e %.1e %.1e\n', ...
    gi(10)/Fs, gi(100)/Fs, gi(1000)/Fs, gf(10)/Fs, gf(100)/Fs, gf(1000)/Fs);
fprintf('worst-case bounds hold: ISTA %d, FISTA %d\n', all(gi <= bI + 1e-12), all(gf <= bF + 1e-12));
k6i = find(gi/Fs < 1e-6, 1); k6f = find(gf/Fs < 1e-6, 1);
if isempty(k6i), k6i = NaN; end
fprintf('iterations to a 1e-6 relative gap: ISTA %g, FISTA %d\n', k6i, k6f);

%% Support, debiasing, and the L2 comparison
supp = abs(xs) > 0.05*max(abs(xs));
xd = zeros(N,1); idx = find(supp);
Ad = zeros(N, numel(idx));
for q = 1:numel(idx), e = zeros(N,1); e(idx(q)) = 1; Ad(:,q) = Aop(e); end
xd(idx) = Ad\b;                                  % least squares on the support
fprintf('L1 support (|x| > 5%% of max): %d spikes at samples %s\n', numel(idx), mat2str(idx' - 1));
fprintf('L1 amplitudes on that support: %s\n', mat2str(round(xs(idx)'*100)/100));
fprintf('true spikes at samples %s\n', mat2str(pos));
fprintf('L1 amplitudes at true positions:       %s\n', mat2str(round(xs(pos + 1)'*100)/100));
fprintf('debiased amplitudes at true positions: %s\n', mat2str(round(xd(pos + 1)'*100)/100));
lams = logspace(-4, 0.5, 200); Et = zeros(size(lams)); B = fft(b);
for q = 1:numel(lams)
    Et(q) = norm(real(ifft(conj(H)./(abs(H).^2 + lams(q)^2).*B)) - x)/norm(x);
end
[et, qt] = min(Et); xt = real(ifft(conj(H)./(abs(H).^2 + lams(qt)^2).*B));
fprintf('relative error: Tikhonov (oracle lambda %.4f) %.3f | L1 %.3f | L1 debiased %.3f\n', ...
    lams(qt), et, norm(xs - x)/norm(x), norm(xd - x)/norm(x));
% (Indices printed 0-based to match the browser demo; MATLAB arrays are 1-based,
%  so x(pos+1) is the spike listed at sample pos.)

%% Figures
tt = (0:N-1)/fs*1e6;
figure('Position', [50 50 1000 650]);
subplot(2,2,1); plot(tt, b, 'Color', [0.8 0.5 0]); xlabel('time (\mus)'); title('RF data b');
subplot(2,2,2); stem(tt(pos + 1), amp, 'k', 'filled'); hold on; plot(tt, xt, 'm');
xlabel('time (\mus)'); title('truth (stems) and Tikhonov (oracle \lambda)');
subplot(2,2,3); stem(tt(pos + 1), amp, 'k', 'filled'); hold on; stem(tt(idx), xs(idx), 'b');
xlabel('time (\mus)'); title('truth (black) and L1 / FISTA (blue)');
subplot(2,2,4); loglog(kk, gi, 'r', kk, gf, 'b', kk, bI, 'r:', kk, bF, 'b:');
xlabel('iteration k'); ylabel('F(x_k) - F^*'); title('ISTA (red) vs FISTA (blue), dotted: bounds');

What to look for: universal μ = 0.0023; both worst-case bounds hold; relative gap after 1000 iterations \(6.5\times10^{-3}\) (ISTA) vs \(10^{-11}\) (FISTA reaches \(10^{-6}\) at 274). Spikes at 30, 70, 76, 120, 215 are found with amplitudes shrunk by 0.02–0.04 and restored by debiasing; the same-sign pair at 170/174 is replaced by opposite-sign spikes (+0.40, +0.36) at 166 and 178 and a residual −0.07 at 170. Relative error: Tikhonov 0.82, L1 0.66, debiased 0.63.

m13_tv_vs_tikhonov.mAverage-to-local sound speed: naive noise growth, Tikhonov with I/D₁/D₂, TV by ADMM; TV contrast-loss rule · change profile, noise, Δz · runtime ≈ 6 s

%% m13_tv_vs_tikhonov.m
% Average-to-local sound speed: naive inversion, Tikhonov (L = I, D1, D2) and TV.
% Purpose : A depth-dependent average sound speed cbar(z) (for example from
%           autofocus at each depth, vertical-path approximation) relates to
%           the local slowness s(z) = 1/c(z) by the running mean
%              1/cbar(z) = (1/z) * integral_0^z s(z') dz'      (linear in s).
%           Naive inversion differentiates z/cbar(z), amplifying noise by
%           about sqrt(2) z/dz. Regularize the slowness deviation
%           x = s - 1/c0 (ns/mm) with
%              Tikhonov : min ||A x - b||^2 + lam^2 ||L x||^2,  L = I, D1 or D2
%              TV       : min 0.5||A x - b||^2 + mu ||D1 x||_1  (ADMM, rho = 0.1)
%           Also checks the 1D TV-denoising rule: an isolated plateau of width w
%           (samples) and height 1 comes back with height 1 - 2 mu / w.
% Change  : profile ('layers' or 'ramp'), noise std sc (m/s), dz.
% Expect  : naive std grows linearly with depth (about 110 m/s at 30 mm for
%           1 m/s noise on cbar); TV beats every Tikhonov on the layers and keeps
%           the 8 mm step sharp; on the ramp segment (z > 24 mm) TV staircases,
%           so compare the 'z > 24 mm' column, not only the total RMSE.
% Runtime : about 10 s (the ADMM parameter scan dominates).
clear; close all; rng(1);

dz = 0.4; n = 100; z = dz*(1:n)'; c0 = 1540;      % mm, samples, mm, m/s
A  = tril(ones(n))*dz./repmat(z, 1, n);           % running-mean operator
D1 = diff(eye(n)); D2 = diff(eye(n), 2);
sv = svd(A); fprintf('running-mean operator: sigma_max %.2f, sigma_min %.4f, kappa %.0f\n', sv(1), sv(end), sv(1)/sv(end));
toS = @(c) (1./c - 1/c0)*1e6;                     % m/s -> slowness deviation (ns/mm)
toC = @(x) 1./(1/c0 + x*1e-6);                    % back to m/s
prof = @(kind) 1450*(z <= 8) + 1580*(z > 8 & z <= 20) + 1560*(z > 20 & z <= 24) + ...
    (1560 + strcmp(kind, 'ramp')*2.0*(z - 24)).*(z > 24);

%% Naive inversion: noise growth with depth
c = prof('layers'); cbar = 1./(A*(1./c)); sc = 1.0; M = 2000;
err = zeros(M, n);
for m = 1:M
    sb = 1./(cbar + sc*randn(n, 1));
    sn = [sb(1); (z(2:end).*sb(2:end) - z(1:end-1).*sb(1:end-1))/dz];
    err(m, :) = (1./sn - c)';
end
iz = [25 50 75];
% first-order theory: std(c) = (c/cbar)^2 * sqrt(z_i^2 + z_(i-1)^2)/dz * sc  ~  sqrt(2) z/dz * sc
fprintf('naive inverse, cbar noise %.1f m/s: std of c at z = 10, 20, 30 mm: %s m/s (theory %s)\n', ...
    sc, mat2str(round(std(err(:, iz)))), mat2str(round((c(iz)'./cbar(iz)').^2.*sqrt(z(iz).^2 + z(iz-1).^2)'/dz*sc)));

%% Regularized inversions
edge = @(ce) z(find(ce > 1450 + 0.9*130 & z > 3, 1)) - z(find(ce > 1450 + 0.1*130 & z > 3, 1));  % 10-90 % width of the 8 mm step
hf = figure('Position', [50 50 1000 380]);
kinds = {'layers', 'ramp'}; lamg = logspace(-4, 1, 60); mug = logspace(-3, 1, 25); rho = 0.1;
for kk = 1:2
    c = prof(kinds{kk}); sc = 2.0;
    b = toS(1./(A*(1./c)) + sc*randn(n, 1));
    Ls = {eye(n), D1, D2}; Ln = {'I', 'D1', 'D2'};
    fprintf('\n%s, cbar noise %.1f m/s (oracle parameters, RMSE in m/s):\n', kinds{kk}, sc);
    for q = 1:3
        L = Ls{q}; e = zeros(size(lamg));
        for j = 1:numel(lamg), e(j) = sqrt(mean((toC((A'*A + lamg(j)^2*(L'*L))\(A'*b)) - c).^2)); end
        [eb, jb] = min(e); ce = toC((A'*A + lamg(jb)^2*(L'*L))\(A'*b));
        fprintf('  Tikhonov L = %-2s  lambda %.3g  RMSE %5.1f  (z > 24 mm: %5.1f)  edge 10-90 %% %.1f mm\n', Ln{q}, lamg(jb), eb, sqrt(mean((ce(z > 24) - c(z > 24)).^2)), edge(ce));
        if q == 2, cT = ce; end
    end
    Mi = inv(A'*A + rho*(D1'*D1)); Atb = A'*b; e = zeros(size(mug));
    for j = 1:numel(mug)
        x = Mi*Atb; zz = D1*x; u = zeros(n-1, 1); th = mug(j)/rho;
        for it = 1:1500
            x = Mi*(Atb + rho*D1'*(zz - u)); w = D1*x + u; zz = sign(w).*max(abs(w) - th, 0); u = w - zz;
        end
        e(j) = sqrt(mean((toC(x) - c).^2));
        if j == 1 || e(j) < min(e(1:j-1)), cV = toC(x); end
    end
    [eb, jb] = min(e);
    fprintf('  TV (ADMM)       mu     %.3g  RMSE %5.1f  (z > 24 mm: %5.1f)  edge 10-90 %% %.1f mm\n', mug(jb), eb, sqrt(mean((cV(z > 24) - c(z > 24)).^2)), edge(cV));
    subplot(1, 2, kk);
    plot(z, c, 'k', z, cT, 'm', z, cV, 'b', 'LineWidth', 1.2);
    xlabel('depth z (mm)'); ylabel('local sound speed (m/s)'); title([kinds{kk} ': truth, Tikhonov-D1 (magenta), TV (blue)']);
    ylim([1420 1640]);
end

%% TV denoising contrast-loss rule (A = I)
m = 128; D = diff(eye(m)); y = zeros(m, 1); y(41:60) = 1;   % plateau width w = 20 samples
Mi = inv(eye(m) + 1*(D'*D));
fprintf('\nTV denoising of a 20-sample plateau of height 1:\n');
for mu = [0.5 1 2 4]
    x = y; zz = D*x; u = zeros(m-1, 1);
    for it = 1:3000
        x = Mi*(y + D'*(zz - u)); w = D*x + u; zz = sign(w).*max(abs(w) - mu, 0); u = w - zz;
    end
    fprintf('  mu = %.1f: plateau %.3f (rule 1 - 2mu/w = %.3f), background next to it %.4f\n', mu, x(50), 1 - 2*mu/20, x(10));
end

What to look for: κ = 310; naive std 39 / 77 / 111 m/s at 10 / 20 / 30 mm (theory 40 / 75 / 109). Layers, 2 m/s noise: Tikhonov RMSE 13.1 / 12.5 / 13.2 m/s (I / D₁ / D₂), TV 6.0 with a 1.2 mm edge. Ramp segment (z > 24 mm): D₂ 7.2, TV 7.2, D₁ 11.4. TV denoising plateau 0.950 / 0.900 / 0.800 / 0.600 for μ = 0.5 / 1 / 2 / 4, exactly 1 − 2μ/20.

m13_sos_tomography_limited_angle.mEcho-mode SoS tomography toy: common-midangle data, Tikhonov-gradient and TV, resolution row vs angle span · change thmax_list, noise_ns, kind · runtime ≈ 7 s

%% m13_sos_tomography_limited_angle.m
% Echo-mode (CUTE-type) sound-speed tomography toy with limited angular aperture.
% Purpose : Unknown: slowness deviation ds = 1/c - 1/c0 on a 24 x 20 grid of
%           1.6 mm pixels (x = +-19.2 mm under a 38.4 mm, 128 x 0.3 mm array;
%           z = 3 to 35 mm). For every pixel r (echoes come from speckle
%           everywhere) and steering angle th, the beamforming delay error is
%           the line integral of ds along the straight path from r to the
%           array at angle th:  L_th(r) = int ds dl.
%           Data = common-midangle differences between transmit/receive angle
%           pairs:  d = L_i + L_j - 2 L_m, with m = (i+j)/2, K = 7 angles in
%           [-thmax, thmax]. The absolute delay is never measured.
%           Inversion: Tikhonov with a gradient penalty (Hansen convention)
%              min ||A ds - d||^2 + lam^2 ||G ds||^2
%           and anisotropic TV by ADMM. Data are generated by fine ray
%           sampling of a smooth-edged truth (0.05 mm steps) and inverted with
%           bilinear 0.4 mm steps on the coarse grid, to avoid the inverse crime.
%           The 0-3 mm band above the grid is lumped into the top row (clamped weights).
%           The model-resolution matrix R = (A'A + lam^2 G'G)^-1 A'A shows
%           what the inversion can see: row k = how pixel k mixes its neighbours.
% Change  : thmax_list, noise_ns, the truth (kind = 'vpair', 'disk', 'hpair').
% Expect  : two inclusions 12 mm apart in depth merge into one vertical streak
%           at 10 deg and separate at 30 deg; the resolution row is several
%           times longer in z than in x; a horizontal layer is only weakly
%           sensed (its data scale as 2(sec th - 1) ~ th^2).
% Runtime : about 10 s in Octave.
clear; close all; rng(1);

c0 = 1540;                                   % reference sound speed (m/s)
Nx = 24; Nz = 20; dx = 1.6;                  % grid (pixels, mm)
xc = ((0:Nx-1) - (Nx-1)/2)*dx;               % lateral centres (mm)
zc = 3 + ((0:Nz-1) + 0.5)*dx;                % depth centres (mm)
Xh = 19.2;                                   % aperture half-width (mm)
np = Nx*Nz; K = 7; noise_ns = 5;             % unknowns, angles, delay noise std (ns)
kind = 'vpair';
sg = @(r, r0) 1./(1 + exp((r - r0)/0.4));    % smooth edge, 0.4 mm
switch kind
    case 'vpair', truth = @(x, z) 1540 + 50*(sg(hypot(x, z - 13), 3) + sg(hypot(x, z - 25), 3));
    case 'hpair', truth = @(x, z) 1540 + 50*(sg(hypot(x - 6, z - 19), 3) + sg(hypot(x + 6, z - 19), 3));
    otherwise,    truth = @(x, z) 1540 + 50*sg(hypot(x, z - 19), 5);
end
dsf = @(c) (1./c - 1/c0)*1e6;                % slowness deviation in ns/mm
[XX, ZZ] = meshgrid(xc, zc); ctrue = truth(XX, ZZ); ctrue = reshape(ctrue', [], 1);  % row-major: index = iz*Nx + ix

% gradient operator G (horizontal and vertical neighbour differences)
I = speye(np); G = [];
for iz = 0:Nz-1, for ix = 0:Nx-2, k = iz*Nx + ix + 1; G = [G; I(k + 1, :) - I(k, :)]; end, end %#ok<AGROW>
for iz = 0:Nz-2, for ix = 0:Nx-1, k = iz*Nx + ix + 1; G = [G; I(k + Nx, :) - I(k, :)]; end, end %#ok<AGROW>
GG = full(G'*G);                              % G stays sparse for fast products

fprintf('Layer sensitivity 2(sec th - 1) for the symmetric pair at +-th: 10 deg %.3f, 20 deg %.3f, 30 deg %.3f\n', ...
    2*(sec(pi/18) - 1), 2*(sec(pi/9) - 1), 2*(sec(pi/6) - 1));
kz0 = find(abs(zc - 19) == min(abs(zc - 19)), 1); kx0 = find(abs(xc - 0.8) == min(abs(xc - 0.8)), 1);
kpix = (kz0 - 1)*Nx + kx0;                   % pixel whose resolution row we inspect (x = 0.8, z = 19.0 mm)
fwhm = @(p, dd) dd*0.01*sum(interp1(1:numel(p), p, 1:0.01:numel(p)) >= max(p)/2);
thmax_list = [10 20 30];
figure('Position', [40 40 1100 800]); colormap(jet);
for q = 1:numel(thmax_list)
    tic;
    th = linspace(-thmax_list(q), thmax_list(q), K)*pi/180;
    Lr = zeros(np, np, K); Lt = zeros(np, K); ok = false(np, K);
    for a = 1:K
        S = []; W = []; Rw = [];
        for iz = 1:Nz
            rows = (iz - 1)*Nx + (1:Nx)'; z = zc(iz); x = xc(:);              % one depth row at a time
            ok(rows, a) = abs(x - z*tan(th(a))) <= Xh;
            Lp = z/cos(th(a));
            ns = ceil(Lp/0.4); tt = ((1:ns) - 0.5)*Lp/ns; w = Lp/ns;          % coarse model, 0.4 mm steps
            gx = min(max((x - tt*sin(th(a)) - xc(1))/dx, 0), Nx - 1 - 1e-9);  % Nx x ns grid coordinates
            gz = min(max((z - tt*cos(th(a)) - zc(1))/dx, 0), Nz - 1 - 1e-9);
            gz = repmat(gz, Nx, 1);
            i0 = floor(gx); j0 = floor(gz); fx = gx - i0; fz = gz - j0;
            i1 = min(i0 + 1, Nx - 1); j1 = min(j0 + 1, Nz - 1);
            S  = [S;  j0(:)*Nx + i0(:) + 1; j0(:)*Nx + i1(:) + 1; j1(:)*Nx + i0(:) + 1; j1(:)*Nx + i1(:) + 1]; %#ok<AGROW>
            W  = [W;  w*[(1 - fx(:)).*(1 - fz(:)); fx(:).*(1 - fz(:)); (1 - fx(:)).*fz(:); fx(:).*fz(:)]];    %#ok<AGROW>
            Rw = [Rw; repmat(repmat(rows, ns, 1), 4, 1)];                                                      %#ok<AGROW>
            nf = ceil(Lp/0.05); tf = ((1:nf) - 0.5)*Lp/nf;                     % fine truth, 0.05 mm steps
            Lt(rows, a) = sum(dsf(truth(x - tf*sin(th(a)), z - tf*cos(th(a)))), 2)*Lp/nf;
        end
        Lr(:, :, a) = accumarray([Rw S], W, [np np]);
    end
    A = []; d0 = [];
    for i = 1:K
        for j = i + 2:2:K
            m = (i + j)/2; use = ok(:, i) & ok(:, j) & ok(:, m);
            A  = [A;  Lr(use, :, i) + Lr(use, :, j) - 2*Lr(use, :, m)];   %#ok<AGROW>
            d0 = [d0; Lt(use, i) + Lt(use, j) - 2*Lt(use, m)];            %#ok<AGROW>
        end
    end
    d = d0 + noise_ns*randn(size(d0));
    AtA = A'*A; Atd = A'*d; s1 = sqrt(max(eig(AtA)));
    lr = logspace(-3, 0, 31); e = zeros(size(lr));
    for j = 1:numel(lr)
        x = (AtA + (lr(j)*s1)^2*GG)\Atd; e(j) = sqrt(mean((1./(1/c0 + x*1e-6) - ctrue).^2));
    end
    [eb, jb] = min(e); M = AtA + (lr(jb)*s1)^2*GG;
    x = M\Atd; cT = 1./(1/c0 + x*1e-6);
    Rm = M\AtA; row = reshape(Rm(kpix, :), Nx, Nz)';            % resolution row as an Nz x Nx image
    fprintf('\nthmax = %2d deg: %d data, rms data %.0f ns, sigma_max %.0f\n', thmax_list(q), size(A, 1), sqrt(mean(d0.^2)), s1);
    fprintf('  Tikhonov-gradient, oracle lambda = %.4f sigma_max: RMSE %.1f m/s, max c %.0f m/s (true 1590), trace(R) = %.0f of %d\n', ...
        lr(jb), eb, max(cT), trace(Rm), np);
    fprintf('  resolution row at (x, z) = (%.1f, %.1f) mm: FWHM axial %.1f mm, lateral %.1f mm\n', ...
        xc(kx0), zc(kz0), fwhm(row(:, kx0), dx), fwhm(row(kz0, :), dx));
    % ---- TV by ADMM: min 0.5||A x - d||^2 + mu ||G x||_1
    rho = 1e-2*s1^2; Rc = chol(AtA + rho*GG); mus = logspace(0, 2.5, 8); Gt = G'; ev = zeros(size(mus));
    for j = 1:numel(mus)
        zz = zeros(size(G, 1), 1); u = zz;
        for it = 1:200
            x = Rc\(Rc'\(Atd + rho*Gt*(zz - u))); w = G*x + u;
            zz = sign(w).*max(abs(w) - mus(j)/rho, 0); u = w - zz;
        end
        ev(j) = sqrt(mean((1./(1/c0 + x*1e-6) - ctrue).^2));
        if j == 1 || ev(j) < min(ev(1:j-1)), cV = 1./(1/c0 + x*1e-6); end
    end
    [evb, jv] = min(ev);
    fprintf('  TV (ADMM, 200 it), oracle mu = %.3g: RMSE %.1f m/s, max c %.0f m/s   [%.1f s]\n', mus(jv), evb, max(cV), toc);
    subplot(3, 4, 4*(q - 1) + 1); imagesc(xc, zc, reshape(ctrue, Nx, Nz)', [1520 1600]); axis image; title('truth c (m/s)'); ylabel(sprintf('%d deg: z (mm)', thmax_list(q)));
    subplot(3, 4, 4*(q - 1) + 2); imagesc(xc, zc, reshape(cT, Nx, Nz)', [1520 1600]); axis image; title('Tikhonov-gradient');
    subplot(3, 4, 4*(q - 1) + 3); imagesc(xc, zc, reshape(cV, Nx, Nz)', [1520 1600]); axis image; title('TV');
    subplot(3, 4, 4*(q - 1) + 4); imagesc(xc, zc, row); axis image; title('resolution row R(k,:)'); colorbar;
end
xlabel('x (mm)');

What to look for: layer sensitivity 0.031 / 0.128 / 0.309. For \(\theta_{\max}\) = 10 / 20 / 30°: resolution-row FWHM axial 9.4 / 4.6 / 3.3 mm vs lateral 1.2–1.3 mm, trace R = 66 / 132 / 185, RMSE Tikhonov 7.7 / 6.7 / 4.8 m/s and TV 7.6 / 5.6 / 2.4 m/s. The figure shows the stacked pair merging into a streak at 10° and separating at 30°. The 0–3 mm band above the grid is lumped into the top row, which is why the top-row pixels light up in the resolution rows.

13.14Self-check

Q1σ = (1, 0.1, 0.01) and λ = 0.03 (penalty \(\lambda^2\lVert x\rVert^2\)). Give the three filter factors. What λ would a code that uses the penalty \(\alpha\lVert x\rVert^2\) need for the same result?

\(f = 1/(1 + 0.0009) = 0.999\); \(0.01/(0.01 + 0.0009) = 0.917\); \(0.0001/(0.0001 + 0.0009) = 0.100\). The other convention needs \(\alpha = \lambda^2 = 9\times10^{-4}\).

Q2With λ = 0.03, what is the largest noise amplification of any component, and at which σ? Compare with the naive gain at σ = 0.01.

\(1/(2\lambda) = 16.7\) at σ = 0.03. At σ = 0.01: Tikhonov \(0.01/(0.0001 + 0.0009) = 10\), naive \(1/0.01 = 100\).

Q3Landweber with \(\omega = 1/\sigma_1^2\) should mimic Tikhonov with \(\lambda = 0.01\,\sigma_1\). About how many iterations? How many would you expect for CGLS?

\(k \approx \ln 2/10^{-4} \approx 6900\). CGLS needs roughly the square root, of order 80 to 100, but check with the discrepancy principle rather than trusting the estimate.

Q4A 10 mm path crosses tissue at 1590 m/s instead of 1540 m/s. What is Δs in ns/mm, the delay change, and the phase change at 5 MHz?

\(\Delta s = (1/1590 - 1/1540)\times10^{6} = -20.4\) ns/mm; delay −204 ns; phase 204 ns × 5 MHz × 360° = 367°, more than a full cycle. Absolute phases wrap; differences between nearby angles do not.

Q5TV denoising (\(\tfrac12\lVert x - y\rVert^2 + \mu\sum\lvert\Delta x\rvert\)) of an isolated plateau 10 samples wide and 1 high with μ = 1.5. What height comes back?

\(1 - 2\mu/w = 1 - 3/10 = 0.7\).

Q6Delay noise \(\sigma_n\) = 2 ns, and the slowness deviation is believed to be zero-mean Gaussian with \(\sigma_x\) = 10 ns/mm per pixel. What λ does the MAP reading give for Tikhonov with L = I?

\(\lambda = \sigma_n/\sigma_x = 2/10 = 0.2\), with units of mm (ns divided by ns/mm), the units of the entries of A. It is meaningful only with A in those units.

Q7Average-to-local inversion with Δz = 0.5 mm and \(\sigma_{\bar c}\) = 0.5 m/s. What noise does the naive inverse have at 40 mm?

\(\sqrt{2}\times(40/0.5)\times 0.5 \approx 57\) m/s (times \((c/\bar c)^2\), close to 1 at that depth).

Q8For a symmetric angle pair at ±15°, what fraction of a layer's vertical delay appears in the common-midangle datum?

\(2(\sec 15° - 1) = 2(1.0353 - 1) = 0.071\), about 7 %.

13.15Figures from the literature

Published figures with real or published-benchmark data for the ideas of this module: common-midangle data, layer sensitivity and prior dependence in echo-mode sound-speed imaging, angle span versus resolution, and FISTA's convergence. Links only. Figure numbers for the two Bern papers follow their arXiv versions (cited), which may differ from the journal layout.

Why common-midangle pairs keep speckle correlated

P. Stähli, M. Kuriakose, M. Frenz, M. Jaeger, "Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics 108, 106168, 2020 (open access, CC BY-NC-ND 4.0); arXiv:1902.10639v2. Fig. 1 (b–d: isochrones for different transmit/receive angle pairs). doi:10.1016/j.ultras.2020.106168 · arXiv

Only pairs with the same mid-angle keep the isochrone spacing, so their beamformed speckle stays correlated and the phase shift between them is trackable. That is the datum \(d_{ij}\) of (13.19), and the reason the toy differences each pair against its own mid-angle reference.

Layered phantoms: the forward model and the assumed c₀ both matter

Stähli et al., Ultrasonics 108, 106168, 2020 (as above). Fig. 4 (two-layer fat-over-liver phantom, reconstructed with each forward model and four a priori sound speeds) and Fig. 5 (four-layer abdominal-wall phantom). doi:10.1016/j.ultras.2020.106168

Rows compare forward models, columns the beamforming sound speed. Watch how much the layer values move with the model: differential data sense layers only through the second-order factor \(2(\sec\theta - 1)\) of (13.20), so small model errors are not small relative to the signal.

In vivo abdominal wall and liver sound speed

Stähli et al., Ultrasonics 108, 106168, 2020 (as above). Fig. 7 (B-mode and sound-speed images of a volunteer's abdominal wall and liver, old and new forward model). doi:10.1016/j.ultras.2020.106168

Real tissue: fat, muscle and liver compartments. Compare the axial sharpness of the layer boundaries with the lateral sharpness of features, and the old with the new forward model: the k-space wedge of Section 13.8 predicts axial blur, and the layer values depend on how well the model handles second-order layer sensitivity.

Point-spread functions of a pulse-echo sound-speed reconstruction

M. Jaeger, P. Stähli, N. Korta Martiartu, P. Salemi Yolgunlu, T. Frappart, C. Fraschini, M. Frenz, "Pulse-echo speed-of-sound imaging using convex probes," Phys. Med. Biol. 67(21), 215016, 2022; arXiv:2202.10633. Fig. 7 (tilted-probe phantom with theoretical point-spread functions). doi:10.1088/1361-6560/ac96c6 · arXiv

The published PSFs are rows of the model-resolution matrix \(\mathbf{R}\) of (13.4), the same object as the resolution row in Demo 13.7. Check the PSF shapes against the k-space wedge: the missing cone follows the ray directions, so the elongation should follow the beam axis of the tilted probe rather than the image axes. The paper regularizes with first-order Tikhonov, much stronger laterally than axially.

Sound speed against the maximum Tx–Rx angle spread

Jaeger et al., Phys. Med. Biol. 67, 215016, 2022 (as above). Fig. 8 (angle coverage and sound-speed results for different maximum angle spreads). doi:10.1088/1361-6560/ac96c6

The experimental counterpart of the \(\theta_{\max}\) table in Section 13.8. The toy predicts that a wider angle spread fills more of k-space and reduces axial smearing, at the cost of fewer valid pairs and more decorrelation; check how far the measured images follow that trend.

ISTA versus FISTA: function-value error over 10,000 iterations

A. Beck, M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," SIAM J. Imaging Sci. 2(1), 183–202, 2009. Fig. 5 (\(F(\mathbf{x}_k) - F(\mathbf{x}^*)\) for ISTA, MTWIST and FISTA, log scale). doi:10.1137/080716542

FISTA reaches about \(10^{-7}\) while ISTA stalls near \(10^{-3}\) on the same problem: the \(O(1/k^2)\) versus \(O(1/k)\) rates of (13.14) on their noiseless 64 × 64 deblurring problem (run with zero L1 weight, so plain least squares with \(F^* = 0\)), the same ordering as what m13_ista_sparse_spikes.m prints for sparse-spike deconvolution.

13.16Go deeper

P. C. Hansen, Discrete Inverse Problems: Insight and Algorithms, SIAM, 2010. Filter factors, Picard, L-curve, GCV, iterative methods; the convention used here.
P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, 1998. General-form Tikhonov and the GSVD.
H. W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer, 1996. Landweber, CG and semi-iterative methods, parameter choice theory.
J. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer, 2005. The Bayesian view and the inverse crime.
G. H. Golub, M. Heath, G. Wahba, "Generalized cross-validation as a method for choosing a good ridge parameter," Technometrics 21, 1979.
P. C. Hansen, D. P. O'Leary, "The use of the L-curve in the regularization of discrete ill-posed problems," SIAM J. Sci. Comput. 14, 1993. C. R. Vogel, "Non-convergence of the L-curve regularization parameter selection method," Inverse Problems 12, 1996.
I. Daubechies, M. Defrise, C. De Mol, "An iterative thresholding algorithm for linear inverse problems with a sparsity constraint," Comm. Pure Appl. Math. 57, 2004. A. Beck, M. Teboulle, "A fast iterative shrinkage-thresholding algorithm for linear inverse problems," SIAM J. Imaging Sci. 2, 2009.
D. L. Donoho, I. M. Johnstone, "Ideal spatial adaptation by wavelet shrinkage," Biometrika 81, 1994.
L. I. Rudin, S. Osher, E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D 60, 1992. S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," Found. Trends Mach. Learn. 3, 2011.
M. Jaeger, G. Held, S. Peeters, S. Preisser, M. Grünig, M. Frenz, "Computed ultrasound tomography in echo mode for imaging speed of sound using pulse-echo sonography: proof of principle," Ultrasound Med. Biol. 41, 2015. P. Stähli, M. Kuriakose, M. Frenz, M. Jaeger, "Improved forward model for quantitative pulse-echo speed-of-sound imaging," Ultrasonics 108, 2020.
S. J. Sanabria, E. Ozkan, M. Rominger, O. Goksel, "Spatial domain reconstruction for imaging speed-of-sound with pulse-echo ultrasound: simulation and in vivo study," Phys. Med. Biol. 63, 2018.
M. Jakovljevic, S. Hsieh, R. Ali, G. Chau Loo Kung, D. Hyun, J. J. Dahl, "Local speed of sound estimation in tissue using pulse-echo ultrasound: model-based approach," J. Acoust. Soc. Am. 144, 2018.
T. Szasz, A. Basarab, D. Kouamé, "Beamforming through regularized inverse problems in ultrasound medical imaging," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 63, 2016.

← 12 Estimation & CRLB 14 Optimization & autodiff →