9. Inverse Problems

Chapter PDF · Complete book

The main references for this chapter are [23, 29, 19].

9.1 Regularization of Inverse Problems¶

We extend the convex regularization introduced in Chapter 8 to linear measurement operators.

We consider a bounded linear map $\Phi : \mathcal{S}\rightarrow \mathcal{H}$, where the signal space $\mathcal{S}$ is a Hilbert space or, more generally, a Banach space. The data space $\mathcal{H}$ is a Hilbert space.

The operator models acquisition: an unknown high-resolution signal $f_0 \in \mathcal{S}$ gives rise to the noisy observation

where $w \in \mathcal{H}$ represents acquisition noise. Here the noise is deterministic; we assume only that $|\!| w |\!|_\mathcal{H}$ is bounded.

An acquisition device records finitely many observations, so applications usually take $\mathcal{H}=\mathbb{R}^P$. The number $P$ may be small relative to the desired reconstruction dimension.

For numerical implementation, we discretize the signal space as $\mathcal{S}=\mathbb{R}^N$, where $N$ is the number of grid points. Usually the grid is much finer than the measurement set, so $N \gg P$.

Infinite-dimensional function spaces remain useful for modeling the unknown $f_0$ and analyzing recovery, particularly when choosing the signal space $\mathcal{S}$.

Direct inversion may be impossible because $\Phi$ is not injective or the data do not lie in its range. Even when an inverse exists on the range, it may have a large norm or be unbounded. Applying it to noisy data amplifies the perturbation: formally, $\Phi^{-1}y=f_0+\Phi^{-1}w$.

We now give a few representative examples of forward operators $\Phi$.

Denoising.

The identity operator $\Phi=\mathrm{Id}_\mathcal{S}$ with $\mathcal{S}=\mathcal{H}$ gives the denoising problem studied in Chapters 7 and 8.

De-blurring and super-resolution.

For a general operator $\Phi$, recovering $f_0$ requires both inversion and denoising. These goals often conflict because inversion amplifies noise.

For example, in deblurring, $\Phi$ is a translation-invariant operator corresponding to low-pass filtering with a kernel $h$

A periodic model takes $\mathcal{S}=\mathcal{H}=L^2(\mathbb{T}^d)$; see Proposition 2.3.

In practice, the blurred signal is sampled on a grid: $\Phi f = \left\{ (f \star h)(x_k) \;;\; 0 \leqslant k< P \right\}$. Figure 9.1, middle, shows the resulting low-resolution image $\Phi f_0$.

Inverting such operators is useful for increasing the resolution of digital photographs and videos.

Interpolation and inpainting.

Inpainting reconstructs missing pixels. We model the masking operation by an operator that is diagonal in the spatial domain:

Here $\Omega$ is the missing region, either a subset of $[0,1]^d$ in the continuous model or a set of pixel indices in the discrete model.

Figure 9.1, right, shows an example of a damaged image $\Phi f_0$.

Medical imaging.

Many medical imaging devices provide indirect access to the signal of interest, with acquisition processes that can be approximated by a linear operator $\Phi$.

For a tomographic scanner, the Radon transform models acquisition. The Fourier slice theorem relates these measurements to Fourier data along radial lines.

An idealized magnetic resonance imaging (MRI) model also uses partial Fourier measurements:

Here $\Omega$ consists of radial lines for tomography or acquisition trajectories such as spirals for MRI.

Electroencephalography (EEG) and magnetoencephalography (MEG) also infer internal activity from measurements at sensors outside the sources. Abstractly, such a linearized model can be written as $(\Phi f)(s)=\int K(s,x)f(x)\,\mathrm{d}x$, where the kernel depends on the geometry and physical model. Unlike a simple image blur, this kernel need not be translation invariant.

Regression for supervised learning.

Linear supervised learning is closely related to the imaging problems studied here, although it uses different notation. Section 13.3 develops the connection between regression and inverse problems.

In statistical learning, the data are $n$ pairs $(x_i,y_i)_{i=1}^n$ with feature vectors $x_i \in \mathbb{R}^p$. A linear prediction model has the form $y_i = \langle \beta,\,x_i\rangle$, with unknown parameter $\beta\in \mathbb{R}^p$. Placing the vectors $x_i$ in the rows of $X \in \mathbb{R}^{n \times p}$ gives the approximate system $X \beta\approx y$. This corresponds to the inverse problem $\Phi f=y$ under the substitutions $\Phi \mapsto X$ and $f \mapsto \beta$, with dimensions $(P,N) \rightarrow (n,p)$.

In statistical learning, the model need not be correctly specified, and the design matrix $X$ is often random. Its sampling fluctuations must then be controlled as $n\to\infty$.

Ridge regression estimates the parameter by solving the normalized problem

The corresponding empirical quantities are therefore $\frac{1}{n}X^*X \sim \Phi^*\Phi$ for the covariance and $\frac{1}{n}X^*y \sim \Phi^*y$ for the data term.

For independent observations in fixed dimension and with suitable finite moments, these empirical quantities fluctuate around their expectations at order $n^{-1/2}$. In growing dimension, the relevant norm bounds also depend on the dimension and tail assumptions.

9.2 Warmup: Oracle Linear Inversion¶

We first seek an oracle estimator among inversion methods that act diagonally in a fixed basis, allowing the gains to depend on the unknown clean signal. The optimization then separates into one scalar problem per coefficient. Later we will obtain such estimators from a single optimization problem on the signal space.

Assume that the forward operator is diagonal in an orthonormal basis $(\psi_k)$:

where the coefficients $\varphi_k$ typically represent attenuation.

The observation model is

where $w$ is Gaussian white noise with variance $\sigma^2$. In a real orthonormal basis, its coefficients are independent centered Gaussians. In a complex Fourier basis they still have mean zero and second moment $\sigma^2$, which suffices for the risk calculation below. We work in finite dimension; infinite-dimensional white noise requires a generalized random-process interpretation.

For instance, if $(\psi_k)_k$ is the Fourier basis, this corresponds to a deconvolution problem $\Phi f = h \star f$ where $\hat h_k = \varphi_k$. We consider the same diagonal estimator $\tilde{f}$ as in the denoising case:

We minimize the mean squared error with respect to $\lambda$:

Writing $c_k=\langle f_0,\psi_k\rangle$ and using the zero mean and variance of the noise,

The optimal coefficient of the oracle estimator solves

i.e., it satisfies

which gives

For $\varphi_k \neq 0$ and $c_k\neq0$, the oracle gains approach direct inversion as the noise vanishes:

As $\sigma$ increases, inversion is attenuated. For a simple Fourier-domain example, let $k\geqslant 1$, $|c_k|^2=k^{-\alpha}$ and $\varphi_k=k^{-\beta/2}$, with $\alpha,\beta>0$. Then

At low frequencies this filter approximately compensates for attenuation, whereas at high frequencies it suppresses noise. Treating $k$ as a positive real variable, its maximum occurs at

On the discrete range $k\geqslant 1$, the maximum is attained near this value, or at the lowest frequency if $k^\star<1$.

9.3 Theoretical Study of Quadratic Regularization¶

We outline a standard approach to recovery guarantees in Hilbert spaces. The emphasis is on modeling assumptions expressed through source conditions, and on limitations of the squared-norm penalty, notably saturation of convergence rates and loss of components in $\ker(\Phi)$.

9.3.1 Singular Value Decomposition¶

Finite dimension.

We begin with the finite-dimensional case $\Phi \in \mathbb{R}^{P \times N}$ so that $\mathcal{S}=\mathbb{R}^N$ and $\mathcal{H}=\mathbb{R}^P$ are Hilbert spaces.

In this case, the singular value decomposition (SVD) provides a precise description of the operator and of linear inversion.

The theorem also holds for complex matrices, with $^\top$ replaced by $^*$.

Expression (9.4) describes $\Phi$ as a sum of rank-1 matrices $u_m v_m^\top$.

One usually orders the singular values $(\sigma_m)_m$ in decreasing order $\sigma_1 \geqslant\ldots \geqslant\sigma_R$. If the positive singular values are distinct, the real reduced SVD is unique up to simultaneous sign changes of each pair $(u_m,v_m)$. In the complex case, simultaneous multiplication by a unit-modulus scalar replaces the sign change.

The columns of $U$ form an orthonormal basis of $\mathop{\mathrm{Im}}(\Phi)$, while the columns of $V$ form an orthonormal basis of $\mathop{\mathrm{Im}}(\Phi^\top)=\ker(\Phi)^\bot$.

The decomposition (9.4) is called the reduced SVD because it retains only the $R$ nonzero singular values. The full SVD completes the columns of $U$ and $V$ to orthonormal bases of $\mathbb{R}^P$ and $\mathbb{R}^N$, respectively. The diagonal factor $\Sigma$ then has size $P\times N$.

For a periodic discrete convolution $\Phi f=h\star f$, let $F$ be the unitary discrete Fourier matrix. Then

where $\hat h_m$ denotes the convolution multiplier (the unnormalized DFT of $h$ for the usual discrete convolution). If $q_m=F^*e_m$ and $\hat h_m\neq0$, one may take $v_m=q_m$, $u_m=(\hat h_m/|\hat h_m|)q_m$, and $\sigma_m=|\hat h_m|$.

Computing the SVD of a dense matrix $\Phi \in \mathbb{R}^{N \times N}$ typically costs $O(N^3)$ operations.

Compact operators.

The SVD extends to compact operators $\Phi : \mathcal{S}\rightarrow \mathcal{H}$ between separable Hilbert spaces. Compactness means that $\Phi B_1$ is relatively compact, where $B_1= \left\{ s \in \mathcal{S} \;;\; |\!| s |\!|\leqslant 1 \right\}$ is the unit ball. Equivalently, every sequence $(\Phi s_k)_k$ with $s_k \in B_1$ has a convergent subsequence. In infinite dimension, the identity operator $\Phi : \mathcal{S}\rightarrow \mathcal{S}$ is not compact.

Equivalently, compact operators $\Phi$ admit an expansion analogous to (9.4):

Here $(u_m)$ and $(v_m)$ are orthonormal systems in $\mathcal{H}$ and $\mathcal{S}$, respectively, and the positive singular values decrease to zero; finite-rank operators give a finite sum.

Convergence in (9.6) holds in the operator norm induced by the norms on $\mathcal{S}$ and $\mathcal{H}$

For a nonzero operator $\Phi$ with the expansion (9.6), $|\!| \Phi |\!|_{\mathcal{L}(\mathcal{S},\mathcal{H})}=\sigma_1$.

If only $R$ singular values are positive, then $\Phi$ has finite rank $R=\dim(\mathop{\mathrm{Im}}(\Phi))$. Squared-norm regularization restricts the reconstruction to $\ker(\Phi)^\bot$, a space of dimension $R$, so the problem reduces to finite dimension. Nonlinear methods can also recover components in $\ker(\Phi)$ under suitable priors; this is sometimes called super-resolution.

Between Hilbert spaces, compact operators are exactly the operator-norm limits of finite-rank operators. This approximation property need not hold for arbitrary Banach spaces, although the definition of compactness itself extends to them.

Examples include integral operators on $L^2(\Omega)$ with a continuous kernel $k(x,y)$ defined for $(x,y) \in \Omega\times\Omega$, where $\Omega$ is a compact subset of $\mathbb{R}^d$ or the torus $\mathbb{T}^d$:

where $\mathrm{d}y$ is the Lebesgue measure.

Convolution $\Phi f = f \star h$ on $\mathbb{T}^d=(\mathbb{R}/2\pi\mathbb{Z})^d$ generalizes (9.5). Its kernel $k(x,y) = h(x-y)$ is translation invariant, and the Fourier functions $q_m(x)=(2\pi)^{-d/2}e^{\mathrm{i}m\cdot x}$ diagonalize the operator. Its singular values are $\sigma_m=|\hat h_m|$, where $\hat h_m=\int_{\mathbb{T}^d}h(x)e^{-\mathrm{i}m\cdot x}\,\mathrm{d}x$ denotes the convolution multiplier.

Another example on $\Omega=[0,1]$ is the integration operator $(\Phi f)(x) = \int_{0}^x f(y) \mathrm{d}y$, which has the square-integrable kernel $k(x,y)=1_{y\leqslant x}$ and is also compact.

Pseudoinverse.

We first work in finite dimension, or more generally assume that $\mathop{\mathrm{Im}}(\Phi)$ is closed. Even without noise, $\Phi f=y$ may be inconsistent, and a nontrivial kernel makes a solution nonunique. The Moore–Penrose pseudoinverse selects the least-squares solution of minimum norm

The next proposition expresses the pseudoinverse through the SVD. Under injectivity or surjectivity, it also gives formulas involving linear systems with $\Phi\Phi^*$ or $\Phi^*\Phi$.

In infinite dimension, a compact operator of infinite rank has nonclosed range. Its pseudoinverse is unbounded and is defined only for data satisfying the corresponding square-summability condition on the inverse-weighted SVD coefficients. The projections and inverse formulas above cannot then be used on arbitrary data.

For convolution operators $\Phi f = f\star h$, the spectral formula gives

9.3.2 Tikhonov Regularization¶

Regularized inverse.

For an ill-conditioned operator, applying the unregularized inverse to noisy data can be unstable, since

The recovery error is therefore $|\!| \Phi^+ y -f_0^+ |\!| = |\!| \Phi^+ w |\!|$. It reaches $|\!| w |\!|/\sigma_R$ when $w \propto u_R$, so a small singular value produces a large amplification factor $1/\sigma_R$. In the infinite-rank case $R=+\infty$, the compact operator’s pseudoinverse is unbounded. To control this amplification, replace $\Phi^+$ by a regularized inverse of the form

where the spectral filter $\mu_\lambda$ depends on a parameter $\lambda>0$ and satisfies the boundedness and consistency conditions

Figure 9.2, left, illustrates a capped reciprocal filter, which limits the amplification of small singular values.

Variational regularization.

A typical regularized inverse is obtained by solving a penalized least-squares problem

where $J$ is a regularization functional. For general convex priors, lower semicontinuity, rather than continuity, is the natural assumption; existence also requires compactness or coercivity conditions. The simplest example is the quadratic norm $J=|\!| \cdot |\!|_\mathcal{S}^2$,

Proposition 9.3 shows that this estimator is a special case of (9.8). Its normal equations give

This shows that $f_\lambda\in \mathop{\mathrm{Im}}(\Phi^*)\subset\ker(\Phi)^\bot$, and that it depends linearly on $y$.

For convolution $\Phi f = f \star h$, the FFT computes the regularized inverse in $O(N \log(N))$ operations:

Figure 9.2 compares quadratic regularization (9.11) (right) with a capped reciprocal filter (left).

We choose $\lambda$ according to the noise level and seek a rate for $f_\lambda\to f_0$. This requires $f_0=f_0^+$, equivalently $f_0\in\overline{\mathop{\mathrm{Im}}(\Phi^*)}=\ker(\Phi)^\bot$: squared-norm regularization always produces $f_\lambda\in\mathop{\mathrm{Im}}(\Phi^*)\subset\ker(\Phi)^\bot$, so it cannot recover a kernel component of $f_0$. Nonlinear priors can recover such components when additional structure makes the signal identifiable.

Source condition.

To obtain a convergence rate, we impose a source condition of order $\beta$, which reads

For a fixed bound on the source vector, larger $\beta$ imposes stronger decay along directions with small singular values and makes stable inversion easier. The condition asserts the existence of $z\in\mathcal{S}$ such that $f_0 = V \mathop{\mathrm{diag}}( \sigma_m^{\beta} ) V^* z$. Choose the source vector of minimum norm, $z=V\mathop{\mathrm{diag}}(\sigma_m^{-\beta})V^*f_0$, and impose $|\!| z |\!| \leqslant\rho$ for some $\rho>0$. This is equivalent to the coefficient bound

The assumptions $\beta>0$ and $f_0\in\ker(\Phi)^\bot$ are part of the source condition; the coefficient bound alone leaves kernel components unconstrained. The bound resembles the Sobolev constraint in 7.6. For example, consider a low-pass convolution $\Phi f = f \star h$ where the kernel $h$ has a polynomially decaying multiplier: $|\hat h_m| \sim 1/m^{\alpha}$ for large $m$. Since the vectors $v_m$ are Fourier modes, ($S_{\beta ,\rho }$) gives, up to constants and the Fourier normalization, a bound of the form

This is a Sobolev-type coefficient bound of smoothness order $\alpha\beta$.

Sublinear convergence speed.

The following theorem gives the convergence rate implied by this source condition.

With the noise bound $|\!| w |\!| \leqslant\delta$, recovery depends on the parameters $(\delta,\rho,\beta)$. Assuming $f_0 \in \ker(\Phi)^\bot$, we study the convergence of $f_\lambda$ to $f_0$ for data $y=\Phi f_0+w$ as $\delta\rightarrow 0$. The analysis also determines how $\lambda$ should depend on $\delta$.

Larger source orders $\beta\leqslant 2$ therefore give faster convergence as $|\!| w |\!|$ tends to zero. The rate in Theorem 9.4 saturates: choosing $\beta> 2$ gives the same general worst-case rate as $\beta=2$. The best rate obtainable in this way is

Alternative spectral filters $\mu_\lambda$, together with sufficiently large $\beta$, can give the rate $|\!| f_\lambda- f_0 |\!| = O( \delta^{1-\kappa} )$ for arbitrarily small $\kappa>0$.

The capped inverse in Figure 9.2, left, avoids this finite-order saturation, whereas Tikhonov regularization is limited to source orders $\beta\leqslant 2$. Quadratic regularization is easier to implement: its variational formulation leads to a linear system and does not require an SVD.

For compact operators with nonclosed range, no uniform linear rate follows over these source balls from finite-order source conditions alone. Linear rates are possible in finite dimension with a fixed smallest positive singular value, or under different structural assumptions, such as those used for sparse $\ell^1$ regularization in Chapter 10.

9.4 Quadratic Regularization¶

We now turn from infinite-dimensional recovery theory to numerical methods in finite dimension.

Convex regularization.

As in (9.9), we reconstruct a high-resolution image $f_0 \in \mathbb{R}^N$ from noisy measurements $y = \Phi f_0 + w \in \mathbb{R}^P$ by minimizing a convex regularized objective:

Here $|\!| y-\Phi f |\!|^2$ measures data fidelity, with $|\!| \cdot |\!|$ denoting the $\ell^2$ norm on $\mathbb{R}^P$, and $J(f)$ is a convex penalty on $\mathbb{R}^N$.

The regularization parameter $\lambda>0$ balances these two terms and can be difficult to choose in practice. In simulations with a known reference signal $f_0$, the parameter can be calibrated by minimizing the reconstruction error $|\!| f_0-\tilde f |\!|$. Such a reference is generally unavailable in applications.

For noiseless data, $w = 0$, we study small values of $\lambda$. As $\lambda\rightarrow 0$, minimizers of (9.17) approach solutions of a constrained problem under the assumptions below.

We may assume $y \in \mathop{\mathrm{Im}}(\Phi)$ without changing the minimizers, since orthogonality gives

The first term is independent of $f$, so replacing $y$ by $\mathop{\mathrm{Proj}}_{\mathop{\mathrm{Im}}(\Phi)}(y)$ in (9.17) changes only an additive constant.

Let us recall that a function $J$ is coercive if

i.e.

Equivalently, every sublevel set $\left\{ f \;;\; J(f) \leqslant c \right\}$ is bounded. In finite dimension, lower semicontinuity of $J$ then makes these sets compact.

Full coercivity of $J$ can be replaced by boundedness of joint sublevel sets of $J(f)$ and $|\!| \Phi f |\!|$. For seminorm penalties such as TV, this holds when $\ker(\Phi)$ intersects the nullspace of the seminorm only at zero; in particular, the measurements must detect constant images.

Quadratic Regularization.

The simplest prior functionals are quadratic, and can be written as

where $G \in \mathbb{R}^{K \times N}$ and $L=G^*G \in \mathbb{R}^{N \times N}$ is a positive semidefinite matrix.

The special case (9.10) is recovered when setting $G=L=\mathrm{Id}_N$.

Writing down the first-order optimality conditions for (9.17) leads to

hence, if

then (9.17) has a unique minimizer $f_\lambda$, which is obtained by solving a linear system

If $L$ is diagonal in a full right-singular basis of $\Phi$, with eigenvalues $\alpha_m^2$, then, for the positive singular values,

Example of convolution.

Consider the convolution operator

and $G=\nabla$, a discretization of the gradient operator, for example using first-order finite differences (2.16). This corresponds to the discrete Sobolev prior introduced in Section 8.1.2.

Such an operator computes, for a $d$-dimensional signal $f \in \mathbb{R}^N$ (for instance a 1-D signal for $d=1$ or an image when $d=2$), an approximation $\nabla f_n \in \mathbb{R}^d$ of the gradient vector at each sample location $n$. Thus typically, $\nabla : f \mapsto (\nabla f_n)_n \in \mathbb{R}^{N \times d}$ maps to $d$-dimensional vector fields.

Then $-\nabla^* : \mathbb{R}^{N \times d} \rightarrow \mathbb{R}^N$ is a discretized divergence operator.

In this case, $\Delta=-G^*G$ is a discretization of the Laplacian, which is itself a convolution operator. One then has

where $\hat d_2$ is the Fourier transform of the filter $d_2$ corresponding to the Laplacian. For instance, in dimension 1, using first-order finite differences, the expression for $\hat d_{2,m}$ is given in (2.18).

9.4.1 Solving Linear Systems¶

If $\Phi$ and $L$ are not simultaneously diagonalizable in the required bases, the coefficient formula (9.21) is unavailable. We instead solve (9.20), written as

For moderate $N$, up to a few thousand, direct methods solve this system at a cost of $O(N^3)$ for a generic $A$. Since $A$ is symmetric positive definite, Cholesky factorization $A=BB^*$ is a natural choice, with $B$ lower triangular. Reordering the rows and columns can help preserve sparsity in the factors of $A$.

For large $N$, iterative solvers access $A$ only through matrix-vector products. In imaging, operator structure often makes these products much cheaper than a dense $O(N^2)$ computation. For a sufficiently sparse matrix $A$, the cost is typically $O(N)$.

When $A$ is symmetric positive definite, as assumed here, the conjugate gradient method solves the equivalent quadratic minimization problem

This objective agrees with (9.17) up to a constant. Starting from any $f^{(0)}$, conjugate gradients computes $f^{(\ell+1)}$ by minimizing over the affine space generated by all gradients obtained so far, rather than taking a single gradient step as in Section 19.1:

This subspace minimization requires only one matrix-vector product per iteration. For $\ell\geqslant 1$, with initialization $d^{(0)} = \nabla \mathcal{E}(f^{(0)}) = Af^{(0)}-b$, the recurrence is

where $g^{(\ell)}\mathrel{:=}\nabla\mathcal{E}(f^{(\ell)})=Af^{(\ell)}-b$. The iteration stops as soon as $g^{(\ell)}=0$.

It can also be shown that the directions $d^{(\ell)}$ are $A$-orthogonal (conjugate), so that after at most $N$ iterations in exact arithmetic, the conjugate gradient method computes the unique solution $f^{(\ell)}$ of the linear system $A f = b$. In practice, it is usually used as an approximate solver, with far fewer than $N$ iterations.

Nonlinear conjugate-gradient methods extend this idea to smooth objectives, but require a line search and additional safeguards. The linear-system step-size formula and finite-termination property above do not carry over unchanged.

9.5 Non-Quadratic Regularization¶

9.5.1 Total Variation Regularization¶

Quadratic regularization (9.19) often produces a blurred reconstruction $f_\lambda$. This can poorly approximate a signal $f_0$ with sharp transitions, such as image edges.

In the convolutional case (9.23), the restoration operator is a filter, which tends to smooth sharp features.

For the squared-norm penalty, the restored signal belongs to $\ker(\Phi)^\bot$, so missing components cannot be recovered. This statement does not hold for every quadratic penalty: a different matrix $L$ can couple observed and unobserved components. Nonquadratic priors, particularly nonsmooth ones such as TV and $\ell^1$, encode structures such as edges and sparse coefficients that quadratic penalties do not model well.

Total variation.

Total variation is a widely used nonquadratic, nonsmooth prior. For smooth functions $f : \mathbb{R}^d \mapsto \mathbb{R}$, it replaces the Sobolev or Dirichlet energy

where $\nabla f(x) = (\partial_{x_1} f(x), \ldots, \partial_{x_d} f(x))^\top$ is the gradient, by the (vectorial) $L^1$ norm of the gradient

Section 8.1.1 gives further background on these priors.

Removing the square $^2$ from the integrand changes which functions have finite energy.

For a nontrivial bounded set $\Omega$ with a smooth boundary, $J_{\text{Sob}}(1_\Omega)=+\infty$, whereas $J_{\text{TV}}(1_\Omega)=\operatorname{Per}(\Omega)<+\infty$. More generally, an integrable function $f$ is of bounded variation when its distributional gradient $Df$ is a finite vector-valued Radon measure. Its total variation is

For a finite vector measure $m$, its total mass is equivalently

The coarea formula for $f\in\mathrm{BV}(\mathbb{R}^d)$ reads

For smooth functions, this agrees, for almost every level, with integrating the $(d-1)$-dimensional Hausdorff measure of the level sets. For general BV functions, superlevel-set perimeters, rather than measures of the pointwise sets $\{f=t\}$, are required.

Discretized Total variation.

For discrete data $f \in \mathbb{R}^N$, the construction in Section 8.1.2 generalizes (8.6) to any spatial dimension:

where $\nabla f_n \in \mathbb{R}^d$ approximates the spatial gradient at grid point $n$.

The discrete total variation prior $J_{\text{TV}}(f)$ defined in (8.6) is a convex but nondifferentiable function of $f$, since a term of the form $|\!| \nabla f_n |\!|$ is nondifferentiable if $\nabla f_n=0$.

Chapters 18 and 19 develop nonsmooth convex optimization methods for such functionals.

To use classical gradient descent, we replace the nonsmooth $\ell^2$ norm $|\!| \cdot |\!|$ by a differentiable approximation, for example

This gives the smoothed TV functional already introduced in (8.12):

For fixed $u$, the limits as $\varepsilon\downarrow0$ and $\varepsilon\to+\infty$ are described by

Thus $J_{\text{TV}}^\varepsilon$ connects $J_{\text{TV}}$ to a rescaled version of $J_{\text{Sob}}$, up to an additive constant.

The resulting regularized inverse problem (9.17) thus reads

If $\ker(\Phi)\cap\ker(\nabla)=\{0\}$, the objective is coercive and strictly convex, so it has a unique minimizer. Strict convexity of $|\!| \cdot |\!|_\varepsilon$ alone does not ensure uniqueness after composition with the gradient, whose nullspace contains constants.

9.5.2 Gradient Descent Method¶

The optimization program (9.26) is an example of smooth unconstrained convex optimization of the form

where $\mathcal{E}: \mathbb{R}^N \rightarrow \mathbb{R}$ is a $\mathcal{C}^1$ function. Recall that the gradient $\nabla \mathcal{E}: \mathbb{R}^N \mapsto \mathbb{R}^N$ of this functional (not to be confused with the discretized gradient $\nabla f \in \mathbb{R}^{N\times d}$ of $f$) is defined by the following first-order relation

A locally Lipschitz gradient improves the remainder to $O(|\!| r |\!|^2)$. Mere continuity of the gradient does not give this stronger estimate.

For such a function, the gradient descent algorithm is defined as

The step size $\tau_\ell>0$ must balance stability against progress per iteration.

Section 19.1 analyzes convergence and explains how it depends on the step size $\tau_\ell$.

9.5.3 Examples of Gradient Computation¶

Note that the gradient of a quadratic function $\mathcal{Q}(f)$ of the form (9.24) reads

In particular, the first-order optimality condition $\nabla \mathcal{Q}(f)=0$ is equivalent to the linear system $Af=b$.

For the quadratic fidelity term $\mathcal{G}(f)=\frac{1}{2}|\!| \Phi f-y |\!|^2$, one thus obtains

In the special case of the regularized TV problem (9.26), the gradient of $\mathcal{E}$ reads

For differentiable maps $\mathcal{F}:\mathbb{R}^P\to\mathbb{R}$ and $\mathcal{H}:\mathbb{R}^N\to\mathbb{R}^P$, the chain rule gives

Apply it to $J_{\text{TV}}^\varepsilon=|\!| \cdot |\!|_{1,\varepsilon}\circ\nabla$, where $|\!| u |\!|_{1,\varepsilon}=\sum_n\sqrt{\varepsilon^2+|\!| u_n |\!|^2}$. The result is

Thus the gradient of the energy is computed by taking finite differences, normalizing the resulting vector field, and applying minus the divergence.

Since $\mathop{\mathrm{div}}=-\nabla^*$, their operator norms are equal. For unscaled forward differences in $d$ dimensions with periodic boundaries, $|\!| \nabla |\!|_{\mathrm{op}}\leqslant 2\sqrt d$. The Hessian of $u\mapsto\sqrt{\varepsilon^2+|\!| u |\!|^2}$ is