2.2.1 Image denoising

The observed data are the ideal solution corrupted by additive noise, so the forward operator in (2.1) is the identity transform ${\mathcal{A}}=\text{id}$ , and we get

(2.2) $$\begin{eqnarray}g=f_{\text{true}}+e,\end{eqnarray}$$

In the simplest case the distribution of the observational noise is known. Furthermore, this distribution may in more advanced problems be correlated, spatially varying and of mixed type.

In Section 7.2 we will discuss bilevel learning of total variation (TV)-type variational models for denoising of data corrupted with mixed noise distributions.

2.2.2 Image deblurring

The observed data are given by convolution with a known filter function $K$ together with additive noise, so (2.1) becomes

(2.3) $$\begin{eqnarray}g=f_{\text{true}}\ast K+e.\end{eqnarray}$$

Any inverse problem of the type (2.1) with a linear forward operator that is translation-invariant will be of this form.

In the absence of noise, the inverse problem (image deconvolution) is exactly solvable by division in the Fourier domain, i.e.  $f_{\text{true}}={\mathcal{F}}^{-1}[{\mathcal{F}}[g]/{\mathcal{F}}[K]]$ , provided that ${\mathcal{F}}[K]$ has infinite support in the Fourier domain. In the presence of noise, the estimated solution is corrupted by noise whose frequency spectrum is the reciprocal of the spectrum of the filter $K$ . The distribution of the observational also has the same considerations as in (2.2). Finally, extensions include the case of a spatially varying kernel and the case where $K$ is unknown (blind deconvolution).

2.2.3 Image in-painting

Here, the observed data represents a noisy observation of the true model parameter $f_{\text{true}}:\unicode[STIX]{x1D6FA}\rightarrow \mathbb{R}$ restricted to a fixed measurable set $\unicode[STIX]{x1D6FA}_{0}\subset \mathbb{R}^{n}$ :

(2.4) $$\begin{eqnarray}g=f_{\text{true}}\big|_{\unicode[STIX]{x1D6FA}_{0}}+e.\end{eqnarray}$$

In the above, $f_{\text{true}}\big|_{\unicode[STIX]{x1D6FA}_{0}}$ is the restriction of $f_{\text{true}}$ to $\unicode[STIX]{x1D6FA}_{0}$ . Solutions take different forms depending on the size of connected components in $\unicode[STIX]{x1D6FA}_{0}$ . Extensions include the case where $\unicode[STIX]{x1D6FA}_{0}$ is unknown or only partially known.

2.2.4 Computed tomography (CT)

The simplest physical model for CT assumes mono-energetic X-rays and disregards scattering phenomena. The model parameter is then a real-valued function $f:\unicode[STIX]{x1D6FA}\rightarrow \mathbb{R}$ defined on a fixed domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ ( $d=2$ for two-dimensional CT and $d=3$ for three-dimensional CT) that has unit mass per volume. The forward operator is the one given by the Beer–Lambert law:

(2.5) $$\begin{eqnarray}{\mathcal{A}}(f)(\unicode[STIX]{x1D714},x)=\text{e}^{-\unicode[STIX]{x1D707}\int _{-\infty }^{\infty }f(x+s\unicode[STIX]{x1D714})\,\text{d}s}.\end{eqnarray}$$

Here, the unit vector $\unicode[STIX]{x1D714}\in S^{d-1}$ and $x\in \unicode[STIX]{x1D714}^{\bot }$ represent the line $\ell :s\mapsto x+s\unicode[STIX]{x1D714}$ along which the X-rays travel, and we also assume $f$ decays fast enough for the integral to exist. In medical imaging, $\unicode[STIX]{x1D707}$ is usually set to a value that approximately corresponds to water at the X-ray energies used. The above represents pre-logarithm (or pre-log) data, and by taking the logarithm (or log) of data, one can recast the inverse problem in CT imaging to one where the forward model is the linear ray transform:

(2.6)

For low-dose imaging, pre-log data are Poisson-distributed with mean ${\mathcal{A}}(f_{\text{true}})$ , where ${\mathcal{A}}$ is given as in (2.5), that is, $g\in Y$ is a sample of $\unicode[STIX]{x1D558}\sim \operatorname{Poisson}({\mathcal{A}}(f_{\text{true}}))$ . Thus, to get rid of the non-linear exponential in (2.5), it is common to take the log of data. With such post-log data the forward operator is linear and given as in (2.6). A complication with such post-log data is that the noise model becomes non-trivial, since one takes the log of a Poisson-distributed random variable (Fu et al. Reference Fu, Lee, Kim, Alessio, Kinahan, Chang, Sauer, Kalra and Man2017). A common approximate noise model for post-log data is (2.1), with observational noise $e$ which is a sample of a Gaussian or Laplace-distributed random variable.

In the case of complete data, that is, where a full angular set of data is measured, an exact inverse is obtained by the (Fourier-transformed) data backprojected on the same lines as used for the measurements and scaled by the absolute value of the spatial frequency, followed by the inverse Fourier transform. Thus, as in deblurring, the noise is amplified, but only linearly in spatial frequency, making the problem mildly ill-posed. Extensions include the emission tomography problem (single photon emission computed tomography (SPECT) and positron emission tomography (PET)) where the line integrals are exponentially attenuated by a function $\unicode[STIX]{x1D707}$ that may be unknown. A major challenge in tomography is to consider incomplete data, in particular the case where only a subset of lines is measured. This problem is much more ill-posed.

See Sections 7.3, 7.4 and 7.6 for instances of CT reconstruction that use deep neural networks in the solution of the inverse problem.

2.2.5 Magnetic resonance imaging (MRI)

The observed data are often considered to be samples of the Fourier transform of the ideal signal, so the MRI image reconstruction problem is an inverse problem of the type (2.1), where the forward operator is given as a discrete sampling operator concatenated with the Fourier transform. A correct description of the problem takes account of the complex-valued nature of the data, which implies that when $e$ is normally distributed then the noise model of $|{\mathcal{F}}^{-1}[g]|$ is Rician. As in CT, the case of under-sampled data is of high practical importance. In MRI, the subsampling operator has to consist of connected trajectories in Fourier space but is not restricted to straight lines.

Extensions include the case of parallel MRI where the Extensions include the case of parallel MRI where the forward operator is combined with (several) spatial sensitivity functions. More exact forward operators take account of other non-linear physical effects and can reconstruct several functions in the solution space.

The data discrepancy.

Here the choice is ideally guided by statistical considerations for the observation noise (Bertero, Lantéri and Zanni Reference Bertero, Lantéri, Zanni and Censor2008). Ideally one selects ${\mathcal{L}}$ as an appropriate affine transform of the negative log-likelihood of data, in which case minimizing $f\mapsto {\mathcal{L}}({\mathcal{A}}(f),g)$ becomes the same as computing an maximum likelihood estimator. Hence, Poisson-distributed data that typically appear in photography (Costantini and Susstrunk Reference Costantini and Susstrunk2004) and emission tomography applications (Vardi, Shepp and Kaufman Reference Vardi, Shepp and Kaufman1985) lead to a data discrepancy given by the Kullback–Leibler divergence (Sawatzky, Brune, Müller and Burger Reference Sawatzky, Brune, Müller, Burger, Jiang and Petkov2009, Hohage and Werner Reference Hohage and Werner2016), while additive normally distributed data, as for Gaussian noise, result in a least-squares fit model.

The regularizer.

As stated in Definition 2.2, ${\mathcal{S}}$ acts as a penalizer and is chosen to enforce stability by encoding a priori information about $f_{\text{true}}$ . How to set the (regularization) parameter $\unicode[STIX]{x1D703}$ reflects noise level in data: see Section 4.1.

Classical Tikhonov regularization (2.10) uses Hilbert-space norms (or semi-norms) to regularize the inverse problems. In more recent years, Banach-space regularizers have become more popular in the context of sparsity-promoting and discontinuity-preserving regularization, which are revisited in Section 2.7. TV regularization was introduced by Rudin, Osher and Fatemi (Reference Rudin, Osher and Fatemi1992) for image denoising due to its edge-preserving properties, favouring images $f$ that have a sparse gradient. Here, the TV regularizer is given as

(2.12)

where $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ is a fixed open and bounded set. The above functional (TV regularizer) uses the total variation measure of the distributional derivative of $f$ defined on $\unicode[STIX]{x1D6FA}$ (Ambrosio, Fusco and Pallara Reference Ambrosio, Fusco and Pallara2000). A drawback of using such a regularization procedure is apparent as soon as the true model parameter not only consists of constant regions and jumps but also possesses more complicated, higher-order structures, e.g. piecewise linear parts. In this case, TV introduces jumps that are not present in the true solution, which is referred to as staircasing (Ring Reference Ring2000). Examples of generalizations of TV for addressing this drawback typically incorporate higher-order derivatives, e.g. total generalized variation (TGV) (Bredies, Kunisch and Pock Reference Bredies, Kunisch and Pock2011) and the infimal-convolution total variation (ICTV) model (Chambolle and Lions Reference Chambolle and Lions1997). These read as

(2.13) $$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{S}}_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}(f) :=\operatorname{ICTV}_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}(f)\nonumber\\ \displaystyle & & \displaystyle \quad =\min _{\substack{ v\in W^{1,1}(\unicode[STIX]{x1D6FA}) \\ \unicode[STIX]{x1D6FB}v\in BV(\unicode[STIX]{x1D6FA})}}\{\unicode[STIX]{x1D6FC}\Vert Df-\unicode[STIX]{x1D6FB}v\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{2})}+\unicode[STIX]{x1D6FD}\Vert D\unicode[STIX]{x1D6FB}v\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{2\times 2})}\},\end{eqnarray}$$

and the second-order TGV (Bredies and Valkonen Reference Bredies and Valkonen2011, Bredies, Kunisch and Valkonen Reference Bredies, Kunisch and Valkonen2013) reads as

(2.14) $$\begin{eqnarray}\displaystyle & & \displaystyle {\mathcal{S}}_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}(f) :=\operatorname{TGV}_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}^{2}(f)\nonumber\\ \displaystyle & & \displaystyle \quad =\min _{w\in B\!D(\unicode[STIX]{x1D6FA})}\{\unicode[STIX]{x1D6FC}\Vert Df-w\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{2})}+\unicode[STIX]{x1D6FD}\Vert Ew\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\operatorname{Sym}^{2}(\mathbb{R}^{2}))}\}.\end{eqnarray}$$

Here

$$\begin{eqnarray}B\!D(\unicode[STIX]{x1D6FA}) :=\{w\in L^{1}(\unicode[STIX]{x1D6FA};\mathbb{R}^{d})\mid \Vert Ew\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{d\times d})}<\infty \}\end{eqnarray}$$

is the space of vector fields of bounded deformation on $\unicode[STIX]{x1D6FA}$ with $E$ denoting the symmetrized gradient and $\text{Sym}^{2}(\mathbb{R}^{2})$ the space of symmetric tensors of order $2$ with arguments in $\mathbb{R}^{2}$ . The parameters $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}$ are fixed positive parameters. The main difference between (2.13) and (2.14) is that we do not generally have $w=\unicode[STIX]{x1D6FB}v$ for any function $v$ . That results in some qualitative differences of ICTV and TGV regularization: see e.g. Benning, Brune, Burger and Müller (Reference Benning, Brune, Burger and Müller2013). One may also consider Banach-space norms other than TV, such as Besov norms (Lassas, Saksman and Siltanen Reference Lassas, Saksman and Siltanen2009), which behave more nicely with respect to discretization (see also Section 3.4). Different TV-type regularizers and their adaption to data by bilevel learning of parameters (e.g.  $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ in ICTV and TGV) will be discussed in more detail in Section 4.3.1 and numerical results will be given in Section 7.2.

Finally, in applied harmonic analysis, $\ell _{p}$ -norms of wavelets have been proposed as regularizers (Daubechies et al. Reference Daubechies1991, Mallat Reference Mallat2009, Unser and Blu Reference Unser and Blu2000, Kutyniok and Labate Reference Kutyniok and Labate2012, Eldar and Kutyniok Reference Eldar and Kutyniok2012, Foucart and Rauhut Reference Foucart and Rauhut2013). Other examples are non-local regularization (Gilboa and Osher Reference Gilboa and Osher2008, Buades, Coll and Morel Reference Buades, Coll and Morel2005), anisotropic regularizers (Weickert Reference Weickert1998) and, in the context of free discontinuity problems, the representation of images as a composition of smooth parts separated by edges (Blake and Zisserman Reference Blake and Zisserman1987, Mumford and Shah Reference Mumford and Shah1989, Carriero, Leaci and Tomarelli Reference Carriero, Leaci, Tomarelli, Serapioni and Tomarelli1996).

2.7.1 Notions of sparsity

Let $X$ be a separable Hilbert space, that is, we will assume that it has a countable orthonormal basis. A popular approach to sparse reconstruction uses the notion of a dictionary $\mathbb{D} :=\{\unicode[STIX]{x1D719}_{i}\}\subset X$ , whose elements are called atoms. Here, $\mathbb{D}$ is either given, i.e. knowledge-driven, or data-driven and derived from a set of realizations $f_{i}\in X$ : see Bruckstein, Donoho and Elad (Reference Bruckstein, Donoho and Elad2009), Daubechies et al. (Reference Daubechies, Defrise and De Mol2004), Rubinstein, Bruckstein and Elad (Reference Rubinstein, Bruckstein and Elad.2010), Lanusse, Starck, Woiselle and Fadili (Reference Lanusse, Starck, Woiselle and Fadili2014) and Chen and Needell (Reference Chen and Needell2016).

A special class of dictionaries is that of frames. A dictionary $\mathbb{D} :=\{\unicode[STIX]{x1D719}_{i}\}$ is a frame if there exists $C_{1},C_{2}>0$ such that

(2.15) $$\begin{eqnarray}C_{1}\Vert f\Vert ^{2}\leq \mathop{\sum }_{i}|\langle f,\unicode[STIX]{x1D719}_{i}\rangle |^{2}\leq C_{2}\Vert f\Vert ^{2}\quad \text{for any }f\in X.\end{eqnarray}$$

A frame is called tight if $C_{1}=C_{2}=1$ ; it is called over-complete or redundant if $\mathbb{D}$ does not form a basis for $X$ . Redundant dictionaries, e.g. translation-invariant wavelets, often work better than non-redundant dictionaries: see Peyré, Bougleux and Cohen (Reference Peyré, Bougleux and Cohen2011) and Elad (Reference Elad2010).

To construct sparse representations, i.e. sparsity-promoting regularizers from such a parametrization, we need the notions of an analysis and a synthesis operator. Given a dictionary $\mathbb{D}$ , the analysis operator ${\mathcal{E}}_{\mathbb{D}}:X\rightarrow \unicode[STIX]{x1D6EF}$ maps an element in $X$ to a sequence in $\unicode[STIX]{x1D6EF}$ , which is typically in $\ell ^{2}$ , such that

$$\begin{eqnarray}\boldsymbol{\unicode[STIX]{x1D709}}={\mathcal{E}}_{\mathbb{D}}(f)\in \unicode[STIX]{x1D6EF}~\text{has components }\unicode[STIX]{x1D709}_{i}=\langle f,\unicode[STIX]{x1D719}_{i}\rangle .\end{eqnarray}$$

The corresponding synthesis operator ${\mathcal{E}}_{\mathbb{D}}^{\ast }:\unicode[STIX]{x1D6EF}\rightarrow X$ is the adjoint of the analysis operator, that is,

(2.16) $$\begin{eqnarray}{\mathcal{E}}_{\mathbb{D}}^{\ast }(\boldsymbol{\unicode[STIX]{x1D709}}) :=\mathop{\sum }_{i}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D719}_{i}.\end{eqnarray}$$

We further define the frame operator as ${\mathcal{E}}_{\mathbb{D}}^{\ast }\circ {\mathcal{E}}_{\mathbb{D}}:X\rightarrow X$ , that is,

(2.17) $$\begin{eqnarray}{\mathcal{E}}_{\mathbb{D}}^{\ast }\circ {\mathcal{E}}_{\mathbb{D}}(f) :=\mathop{\sum }_{i}\langle f,\unicode[STIX]{x1D719}_{i}\rangle \unicode[STIX]{x1D719}_{i}.\end{eqnarray}$$

Now, we can define an $f\in X$ to be $s$ -sparse with respect to $\mathbb{D}$ if

(2.18) $$\begin{eqnarray}\Vert {\mathcal{E}}_{\mathbb{D}}(f)\Vert _{0}=\#\{i|\langle f,\unicode[STIX]{x1D719}_{i}\rangle \neq 0\}\leq s.\end{eqnarray}$$

In most applications, the model parameter is not sparse in the strict sense, which leads to the weaker notion of compressibility. A model parameter $f\in X$ is compressible with respect to $\mathbb{D}$ if the following power decay law holds:

(2.19)

Here, $\widetilde{\boldsymbol{\unicode[STIX]{x1D709}}}$ is a non-increasing rearrangement of the sequence $\boldsymbol{\unicode[STIX]{x1D709}}=\{\langle f,\unicode[STIX]{x1D719}_{i}\rangle \}={\mathcal{E}}_{\mathbb{D}}(f)$ . Note that sparse signals are compressible, and in particular, if $q$ is small in (2.19), then compressibility is equivalent to sparsity from any practical viewpoint. We further define $\boldsymbol{\unicode[STIX]{x1D709}}_{s}$ to consist of the $s$ largest (in magnitude) coefficients of the sequence $\widetilde{\boldsymbol{\unicode[STIX]{x1D709}}}$ .

2.7.2 Sparse recovery

Here we consider solving (2.1) when our prior model assumes $f_{\text{true}}$ to be compressible with respect to a given dictionary $\mathbb{D} :=\{\unicode[STIX]{x1D719}_{i}\}$ . Then sparse recovery of $f$ can be done by either the synthesis approach (i.e. sparse coding) or the analysis approach. In the synthesis approach the reconstruction operator ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ is given as

(2.20) $$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D703}}(g) :={\mathcal{E}}_{\mathbb{D}}^{\ast }(\hat{\unicode[STIX]{x1D709}}),\end{eqnarray}$$

where

with $\unicode[STIX]{x1D703}=\{\mathbb{D},\unicode[STIX]{x1D706}\}$ , i.e.  $\unicode[STIX]{x1D703}$ is the scalar $\unicode[STIX]{x1D706}>0$ and the entire dictionary that defines the synthesis operator ${\mathcal{E}}_{\mathbb{D}}^{\ast }$ . In the corresponding analysis approach, we get

(2.21)

If $\mathbb{D}$ is an orthonormal basis then synthesis and analysis formulations are equivalent.

There are different strategies for numerically computing sparse representations that solve (2.20) and (2.21).

Greedy approaches. These build up an approximation in a greedy fashion, computing one non-zero entry of $\hat{\unicode[STIX]{x1D709}}$ at a time by making locally optimal choices at each step. One example of a greedy approach is iterative (hard) thresholding (Blumensath and Davies Reference Blumensath and Davies2008, Blumensath Reference Blumensath2013, Foucart Reference Foucart2016), where for an initial guess $\boldsymbol{\unicode[STIX]{x1D709}}^{(i0}$ and $i=0,1,\ldots$ one iterates

(2.22) $$\begin{eqnarray}\boldsymbol{\unicode[STIX]{x1D709}}^{(i+1)}=T_{s}(\boldsymbol{\unicode[STIX]{x1D709}}^{(i)}-{\mathcal{A}}_{\mathbb{ D}}^{\ast }({\mathcal{A}}_{\mathbb{ D}}(\boldsymbol{\unicode[STIX]{x1D709}}^{(i)})-g)),\end{eqnarray}$$

where ${\mathcal{A}}_{\mathbb{D}}:={\mathcal{A}}\circ {\mathcal{E}}_{\mathbb{D}}^{\ast }$ . Here, $T_{s}(\boldsymbol{\unicode[STIX]{x1D709}})$ sets all but the largest (in magnitude) $s$ elements of $\boldsymbol{\unicode[STIX]{x1D709}}$ to zero. This is therefore a proximal-gradient method with the proximal of the function being $0$ at $0$ and $1$ everywhere else (see Section A.1). Other examples are matching pursuit (MP) (Mallat and Zhang Reference Mallat and Zhang1993), orthogonal matching pursuit (OMP) (Tropp and Gilbert Reference Tropp and Gilbert2007) and variants thereof such as StOMP (Donoho, Tsaig, Drori and Starck Reference Donoho, Tsaig, Drori and Starck2012), ROMP (Needell and Vershynin Reference Needell and Vershynin2009) and CoSamp (Needell and Tropp Reference Needell and Tropp2009).

Convex relaxation. One of the most common approaches to solving sparse recovery is to replace the $\ell _{0}$ -(semi)norm with the $\ell _{1}$ -norm in (2.20) and (2.21). This leads to basis pursuit (Candès, Romberg and Tao Reference Candès, Romberg and Tao2006), also called Lasso in the statistics literature (Tibshirani Reference Tibshirani1996). The optimization literature for solving the resulting $\ell _{1}$ -type problems is vast. A few examples are interior-point methods (Candès et al. Reference Candès, Romberg and Tao2006, Kim et al. Reference Kim, Koh, Lustig, Boyd and Gorinevsky2007), projected gradient methods (Figueiredo, Nowak and Wright Reference Figueiredo, Nowak and Wright2007), iterative soft thresholding (see Section A.3) (Daubechies et al. Reference Daubechies, Defrise and De Mol2004, Fornasier and Rauhut Reference Fornasier and Rauhut2008) and fast proximal gradient methods (FISTA and variants) (Bubeck Reference Bubeck2015), to name just a few.

Combinatorial algorithms. A third class of approaches for sparse coding is that of combinatorial algorithms. They are particularly suitable when acquiring highly structured samples of the signal so that rapid reconstruction via group testing is efficient. This class of approaches includes Fourier sampling, chaining pursuit and HHS pursuit (Berinde et al. Reference Berinde, Gilbert, Indyk, Karloff and Strauss2008).

The above approaches for solving (2.20) and (2.21) all have their advantages and disadvantages. First of all, greedy methods will generally not give the same solution as convex relaxation. However, if the restricted isometry property (RIP) from Section 2.7.3 holds, then both approaches have the same solution. Convex relaxation has the advantage that it succeeds with a very small number of possibly noisy measurements. However, their numerical solution tends to be computationally burdensome. Combinatorial algorithms, on the other hand, can be extremely fast (sublinear in the length of the target signal) but they require a very specific structure of the forward operator ${\mathcal{A}}$ and a large number of samples. The performance of greedy methods falls in between those of convex relaxation and combinatorial algorithms in their run-time and sampling efficiency.

2.7.3 Error estimates for sparse recovery

When the true model parameter is compressible, then it is possible to estimate the error committed by performing sparse recovery. Such estimates have been derived in the finite-dimensional setting when the matrix representing the linear forward operator satisfies the RIP:

(2.23) $$\begin{eqnarray}(1-\unicode[STIX]{x1D716}_{s})\Vert f\Vert _{2}^{2}\leq \Vert {\mathcal{A}}(f)\Vert _{2}^{2}\leq (1+\unicode[STIX]{x1D716}_{s})\Vert f\Vert _{2}^{2}\quad \text{for all }s\text{-sparse }f\in X,\end{eqnarray}$$

for sufficiently small $\unicode[STIX]{x1D716}_{s}>0$ . Then we have the following error estimate for sparse recovery.

Theorem 2.8 (Candès et al. Reference Candès, Romberg and Tao2006).

Let ${\mathcal{A}}:\mathbb{R}^{n}\rightarrow \mathbb{R}^{m}$ be a linear mapping whose matrix satisfies the RIP. If $g={\mathcal{A}}(f_{\text{true}})+e$ with $\Vert e\Vert \leq \unicode[STIX]{x1D6FF}$ and

(2.24)

then

(2.25) $$\begin{eqnarray}\Vert \hat{f}_{\unicode[STIX]{x1D6FF}}-f_{\text{true}}\Vert _{2}\leq C\biggl[\unicode[STIX]{x1D6FF}+\frac{\Vert f_{\text{true}}-f_{\text{true}}^{(s)}\Vert _{2}}{\sqrt{s}}\biggr].\end{eqnarray}$$

In the above, $f_{\text{true}}^{(s)}\in \mathbb{R}^{n}$ is a vector consisting of the $s$ largest (in magnitude) coefficients of $f_{\text{true}}$ and zeros otherwise.

Examples of matrices satisfying RIP are sub-Gaussian matrices and partial bounded orthogonal matrices (Chen and Needell Reference Chen and Needell2016). Theorem 2.8 states that the reconstruction error is at most proportional to the norm of the noise in the data plus the tail $f_{\text{true}}-f_{\text{true}}^{(s)}$ of the signal. Cohen, Dahmen and DeVore (Reference Cohen, Dahmen and DeVore2009) show that this error bound is optimal (up to the precise value of $C$ ). Moreover, if $f_{\text{true}}$ is $s$ -sparse and $\unicode[STIX]{x1D6FF}=0$ (noise-free data), then $f_{\text{true}}$ can be reconstructed exactly. Furthermore, if $f_{\text{true}}$ is compressible with (2.19), then

(2.26) $$\begin{eqnarray}\Vert \hat{f}_{\unicode[STIX]{x1D6FF}}-f_{\text{true}}\Vert _{2}\leq C(\unicode[STIX]{x1D6FF}+C^{\prime }s^{1/2-1/q}).\end{eqnarray}$$

Finally, error estimates of the above type have been extended to the infinite-dimensional setting in Adcock and Hansen (Reference Adcock and Hansen2016).

The choice of dictionary is clearly a central topic in sparsity-promoting regularization and, as outlined in Section 4.4, the dictionary can be learned beforehand or jointly alongside the signal recovery.

2.7.4 Error estimates for convex relaxation

Applying the convex relaxation to the analysis approach of (2.21) yields the regularized Tikhonov functional

(2.27)

with $1\leq p<2$ . In the context of ill-posed inverse problems this functional was introduced and analysed in the ground-breaking paper by Daubechies et al. (Reference Daubechies, Defrise and De Mol2004). We emphasize that this analysis holds in infinite-dimensional function spaces, and does not depend on the finite-dimensional concepts used in compressive sampling or finite-dimensional sparse recovery concepts.

Since then, this model has been studied intensively, and in particular the case $p=1$ . Similar to Theorem 2.7, error estimates and regularizing properties, such as existence of minimizers, well-posedness, stability, convergence rates and error estimates, have been obtained for linear and non-linear operators: see e.g. Scherzer et al. (Reference Scherzer, Grasmair, Grossauer, Haltmeier and Lenzen2009). Minimizers of this functional with $p=1$ are indeed sparse even if the true solution is not. Using the notion of Bregman distances, such sparsity-promoting approaches have been extended to rather general choices of data discrepancy and regularizers. For a more complete introduction, see Chan and Shen (Reference Chan and Shen2006), Scherzer et al. (Reference Scherzer, Grasmair, Grossauer, Haltmeier and Lenzen2009) and Bredies et al. (Reference Bredies, Kunisch and Pock2011) and the recent survey by Benning and Burger (Reference Benning and Burger2018, Section 2).

3.1.1 The data model

A key assumption is that the conditional distribution of $(\unicode[STIX]{x1D558}\mid \unicode[STIX]{x1D557}=f)\sim \unicode[STIX]{x1D6F1}_{\text{data}}^{f}$ exists (data likelihood) under the joint law $\unicode[STIX]{x1D707}$ for any $f\in X$ . This allows us to define the data model, which is the statistical analogue of the forward operator, as the following $\mathscr{P}_{Y}$ -valued mapping defined on $X$ :

(3.1) $$\begin{eqnarray}f\mapsto \unicode[STIX]{x1D6F1}_{\text{data}}^{f}\quad \text{for any }f\in X.\end{eqnarray}$$

The most common data model is the statistical analogue of (2.1) where the model parameter is allowed to be a random variable (Kallenberg Reference Kallenberg2002, Lemma 1.28 and Corollary 3.12):

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D558}={\mathcal{A}}(\unicode[STIX]{x1D557})+\unicode[STIX]{x1D556}.\end{eqnarray}$$

Here, ${\mathcal{A}}:X\rightarrow Y$ is the same forward operator as in (2.1), which models how data are generated in the absence of noise. Likewise, $\unicode[STIX]{x1D556}\sim \unicode[STIX]{x1D6F1}_{\text{noise}}$ , with $\unicode[STIX]{x1D6F1}_{\text{noise}}\in \mathscr{P}_{Y}$ known, is the random variable that generates the observation noise. If $\unicode[STIX]{x1D556}$ is independent from $\unicode[STIX]{x1D557}$ , then (3.2) amounts to the data model

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\text{data}}^{f}=\unicode[STIX]{x1D6FF}_{{\mathcal{A}}(f)}\circledast \unicode[STIX]{x1D6F1}_{\text{noise}}=\unicode[STIX]{x1D6F1}_{\text{noise}}(\,\cdot \,-{\mathcal{A}}(f))\quad \text{for any }f\in X,\end{eqnarray}$$

where $\circledast$ denotes convolution between measures.

3.1.2 The inverse problem

Following Evans and Stark (Reference Evans and Stark2002), a (statistical) inverse problem requires us to perform the task of recovering the posterior given a single sample (measured data) from the data model with unknown true model parameter. A more precise statement reads as follows.

Definition 3.5. A statistical inverse problem is the task of recovering the conditional distribution $\unicode[STIX]{x1D6F1}_{\text{post}}^{g}\in \mathscr{P}_{Y}$ of $(\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g)$ under $\unicode[STIX]{x1D707}$ from measured data $g\in Y$ , where

(3.4) $$\begin{eqnarray}g\text{ is a single sample of }(\unicode[STIX]{x1D558}\mid \unicode[STIX]{x1D557}=f_{\text{true}})\sim \unicode[STIX]{x1D6F1}_{\text{data}}^{f_{\text{true}}}.\end{eqnarray}$$

Here $f_{\text{true}}\in X$ is unknown while $f\mapsto \unicode[STIX]{x1D6F1}_{\text{data}}^{f}$ , which describes how data are generated, is known.

The conceptual difference that comes from adopting such a statistical view brings with it several potential advantages. The posterior, assuming it exists, describes all possible solutions, so recovering it represents a more complete solution to the inverse problem than recovering an approximation of $f_{\text{true}}$ , which is the goal in functional analytic regularization. This is particularly the case when one seeks to quantify the uncertainty in the recovered model parameter in terms of statistical properties of the data. However, recovering the entire posterior is often not feasible, such as in inverse problems that arise in imaging. As an alternative, one can settle for exploring the posterior by computing suitable estimators (Section 3.3). Some may serve as approximations of $f_{\text{true}}$ whereas others are designed for quantifying the uncertainty.

3.1.3 Bayes’ theorem

A key part of solving the above inverse problem is to utilize a relation between the unknown posterior that one seeks to recover and the known data likelihood. Such a relation is given by Bayes’ theorem.

A very general formulation of Bayes’ theorem is given by Schervish (Reference Schervish1995, Theorem 1.31). For simplicity we consider the formulation that holds for the special case where the data likelihood is given as in (3.2), with ${\mathcal{A}}:X\rightarrow Y$ a measurable map. Assume furthermore that $\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D6F1}_{\text{prior}}$ and $\unicode[STIX]{x1D556}\sim \unicode[STIX]{x1D6F1}_{\text{noise}}$ , with $\unicode[STIX]{x1D556}$ independent of $\unicode[STIX]{x1D557}$ . Then the data model is given as in (3.3), that is, at $f$ it yields the translate of $\unicode[STIX]{x1D6F1}_{\text{noise}}$ by ${\mathcal{A}}(f)$ .

Assume next that $\unicode[STIX]{x1D6F1}_{\text{data}}^{f}\ll \unicode[STIX]{x1D6F1}_{\text{noise}}$ (i.e.  $\unicode[STIX]{x1D6F1}_{\text{data}}^{f}$ is absolutely continuous with respect to $\unicode[STIX]{x1D6F1}_{\text{noise}}$ ) $\unicode[STIX]{x1D6F1}_{\text{noise}}$ -almost surely for all $f\in X$ , so there exists some measurable map $\unicode[STIX]{x1D6F7}:X\times Y\rightarrow \mathbb{R}$ (potential) such that

(3.5) $$\begin{eqnarray}\frac{\,\text{d}\unicode[STIX]{x1D6F1}_{\text{data}}^{f}}{\,\text{d}\unicode[STIX]{x1D6F1}_{\text{noise}}}(g)=\exp (-{\mathcal{L}}(f,g))\quad \text{for all }f\in X.\end{eqnarray}$$

with

$$\begin{eqnarray}\mathbb{E}_{\unicode[STIX]{x1D558}\sim \unicode[STIX]{x1D6F1}_{\text{noise}}}[\exp (-{\mathcal{L}}(f,\unicode[STIX]{x1D558}))]=1.\end{eqnarray}$$

The mapping $f\mapsto -{\mathcal{L}}(f,g)$ is called the (data) log-likelihood for the data $g\in Y$ .

Remark 3.6. The equality for the Radon–Nikodym derivative in (3.5) means that

$$\begin{eqnarray}\mathbb{E}_{\unicode[STIX]{x1D558}\sim \unicode[STIX]{x1D6F1}_{\text{data}}^{f}}[F(\unicode[STIX]{x1D558})]=\mathbb{E}_{\unicode[STIX]{x1D558}\sim \unicode[STIX]{x1D6F1}_{\text{noise}}}[\exp (-{\mathcal{L}}(f,\unicode[STIX]{x1D558}))F(\unicode[STIX]{x1D558})]\end{eqnarray}$$

holds for any measurable $F:Y\rightarrow \mathbb{R}$ .

Finally, assume ${\mathcal{L}}$ is $\unicode[STIX]{x1D707}_{0}$ -measurable, where $\unicode[STIX]{x1D707}_{0}:=\unicode[STIX]{x1D6F1}_{\text{prior}}\otimes \unicode[STIX]{x1D6F1}_{\text{noise}}$ and $\unicode[STIX]{x1D707}\ll \unicode[STIX]{x1D707}_{0}$ , which means in particular that the joint law $(\unicode[STIX]{x1D557},\unicode[STIX]{x1D558})\sim \unicode[STIX]{x1D707}$ can be written as

$$\begin{eqnarray}\frac{\,\text{d}\unicode[STIX]{x1D707}}{\,\text{d}\unicode[STIX]{x1D707}_{0}}(f,g)=\exp (-{\mathcal{L}}(f,g))\quad \text{for }(f,g)\in X\times Y.\end{eqnarray}$$

Bearing in mind the above assumptions, we can now state Bayes’ theorem (Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Theorem 14), which expresses the posterior in terms of the data likelihood and the prior.

Theorem 3.7 (Bayes’ theorem).

The normalization constant $Z:Y\rightarrow \mathbb{R}$ is defined by

(3.6) $$\begin{eqnarray}Z(g) :=\mathbb{E}_{\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D6F1}_{\text{prior}}}[\exp (-{\mathcal{L}}(\unicode[STIX]{x1D557},g))]\end{eqnarray}$$

and we assume $Z(g)>0$ holds $\unicode[STIX]{x1D6F1}_{\text{noise}}$ -almost surely for $g\in Y$ . Then the posterior $\unicode[STIX]{x1D6F1}_{\text{post}}^{g}$ , which is the conditional distribution of $(\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g)$ , exists under $\unicode[STIX]{x1D707}$ and it is absolutely continuous with respect to the prior, i.e.  $\unicode[STIX]{x1D6F1}_{\text{post}}^{g}\ll \unicode[STIX]{x1D6F1}_{\text{prior}}$ . Furthermore,

(3.7) $$\begin{eqnarray}\frac{\,\text{d}\unicode[STIX]{x1D6F1}_{\text{post}}^{g}}{\,\text{d}\unicode[STIX]{x1D6F1}_{\text{prior}}}(f)=\frac{1}{Z(g)}\exp (-{\mathcal{L}}(f,g))\end{eqnarray}$$

holds $\unicode[STIX]{x1D707}_{0}$ -almost surely for $(f,g)\in X\times Y$ .

Bayes’ theorem is the basis for Bayesian inversion, where one seeks to solve the statistical inverse problem assuming access to both a prior and a data likelihood, whereas $f_{\text{true}}\in X$ remains unknown. The data likelihood is given by the data model, which is in turn derived from knowledge about how data are generated. The choice of prior, however, is more subtle: it needs to act as a regularizer, and ideally it also encodes subjective prior beliefs about the unknown model parameter $f_{\text{true}}$ by giving high probability to model parameters similar to $f_{\text{true}}$ and low probability to other ‘unnatural’ model parameters. A brief survey of hand-crafted priors is provided in Section 3.4.

3.2.1 Existence

Existence for Bayesian inversion follows when Bayes’ theorem holds. Below we state the precise existence theorem (Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Theorem 16) for the setting in Section 3.1.3, which covers the case when model parameter and data spaces are infinite-dimensional separable Banach spaces.

Theorem 3.9 (existence for Bayes inversion).

Assume that ${\mathcal{L}}:X\times Y\rightarrow \mathbb{R}$ in (3.5) is continuously differentiable on some $X^{\prime }\subset X$ that contains the support of the prior $\unicode[STIX]{x1D6F1}_{\text{prior}}$ and $\unicode[STIX]{x1D6F1}_{\text{prior}}(X^{\prime }\bigcap B)>0$ for some bounded set $B\subset X$ . Also, assume there exists mappings $M_{1},M_{2}:\mathbb{R}_{+}\times \mathbb{R}_{+}\rightarrow \mathbb{R}_{+}$ that are component-wise monotone, non-decreasing and where the following holds:

(3.8) $$\begin{eqnarray}\begin{array}{@{}c@{}}-{\mathcal{L}}(f,g)\leq M_{1}(r,\Vert f\Vert ),\\ |{\mathcal{L}}(f,g)-{\mathcal{L}}(f,v)|\leq M_{2}(r,\Vert f\Vert )\Vert g-v\Vert \end{array}\end{eqnarray}$$

for $f\in X$ and $g,v\in B_{r}(0)\subset Y$ . Then, $Z$ in (3.6) is finite, i.e.  $0<Z(g)<\infty$ for any $g\in Y$ , and the posterior given by (3.7) yields a well-defined $\mathscr{P}_{X}$ -valued mapping on $Y$ : $g\mapsto \unicode[STIX]{x1D6F1}_{\text{post}}^{g}$ .

Under certain circumstances it is possible to work with improper priors on $X$ , for example by computing posterior distributions that approximate the posteriors one would have obtained using proper conjugate priors whose extreme values coincide with the improper prior.

3.2.2 Stability

One can show that small changes in the data lead to small changes in the posterior distribution (in Hellinger metric) on $\mathscr{P}_{X}$ . The precise formulation given by Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Theorem 16) reads as follows.

Theorem 3.10 (stability for Bayes inversion).

Let the assumptions in Theorem 3.9 and assume in addition that

(3.9) $$\begin{eqnarray}f\mapsto \exp (M_{1}(r,\Vert f\Vert ))(1+M_{2}(r,\Vert f\Vert ))^{2}\quad \text{is }\unicode[STIX]{x1D6F1}_{\text{prior}}\text{-integrable on }X\end{eqnarray}$$

for some fixed $r>0$ . Then

$$\begin{eqnarray}d_{\text{H}}(\unicode[STIX]{x1D6F1}_{\text{post}}^{g},\unicode[STIX]{x1D6F1}_{\text{post}}^{v})\leq C(r)\Vert g+v\Vert\end{eqnarray}$$

for some $C(r)>0$ and $g,v\in B_{r}(0)\subset Y$ . In the above, $d_{\text{H}}:\mathscr{P}_{X}\times \mathscr{P}_{X}\rightarrow \mathbb{R}$ is the Hellinger metric (see Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Definition 4).

Theorem 3.10 holds in particular when the negative log-likelihood of data is locally Hölder-continuous, which is the case for many standard probability distribution functions, for example when the negative log-likelihood is continuously differentiable.

By Theorem 3.10 the posterior is Lipschitz in the data with respect to the Hellinger metric, so rephrasing the inverse problem as the task of recovering the posterior instead of a model parameter acts as a regularization of an inverse problem that is ill-posed in the functional analytic sense. Note, however, that the above does not automatically imply that a particular estimator is continuous with respect to data. However, the Hellinger distance possesses the convenient property that continuity with respect to this metric implies the continuity of moments. Hence, a corollary to Theorem 3.10 is that posterior moments, such as the mean and covariance, are continuous (Sprungk Reference Sprungk2017, Corollary 3.23), that is, these estimators are regularizing.

As a final note, one can also show that small changes in the data log-likelihood ${\mathcal{L}}$ in (3.5) lead to small changes in the posterior distribution, again in the Hellinger metric (Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Theorem 18). This enables one to translate errors arising from inaccurate forward operator into errors in the Bayesian solution of the inverse problem, a topic that is also considered in Section 6.2.

3.2.3 Convergence

Posterior consistency is the Bayesian analogue of the notion of convergence in functional analytic regularization. More precisely, the requirement is that the posterior $\unicode[STIX]{x1D6F1}_{\text{post}}^{g}$ , where $g$ is a sample of $(\unicode[STIX]{x1D558}\mid \unicode[STIX]{x1D557}=f_{\text{true}})$ , concentrates in any small neighbourhood of the true model parameter $f_{\text{true}}\in X$ as information in data $g$ increases indefinitely.Footnote ¹ Intuitively, this means our knowledge about the model parameter becomes more accurate and precise as the amount of data increases indefinitely.

In the finite-dimensional setting, consistency of the posterior holds if and only if $f_{\text{true}}\in X$ is contained in the support of the prior $\unicode[STIX]{x1D6F1}_{\text{prior}}$ (Freedman Reference Freedman1963, Schwartz Reference Schwartz1965) (provided the posterior is smooth with respect to the model parameter).

The situation is vastly more complex in the infinite-dimensional (non-parametric) setting. It is known that the posterior is consistent at every model parameter except possibly on a set of measure zero (Doob Reference Doob1948, Breiman, Le Cam and Schwartz Reference Breiman, Le Cam and Schwartz1965), that is, Bayesian inversion is almost always consistent in the measure-theoretic sense. However, the situation changes if ‘smallness’ is measured in a topological sense, as shown by a classical counter-example in Freedman (Reference Freedman1963) involving the simplest non-parametric problem where consistency fails. This is not a pathological counter-example: it is a generic property in the sense that most priors are ‘bad’ in a topological sense (Freedman Reference Freedman1965). In fact, consistency may fail for non-parametric models for very natural priors satisfying the support condition, which means even an infinite amount of data may not be sufficient to ‘correct’ for errors introduced by a prior. Hence, unlike the finite-dimensional setting, many priors do not ‘wash out’ as the information in the data increases indefinitely, so the prior may have a large influence on the corresponding posterior even in the asymptotic setting. Examples of posterior consistency results for Bayesian inversion are given in Ghosal, Ghosh and Ramamoorthi (Reference Ghosal, Ghosh and Ramamoorthi1999), Ghosal, Ghosh and van der Vaart (Reference Ghosal, Ghosh and van der Vaart2000), Neubauer and Pikkarainen (Reference Neubauer and Pikkarainen2008), Bochkina (Reference Bochkina2013), Agapiou, Larsson and Stuart (Reference Agapiou, Larsson and Stuart2013), Stuart and Teckentrup (Reference Stuart and Teckentrup2018), Kekkonen, Lassas and Siltanen (Reference Kekkonen, Lassas and Siltanen2016) and Kleijn and Zhao (Reference Kleijn and Zhao2018).

To summarize, there are many reasonable priors for which posterior consistency holds at every point of the model parameter space $X$ . A general class of such ‘good’ priors is that of tail-free priors (Freedman Reference Freedman1963) and neutral-to-the-right priors (Doksum Reference Doksum1974); see also Le Cam (Reference Le Cam1986, Section 17.7). The fact that there are too many ‘bad’ priors may therefore not be such a serious concern since there are also enough good priors approximating any given subjective belief. It is therefore important to have general results providing sufficient conditions for the consistency given a pair of model parameter and prior as developed, and Schwartz (Reference Schwartz1965) is an example in this direction.

3.2.4 Convergence rates

Posterior consistency is a weak property shared by many different choices of priors. More insight into the performance of Bayesian inversion under the choice of different priors requires characterizing other properties of the posterior, such as quantifying how quickly it converges to the true solution as the observational noise goes to 0. This leads to results about contraction rates, which is the Bayesian analogue of convergence rate theorems in functional analytic regularization.

Formally, consider the setting in (3.2) where the observational noise tends to zero as some scalar $\unicode[STIX]{x1D6FF}\rightarrow 0$ , that is, for some fixed $f_{\text{true}}\in X$ we have

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D558}={\mathcal{A}}(f_{\text{true}})+\unicode[STIX]{x1D556}_{\unicode[STIX]{x1D6FF}}\quad \text{where }\Vert \unicode[STIX]{x1D556}_{\unicode[STIX]{x1D6FF}}\Vert \rightarrow 0\text{ as }\unicode[STIX]{x1D6FF}\rightarrow 0.\end{eqnarray}$$

A contraction rate theorem seeks to find the base rate for $\unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF})\rightarrow 0$ such that

$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{\text{post}}^{g^{\unicode[STIX]{x1D6FF}}}(\{f\in X:\ell _{X}(f,f_{\text{true}})\geq \unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF})\})\rightarrow 0\end{eqnarray}$$

$\unicode[STIX]{x1D6F1}_{\text{post}}^{g_{0}}$ -almost surely as $\unicode[STIX]{x1D716}\rightarrow 0$ . Here, $g^{\unicode[STIX]{x1D6FF}}\in Y$ is a single sample of $\unicode[STIX]{x1D558}$ in (3.10), $g_{0}={\mathcal{A}}(f_{\text{true}})$ is the corresponding ideal data, $\{\unicode[STIX]{x1D6F1}_{\text{prior}}^{\unicode[STIX]{x1D6FF}}\}_{\unicode[STIX]{x1D6FF}}\subset \mathscr{P}_{X}$ is a sequence of prior distributions with $f\sim \unicode[STIX]{x1D6F1}_{\text{prior}}^{\unicode[STIX]{x1D6FF}}$ , and $\ell _{X}:X\times X\rightarrow \mathbb{R}$ is a (measurable) distance function on $X$ . Research in this area has mainly focused on (a) obtaining contraction rates for Bayesian inversion where the prior is from some large class, or (b) improving the rates by changing the parameters of the prior depending on the level of noise (and even the data). Most of the work is done for additive observational white Gaussian noise, that is, the setting in (3.10) with $\unicode[STIX]{x1D556}_{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}W$ , where $W$ is a centred Gaussian white noise process that can be defined by its action on a separable Hilbert space $Y$ . See Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Sections A.3 and A.4) for a survey related to Gaussian measures and Wiener processes on separable Banach spaces.

As one might expect, initial results were for linear forward operators. The first results restricted attention to conjugate priors where the posterior has an explicit expression (it is in the same family as the prior) (Liese and Miescke Reference Liese and Miescke2008, Section 1.2). However, Bayesian non-parametric statistics often carries over to the inverse setting via the singular value decomposition of the linear forward operator. This allowed one to prove contraction rate results for Bayesian inversion with non-conjugate priors (Knapik, van der Vaart and van Zanten Reference Knapik, van der Vaart and van Zanten2011, Knapik, van der Vaart and van Zanten Reference Knapik, van der Vaart and van Zanten2013, Ray Reference Ray2013, Agapiou et al. Reference Agapiou, Larsson and Stuart2013, Agapiou, Stuart and Zhang Reference Agapiou, Stuart and Zhang2014). In particular, if the Gaussian prior is non-diagonal for the singular value decomposition (SVD), then the posterior is still Gaussian and its contraction rate will be driven by the convergence rate of its posterior mean (since the variance does not depend on the data). Furthermore, in the linear setting, the posterior mean is the MAP, so the convergence rate of the MAP will be a posterior contraction rate.

Unfortunately, not many of the methods developed for proving contraction rates in the linear setting carry over to proving contraction rates for general, non-linear, inverse problems. A specific difficulty in the Bayesian setting is that the noise term often does not take values in the natural range of the forward operator. For example, consider the data model in (3.3) with observational white Gaussian noise, i.e. (3.10) with $\unicode[STIX]{x1D556}_{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}W$ , where $W$ is a centred Gaussian white noise process that can be defined by its action on a separable Hilbert space $Y$ . If we have a non-linear forward operator ${\mathcal{A}}$ given as the solution operator to an elliptic PDE, then the noise process $W$ does not define a proper random element in $Y=L^{2}(\mathbb{M})$ for $\mathbb{M}\subset \mathbb{R}^{d}$ : instead it defines a random variable only in a negative Sobolev space $W^{-\unicode[STIX]{x1D6FD}}$ with $\unicode[STIX]{x1D6FD}>d/2$ . Nevertheless, there are also some results in the non-linear setting: for example, Kekkonen et al. (Reference Kekkonen, Lassas and Siltanen2016) derive contraction rate results for Bayesian inversion of inverse problems under Gaussian conjugate priors and where the forward operator is a linear hypoelliptic pseudodifferential operator. Recent developments that make use of techniques from concentration of measure theory have resulted in contraction rate theorems outside the conjugate setting. For example, Nickl and Söhl (Reference Nickl and Söhl2017) and Nickl (Reference Nickl2017a ) derive contraction rates for Bayesian inversion for parameter estimation involving certain classes of elliptic PDEs that are minimax-optimal in prediction loss: see e.g. Nickl (Reference Nickl2017a , Theorem 28) for an example of a general contraction theorem. The case of a (possibly) non-linear forward operator and a Gaussian prior is considered in Nickl, van de Geer and Wang (Reference Nickl, van de Geer and Wang2018), which studies properties of the MAP estimator. Finally, Gugushvili, van der Vaart and Yan (Reference Gugushvili, van der Vaart and Yan2018) derive contraction rates for ‘general’ priors expressed by scales of smoothness classes. The precise conditions are checked only for an elliptic example: it is not clear whether it works in other examples, such as the ray transform.

3.2.5 Characterization of the posterior for uncertainty quantification

Posterior consistency and contraction rates are relevant results, but assessing the performance of Bayesian inversion methods in uncertainty quantification requires a more precise characterization of the posterior. The aim is to characterize the fluctuations of $(\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g)$ near $f_{\text{true}}$ when scaled by some inverse contraction rate.

One approach is to derive Bernstein–von Mises type theorems, which characterize the posterior distribution in terms of a canonical Gaussian distribution in the small-noise or large-sample limit. To better illustrate the role of such theorems, we consider the special case with Gaussian prior on the Hilbert space $X$ and observational noise $\unicode[STIX]{x1D556}_{\unicode[STIX]{x1D6FF}}=\unicode[STIX]{x1D6FF}\mathbb{W}$ in (3.10) (white noise model) where $\mathbb{W}$ denotes a Gaussian white noise process in the Hilbert space $Y$ , that is,

$$\begin{eqnarray}\unicode[STIX]{x1D558}^{\unicode[STIX]{x1D6FF}}={\mathcal{A}}(f_{\text{true}})+\unicode[STIX]{x1D6FF}\mathbb{W}.\end{eqnarray}$$

Many of the results on contraction rates (Section 3.2.4) for Bayesian inversion in this setting are obtained for the distance function induced by the $L^{2}$ -norm, so it is natural to initially consider the statistical fluctuations in $L^{2}$ of the random variable

$$\begin{eqnarray}\unicode[STIX]{x1D56B}_{\unicode[STIX]{x1D6FF}} :=\frac{1}{\unicode[STIX]{x1D716}_{\unicode[STIX]{x1D6FF}}}(\mathbb{E}[\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g^{\unicode[STIX]{x1D6FF}}]-f_{\text{true}}).\end{eqnarray}$$

Now, it turns out that there is no Bernstein–von Mises type asymptotics for $\unicode[STIX]{x1D56B}_{\unicode[STIX]{x1D6FF}}$ as the noise level $\unicode[STIX]{x1D6FF}$ tends to zero, that is, there is no Gaussian process $(\mathbb{G}(\unicode[STIX]{x1D719}))_{\unicode[STIX]{x1D719}\in C^{\infty }}$ such that

(3.11) $$\begin{eqnarray}\biggl(\frac{1}{\unicode[STIX]{x1D6FF}}\langle \unicode[STIX]{x1D56B}_{\unicode[STIX]{x1D6FF}}-\mathbb{E}[\unicode[STIX]{x1D56B}_{\unicode[STIX]{x1D6FF}}],\unicode[STIX]{x1D719}\rangle _{L^{2}}\biggr)_{\unicode[STIX]{x1D719}\in C^{\infty }}\rightarrow (\mathbb{G}(\unicode[STIX]{x1D719}))_{\unicode[STIX]{x1D719}\in C^{\infty }}\quad \text{weakly as }\unicode[STIX]{x1D6FF}\rightarrow 0.\end{eqnarray}$$

To sidestep this difficulty, Castillo and Nickl (Reference Castillo and Nickl2013, Reference Castillo and Nickl2014) seek to determine maximal families $\unicode[STIX]{x1D6F9}$ that replace $C^{\infty }$ in (3.11) and where such an asymptotic characterization holds. This leads to non-parametric Bernstein–von Mises theorems, and while Castillo and Nickl (Reference Castillo and Nickl2013, Reference Castillo and Nickl2014) considered ‘direct’ problems in non-parametric regression and probability density estimation, recent papers have obtained non-parametric Bernstein–von Mises theorems for certain classes of inverse problems. For example, Monard, Nickl and Paternain (Reference Monard, Nickl and Paternain2019) consider the case of inverting the (generalized) ray transform, whereas Nickl (Reference Nickl2017a ) considers PDE parameter estimation problems. The case with general linear forward problem is treated by Giordano and Kekkonen (Reference Giordano and Kekkonen2018), who build upon the techniques of Monard et al. (Reference Monard, Nickl and Paternain2019).

To give a flavour of the type of results obtained, we consider Theorem 2.5 of Monard et al. (Reference Monard, Nickl and Paternain2019), which is relevant to tomographic inverse problems involving inversion of the ray transform for recovering a function (images) defined on $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ , i.e.  $X\subset L^{2}(\unicode[STIX]{x1D6FA})$ . This theorem states that

(3.12) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D6FF}}\langle (\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g^{\unicode[STIX]{x1D6FF}})-\mathbb{E}[\unicode[STIX]{x1D557}|\unicode[STIX]{x1D558}=g^{\unicode[STIX]{x1D6FF}}],\unicode[STIX]{x1D719}\rangle _{L^{2}}\rightarrow \operatorname{N}(0,\Vert {\mathcal{A}}\circ ({\mathcal{A}}^{\ast }\circ {\mathcal{A}})^{-1}(\unicode[STIX]{x1D719})\Vert _{Y})\end{eqnarray}$$

as $\unicode[STIX]{x1D6FF}\rightarrow 0$ for any $\unicode[STIX]{x1D719}\in C^{\infty }(\unicode[STIX]{x1D6FA})$ .

The convergence is in $\unicode[STIX]{x1D6F1}_{\text{prior}}$ -probability and the $Y$ -norm on the right-hand side is a natural $L^{2}$ -norm on the range of the ray transform. The limiting covariance is also shown to be minimal, that is, it attains the semi-parametric Cramér–Rao lower bound (or ‘inverse Fisher information’) for estimating $\langle f,\unicode[STIX]{x1D719}\rangle _{L^{2}}$ near $f_{\text{true}}$ . A key step in the proof is to show smoothness of the ‘Fisher information’ operator $({\mathcal{A}}^{\ast }\circ {\mathcal{A}})^{-1}:X\rightarrow X$ , which is done using techniques from microlocal analysis (Monard et al. Reference Monard, Nickl and Paternain2019, Theorem 2.2). The existence and mapping properties of this inverse Fisher information operator in (3.12) also plays a crucial role in proving that Bernstein–von Mises theorem for other inverse problems. For those with a non-linear forward operator, the information operator to be inverted is found after linearization as shown in Nickl (Reference Nickl2017a ) for parameter estimation in an elliptic PDE.

3.3.1 Some point estimators

Here we briefly define the most common point estimators that are commonly considered in solving many inverse problems. Overall, estimators that include an integration over $X$ are computationally demanding. These include the conditional mean and many estimators associated with uncertainty quantification.

Maximum likelihood estimator.

This estimator maximizes the data likelihood. Under the data model in (3.5), the (maximum likelihood) reconstruction operator ${\mathcal{R}}:Y\rightarrow X$ becomes

The advantage of the maximum likelihood estimator is that it does not involve any integration over $X$ , so it is computationally feasible to use on large-scale problems. Furthermore, it only requires access to the data likelihood (no need to specify a prior), which is known. On the other hand, it does not act as a regularizer, so it is not suitable for ill-posed problems.

Maximum a posteriori (MAP) estimator.

This estimator maximizes the posterior probability, that is, it is the ‘most likely’ model parameter given measured data $g$ . In the finite-dimensional setting, the prior and posterior distribution can typically be described by densities with respect to the Lebesgue measure, and the MAP estimator is defined as the reconstruction operator ${\mathcal{R}}:Y\rightarrow X$ given by

The second equality above assumes a data model as in (3.5).

For many of the infinite-dimensional spaces there exists no analogue of the Lebesgue measure, which makes it difficult to define a MAP estimator through densities. One way to work around this technical problem is to replace the Lebesgue measure with a Gaussian measure on $X$ . Hence, we assume the posterior and prior have densities with respect to some fixed centred (mean-zero) Gaussian measure $\unicode[STIX]{x1D707}_{0}$ and $E$ denotes its Cameron–Martin space.Footnote ³ Now, following Dashti, Law, Stuart and Voss (Reference Dashti, Law, Stuart and Voss2013), we consider the centre of a small ball in $X$ with maximal probability, and then study the limit of this centre as the radius of the ball shrinks to zero. Stated precisely, given fixed data $g\in Y$ , assume there is a functional ${\mathcal{J}}:E\rightarrow \mathbb{R}$ that satisfies

$$\begin{eqnarray}\lim _{r\rightarrow 0}\frac{\unicode[STIX]{x1D6F1}_{\text{post}}^{g}(B_{r}(f_{1}))}{\unicode[STIX]{x1D6F1}_{\text{post}}^{g}(B_{r}(f_{2}))}=\exp ({\mathcal{J}}(f_{1})-{\mathcal{J}}(f_{2})).\end{eqnarray}$$

Here, $B_{r}(f)\subset X$ is the open ball of radius $r>0$ centred at $f\in X$ and ${\mathcal{J}}$ is the Onsager–Machlup functional (Ikeda and Watanabe Reference Ikeda and Watanabe1989, p. 533).

For any fixed $f_{1}\in X$ , a model parameter $f_{2}\in X$ for which the above limit is maximal is a natural candidate for the MAP estimator and is clearly given by minimizers of the Onsager–Machlup functional. The advantage of the MAP estimator is that it does not involve any integration over $X$ . Furthermore, the prior often acts as a regularizer, so MAP can be useful for solving an ill-posed inverse problem. A disadvantage is that a MAP may not always exist (Section 3.3.2): one needs to provide an explicit prior, and many priors result in a non-smooth optimization problem. Due to the latter, MAP estimation is computationally more challenging than maximum likelihood estimation.

Conditional (posterior) mean.

The reconstruction operator ${\mathcal{R}}:Y\rightarrow X$ is here defined as

$$\begin{eqnarray}{\mathcal{R}}(g):=\mathbb{E}[\unicode[STIX]{x1D557}\mid \unicode[STIX]{x1D558}=g]=\int _{X}f\,\,\text{d}\unicode[STIX]{x1D6F1}_{\text{post}}^{g}(f).\end{eqnarray}$$

This estimator involves integration over $X$ , making it challenging to use on even small- to mid-scale problems. It also assumes access to the posterior, which in turn requires a prior. It does, however, act as a regularizer (Sprungk Reference Sprungk2017, Corollary 3.23) and therefore it is suitable for ill-posed problems.

Bayes estimator.

One starts by specifying a function $\ell _{X}:X\times X\rightarrow \mathbb{R}$ that quantifies proximity in $X$ (it does not have to be a metric). Then, a Bayes estimator minimizes the expected loss with respect to the prior $\unicode[STIX]{x1D6F1}_{\text{prior}}$ , that is, ${\mathcal{R}}:Y\rightarrow X$ is defined as

(3.14)

This estimator is challenging to compute even for small- to mid-scale problems, due to the integration over $X\times Y$ and the minimization over all non-randomized decision rules. It also requires access to the joint distribution $\unicode[STIX]{x1D707}$ , which by the law of total probability is expressible as $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D6F1}_{\text{prior}}\otimes \unicode[STIX]{x1D6F1}_{\text{data}}^{f}$ with known data likelihood $\unicode[STIX]{x1D6F1}_{\text{data}}^{f}$ . An important property of the Bayes estimator is that it is a regularizer, so it is suitable for ill-posed problems. It is also equivalent to the conditional mean when the loss is taken as the square of the $L^{2}$ -norm (Helin and Burger Reference Helin and Burger2015) and the conditional median when the loss is the $L^{1}$ -norm. Furthermore, Theorem 1 of Burger and Lucka (Reference Burger and Lucka2014) shows that a MAP estimator with a Gibbs-type prior where the energy functional is Lipschitz-continuous and convex equals a Bayes estimator where the loss is the Bregman distance of the aforementioned energy functional. Finally, in the finite-dimensional setting, one can show (Banerjee, Guo and Wang Reference Banerjee, Guo and Wang2005) that the Bayes estimator is the same as the conditional mean if and only if $\ell _{X}$ is the Bregman distance of a strictly convex non-negative differentiable functional.Footnote ⁴

3.3.2 Relation to functional analytic regularization

It is quite common to interpret a variational method as in (2.7) as a MAP estimators. This is especially the case when one chooses the data discrepancy ${\mathcal{L}}$ so that minimizing $f\mapsto {\mathcal{L}}({\mathcal{A}}(f),g)$ corresponds to computing a maximum likelihood estimator.

Such an interpretation is almost always possible in the finite-dimensional setting since, as mentioned in Definition 2.2, the regularization functional ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ can be interpreted as a Gibbs-type prior $\unicode[STIX]{x1D70C}_{\text{prior}}(f)\propto \exp (-{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f))$ (Kaipio and Somersalo Reference Kaipio and Somersalo2005) where $\unicode[STIX]{x1D70C}_{\text{prior}}$ is the density for the prior. The situation is more complicated in the infinite-dimensional setting since a MAP estimator does not always exist, and if it exists then there is no general scheme for connecting the topological description of a MAP estimate to a variational problem. The main reason is that the posterior no longer has a natural density representation, which significantly complicates the definition and study of the underlying conditional probabilities.

Assume the prior measure is specified by a Gaussian random field, and the likelihood satisfies conditions in Theorem 3.9 that are necessary for the existence of a well-posed posterior measure. Then the MAP estimator is well-defined as the minimizer of an Onsager–Machlup functional defined on the Cameron–Martin space of the prior. If one has Gaussian noise, then this becomes a least-squares functional penalized by a regularization functional that is the Cameron–Martin norm of the Gaussian prior (Dashti et al. Reference Dashti, Law, Stuart and Voss2013); see also Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Section 4.3). To handle the case with non-Gaussian priors, Helin and Burger (Reference Helin and Burger2015) introduce the notion of weak MAP estimate (wMAP) and show that a wMAP estimate can be connected to a variational formulation also for non-Gaussian priors. Furthermore, any MAP estimate in the sense of Dashti et al. (Reference Dashti, Law, Stuart and Voss2013) is a wMAP estimate, so this is a generalization in the strict sense. The wMAP approach, however, fails when the prior does not admit continuous densities. In order to handle this, Clason, Helin, Kretschmann and Piiroinen (Reference Clason, Helin, Kretschmann and Piiroinen2018) introduce the notion of a generalized mode of a probability measure (Clason et al. Reference Clason, Helin, Kretschmann and Piiroinen2018, Definition 2.3) and define a (generalized) MAP estimator as a generalized mode of the posterior measure. Generalized MAP reduces to the classical MAP in certain situations that include Gaussian priors but the former also exist for a more general class of priors, such as uniform priors. The main result in Clason et al. (Reference Clason, Helin, Kretschmann and Piiroinen2018) is that one can characterize a generalized MAP estimator as a minimizer of an Onsager–Machlup functional even in cases when the prior does not admit a continuous density.

As a final remark, while a MAP estimate with a Gaussian noise model does lead to an optimization problem with a quadratic data-fidelity term, Gribonval and Nikolova (Reference Gribonval and Nikolova2018) show via explicit examples that the converse is not true. They characterize those data models of the type in (3.2) where Bayes estimators can be expressed as the solution of a penalized least-squares optimization problem. One example is denoising in the presence of additive Gaussian noise and an arbitrary prior; another is a data model with (a variant of) Poisson noise and any prior probability on the unknowns. In these cases, the variational approach is rich enough to build all possible Bayes estimators via a well-chosen penalty.

Priors on digitized images.

A wide range of priors for Bayesian inversion have been suggested when the unknown model parameter represents a digitized image, the latter given by a real-valued function in some suitable class defined on a domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{n}$ .

One such class is that of smoothness priors, which reflect the belief that values of the model parameter at a point are close to the average of its values in a neighbourhood of that point. Another is structural priors, which allow for more abrupt (discontinuous) changes in the values of the unknown model parameter at specific locations. Yet another is sparsity-promoting priors (see Section 2.7), which encode the a priori belief that the unknown model parameter is compressible with respect to some underlying dictionary, that is, it can be transformed into a linear combination of dictionary elements where most coefficients vanish (or are small). Finally, there are hierarchical priors which are formed by combining other priors hierarchically into an overall prior. This is typically done in a two-step process where one first specifies some underlying priors, often taken as natural conjugate priors, and then mixes them in a second stage over hyper-parameters. Recently, Calvetti, Somersalo and Strang (Reference Calvetti, Somersalo and Strang2019) have reformulated the question of sparse recovery as an inverse problem in the Bayesian framework, and expressed the sparsity criteria by means of a hierarchical prior mode. More information and further examples of priors in the finite-dimensional setting are given by Kaipio and Somersalo (Reference Kaipio and Somersalo2005), Calvetti and Somersalo (Reference Calvetti and Somersalo2008, Reference Calvetti and Somersalo2017) and Calvetti et al. (Reference Calvetti, Somersalo and Strang2019).

Priors on function spaces.

Defining priors when $X$ is an infinite-dimensional function space is somewhat involved. A common approach is to consider a convergent series expansion and then let the coefficients be generated by a random variable.

More precisely, consider the case when $X$ is a Banach space of real-valued functions on some fixed domain $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{d}$ . Following Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Section 2.1), let $\{{\unicode[STIX]{x1D719}_{i}\}}_{i}\subset X$ be a countable sequence whose elements are normalized, i.e.  $\Vert \unicode[STIX]{x1D719}_{i}\Vert _{X}=1$ . Now consider model parameters $f\in X$ of the form

$$\begin{eqnarray}f=f_{0}+\mathop{\sum }_{i}\unicode[STIX]{x1D6FC}_{i}\unicode[STIX]{x1D719}_{i},\end{eqnarray}$$

where $f_{0}\in X$ is not necessarily normalized to 1. A probability distribution on the coefficients $\unicode[STIX]{x1D6FC}_{i}$ renders a real-valued random function on $\unicode[STIX]{x1D6FA}$ : simply define the deterministic sequence $\{{\unicode[STIX]{x1D6FE}_{i}\}}_{i}\subset \mathbb{R}$ and the i.i.d. random sequence $\{{\unicode[STIX]{x1D709}_{i}\}}_{i}\subset \mathbb{R}$ and set $\unicode[STIX]{x1D6FC}_{i}=\unicode[STIX]{x1D6FE}_{i}\unicode[STIX]{x1D709}_{i}$ . This generates a probability measure $\unicode[STIX]{x1D6F1}_{\text{prior}}$ on $X$ by taking the pushforward of the measure on the i.i.d. random sequence $\{{\unicode[STIX]{x1D709}_{i}\}}_{i}\subset \mathbb{R}$ under the map which takes the sequence into the random function.

Using the above technique, one can construct a uniform prior that can be shown to generate random functions that are all contained in a subset $X\subset L^{\infty }(\unicode[STIX]{x1D6FA})$ , which can be characterized (Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Theorem 2). It is likewise possible to define Gaussian priors where the random function exists as an $L^{2}$ limit in Sobolev spaces of sufficient smoothness (Dashti and Stuart Reference Dashti, Stuart and Ghanem2017, Theorem 8). For imaging applications, much effort has been devoted to constructing edge-preserving priors. An obvious candidate is the TV-prior, but it is not ‘discretization-invariant’ as its edge-preserving property is lost when the discretization becomes finer. This prompted the development of other edge-preserving priors, such as Besov space priors, which are discretization-invariant (Lassas et al. Reference Lassas, Saksman and Siltanen2009, Kolehmainen, Lassas, Niinimäki and Siltanen Reference Kolehmainen, Lassas, Niinimäki and Siltanen2012). See also Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Theorem 5), which characterizes Besov priors that generate random functions contained in separable Banach spaces. However, edge-preserving inversion using a Besov space prior often relies on the Haar wavelet basis. Due to the structure of the Haar basis, discontinuities are preferred on an underlying dyadic grid given by the discontinuities of the basis functions. As an example, on $[0,1]$ a Besov space prior prefers a discontinuity at $x=1/4$ over $x=1/3$ . Thus, in most practical cases, Besov priors rely on both a strong and unrealistic assumption. For this reason, Markkanen, Roininen, Huttunen and Lasanen (Reference Markkanen, Roininen, Huttunen and Lasanen2019) propose another class of priors for edge-preserving Bayesian inversion, the Cauchy difference priors. Starting from continuous one-dimensional Cauchy motion, its discretized version, Cauchy random walk, can be used as a non-Gaussian prior for edge-preserving Bayesian inversion. As shown by Markkanen et al. (Reference Markkanen, Roininen, Huttunen and Lasanen2019), one can also develop a suitable posterior distribution sampling algorithm for computing the conditional mean estimates with single-component Metropolis–Hastings. The approach is applied to CT reconstruction problems in materials science.

The above constructions of random functions through randomized series can be linked to each other through the notion of random fields as shown in Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Section 2.5); see also Ghosal and van der Vaart (Reference Ghosal and van der Vaart2017, Chapters 2 and 10) for further examples. These constructions also extend straightforwardly to $\mathbb{R}^{n}$ - or $\mathbb{C}^{n}$ -valued functions.

3.5.1 Choosing a good prior

Current theory for choosing a ‘good’ prior mainly emphasizes the regularizing function of the prior (Section 3.4). In particular, one seeks a prior that ensures posterior consistency (Section 3.2.3) and good contraction rates (Section 3.2.4) and, if possible, also allows for an asymptotic characterization of the posterior (Section 3.2.5).

This theory, however, is for the asymptotic setting as the noise level in data tends to zero, whereas there seems to be no theory for Bayesian inversion that deals with the setting when the data have a fixed noise level. Here, the prior has a role that goes beyond acting as a regularizer, and its choice may have a large influence on the posterior. For example, it is far from clear whether a prior that provides ‘optimal’ contraction rates is the most suitable one when the data are fixed with a given noise level. Another difficulty is that norms used to quantify distance in the mathematical theorems have little to do with what is meant by a ‘good’ estimate of a model parameter. An archetypal example is the difficulty in quantifying the notion of ‘image quality’ in radiology. This is very difficult since the notion of image quality depends on the task motivating the imaging application. Hand-crafted priors surveyed in Section 3.4 have limitations in this regard, and in Sections 4.3, 4.4, 4.7 and 4.10 we survey work that uses data-driven modelling to obtain a prior.

3.5.2 Computational feasibility

The focus here is on techniques that are not based on deep learning for sampling from a very high-dimensional posterior distribution that lacks an explicit expression.

A well-known class of methods is based on Markov chain Monte Carlo (MCMC), where the aim is to define an iterative process whose stationary distribution coincides with the target distribution, which in Bayesian inversion is the posterior. MCMC techniques come in many variants, and one common variant is MCMC sampling with Metropolis–Hastings dynamics (Minh and Le Minh Reference Minh and Le Minh2015), which generates a Markov chain with equilibrium distribution that coincides with the posterior in the limit. Other variants use Gibbs sampling, which reduces the autocorrelation between samples. Technically, Gibbs sampling can be seen as a special case of Metropolis–Hastings dynamics and it requires computation of conditional distributions. Further variants are auxiliary variable MCMC methods, such as slice sampling (Neal Reference Neal2003), proximal MCMC (Green, Łatuszysński, Pereyra and Robert Reference Green, Łatuszysński, Pereyra and Robert2015, Durmus, Moulines and Pereyra Reference Durmus, Moulines and Pereyra2018, Repetti, Pereyra and Wiaux Reference Repetti, Pereyra and Wiaux2019) and Hamiltonian Monte Carlo (Girolami and Calderhead Reference Girolami and Calderhead2011, Betancourt Reference Betancourt2017). See also Dashti and Stuart (Reference Dashti, Stuart and Ghanem2017, Section 5) for a nice abstract description of MCMC in the context of infinite-dimensional Bayesian inversion.

An alternative approach to MCMC seeks to approximate the posterior with more tractable distributions (deterministic inference), for example in variational Bayes inference (Fox and Roberts Reference Fox and Roberts2012, Blei, Küçükelbir and McAuliffe Reference Blei, Küçükelbir and McAuliffe2017) and expectation propagation (Minka Reference Minka, Breese and Koller2001). Variational Bayes inference has indeed emerged as a popular alternative to the classical MCMC methods for sampling from a difficult-to-compute probability distribution, which in Bayesian inversion is the posterior distribution. The idea is to start from a fixed family of probability distributions (variational family) and select the one that best approximates the target distribution under some similarity measure, such as the Kullback–Leibler divergence.

Blei et al. (Reference Blei, Küçükelbir and McAuliffe2017, p. 860) try to provide some guidance on when to use MCMC and when to use variational Bayes. MCMC methods tend to be more computationally intensive than variational inference, but they also provide guarantees of producing (asymptotically) exact samples from the target density (Robert and Casella Reference Robert and Casella2004). Variational inference does not enjoy such guarantees: it can only find a density close to the target but tends to be faster than MCMC. A recent development is the proof of a Bernstein–von Mises theorem (Wang and Blei Reference Wang and Blei2017, Theorem 5), which shows that the variational Bayes posterior is asymptotically normal around the variational frequentist estimate. Hence, if the variational frequentist estimate is consistent, then the variational Bayes posterior converges to a Gaussian with a mean centred at the true model parameter. Furthermore, since variational Bayes rests on optimization, variational inference easily takes advantage of methods such as stochastic optimization (Robbins and Monro Reference Robbins and Monro1951, Kushner and Yin Reference Kushner and Yin1997) and distributed optimization (though some MCMC methods can also exploit these innovations (Welling and Teh Reference Welling and Teh2011, Ahmed et al. Reference Ahmed, Aly, Gonzalez, Narayanamurthy and Smola2012)). Thus, variational inference is suited to large data sets and scenarios where we want to quickly explore many models; MCMC is suited to smaller data sets and scenarios where we are happy to pay a heavier computational cost for more precise samples. Another factor is the geometry of the posterior distribution. For example, the posterior of a mixture model admits multiple modes, each corresponding to label permutations of the components. Gibbs sampling, if the model permits, is a powerful approach to sampling from such target distributions; it quickly focuses on one of the modes. For mixture models where Gibbs sampling is not an option, variational inference may perform better than a more general MCMC technique (e.g. Hamiltonian Monte Carlo), even for small datasets (Küçükelbir, Ranganath, Gelman and Blei Reference Küçükelbir, Ranganath, Gelman and Blei2017).

4.3.1 Learning of TV-type regularizers and data fidelity terms

We start with a simple but theoretically and conceptually important example, namely the learning of total variation (TV)-type regularization models as proposed in De los Reyes and Schönlieb (Reference De los Reyes and Schönlieb2013), Kunisch and Pock (Reference Kunisch and Pock2013) and De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2017).

Let $X=BV(\unicode[STIX]{x1D6FA})$ , where $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{n}$ is a fixed open bounded set with Lipschitz boundary, and $Y=L^{2}(\mathbb{M},\mathbb{R})$ , where $\mathbb{M}$ is a manifold given by the data acquisition geometry. Next, let $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D706},\unicode[STIX]{x1D6FC})$ , where $\unicode[STIX]{x1D706}=(\unicode[STIX]{x1D706}_{1},\ldots ,\unicode[STIX]{x1D706}_{M})$ and $\unicode[STIX]{x1D6FC}=(\unicode[STIX]{x1D6FC}_{1},\ldots ,\unicode[STIX]{x1D6FC}_{N})$ are non-negative scalar parameters. We also assume that the loss $\ell _{X}:X\times X\rightarrow \mathbb{R}$ is convex, proper and weak* lower semicontinuous. We next study the bilevel problem

(4.4)

In the above, ${\mathcal{L}}_{\unicode[STIX]{x1D703}}:Y\times Y\rightarrow \mathbb{R}$ and ${\mathcal{S}}_{\unicode[STIX]{x1D703}}:X\rightarrow \mathbb{R}$ are defined for $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D706}_{i},\unicode[STIX]{x1D6FC}_{i})$ as

(4.5) $$\begin{eqnarray}{\mathcal{L}}_{\unicode[STIX]{x1D703}}({\mathcal{A}}(f),g)=\mathop{\sum }_{i=1}^{M}\unicode[STIX]{x1D706}_{i}{\mathcal{L}}_{i}({\mathcal{A}}(f),g)\quad \text{and}\quad {\mathcal{S}}_{\unicode[STIX]{x1D703}}(f)=\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\Vert {\mathcal{J}}_{j}(f)\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}\end{eqnarray}$$

where ${\mathcal{L}}_{i}:Y\times Y\rightarrow \mathbb{R}$ and ${\mathcal{J}}_{j}\in {\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})$ are given. Hence, in (4.4) the variational model is parametrized in terms of sums of different fidelity terms ${\mathcal{L}}_{i}$ and TV-type regularizers $\Vert {\mathcal{J}}_{j}(f)\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}$ , weighted against each other with parameters $\unicode[STIX]{x1D706}_{i}$ and $\unicode[STIX]{x1D6FC}_{j}$ (respectively). For $N=1$ and $S_{1}=D$ the distributional derivative, then

$$\begin{eqnarray}\Vert S_{1}(f)\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{1}})}=\operatorname{TV}(f).\end{eqnarray}$$

This framework is the basis for the analysis of the learning model, in which convexity of the variational model and compactness properties in the space of functions of bounded variation are crucial for proving existence of an optimal solution: see De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2016). Richer parametrizations for bilevel learning are discussed in Chen et al. (Reference Chen, Pock and Bischof2012, Reference Chen, Yu and Pock2015), for example, where non-linear functions and convolution kernels are learned. Chen et al., however, treat the learning model in finite dimensions, and a theoretical investigation of these more general bilevel learning models in a function space setting is a matter for future research.

In order to derive sharp optimality conditions for optimal parameters of (4.4) more regularity on the lower-level problem is needed. For shifting the problem (4.4) into a more regular setting, the Radon norms are regularized with Huber regularization and a convex, proper and weak* lower-semicontinuous smoothing functional $\operatorname{H}:X\rightarrow [0,\infty ]$ is added to the lower-level problem, typically $\operatorname{H}(f)=\frac{1}{2}\Vert \unicode[STIX]{x1D6FB}f\Vert ^{2}$ . In particular, the former is required for the single-valued differentiability of the solution map $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D6FC})\mapsto f_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D706}}$ , required by current numerical methods, irrespective of whether we are in a function space setting (see e.g. Rockafellar and Wets Reference Rockafellar and Wets1998, Theorem 9.56, for the finite-dimensional case). For parameters $\unicode[STIX]{x1D707}\geq 0$ and $\unicode[STIX]{x1D6FE}\in (0,\infty ]$ , the lower-level problem in (4.4) is then replaced by

(4.6)

with

$$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D703}}^{\unicode[STIX]{x1D6FE}}(f):=\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D6FC}_{j}\Vert {\mathcal{J}}_{j}(f)\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}^{\unicode[STIX]{x1D6FE}}.\end{eqnarray}$$

Here, $f\mapsto \Vert {\mathcal{J}}_{j}(f)\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}^{\unicode[STIX]{x1D6FE}}$ is the Huberized TV measure, given as follows.

Definition 4.1. Given $\unicode[STIX]{x1D6FE}\in (0,\infty ]$ , the Huber regularization for the norm $\Vert \,\cdot \,\Vert _{2}$ on $\mathbb{R}^{n}$ is defined by

$$\begin{eqnarray}\Vert g\Vert _{\unicode[STIX]{x1D6FE}}=\left\{\begin{array}{@{}ll@{}}\displaystyle \Vert g\Vert _{2}-\frac{1}{2\unicode[STIX]{x1D6FE}},\quad & \Vert g\Vert _{2}\geq \displaystyle \frac{1}{\unicode[STIX]{x1D6FE}},\\[9.0pt] \displaystyle \frac{\unicode[STIX]{x1D6FE}}{2}\Vert g\Vert _{2}^{2},\quad & \Vert g\Vert _{2}<\displaystyle \frac{1}{\unicode[STIX]{x1D6FE}}.\end{array}\right.\end{eqnarray}$$

Then, for $\unicode[STIX]{x1D707}\in {\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})$ with Lebesgue decomposition $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D708}L^{n}+\unicode[STIX]{x1D707}^{s}$ we have the Huber-regularized total variation measure,

$$\begin{eqnarray}|\unicode[STIX]{x1D707}|_{\unicode[STIX]{x1D6FE}}(V):=\int _{v}|\unicode[STIX]{x1D708}(x)|_{\unicode[STIX]{x1D6FE}}\text{d}x+|\unicode[STIX]{x1D707}^{s}|(V)\quad (V\subset \unicode[STIX]{x1D6FA}~\text{Borel-measurable}),\end{eqnarray}$$

and finally its Radon norm,

$$\begin{eqnarray}\Vert \unicode[STIX]{x1D707}\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}^{\unicode[STIX]{x1D6FE}}:=\Vert |\unicode[STIX]{x1D707}|_{\unicode[STIX]{x1D6FE}}\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{m_{j}})}.\end{eqnarray}$$

In all of these, we interpret the choice $\unicode[STIX]{x1D6FE}=\infty$ to give back the standard unregularized total variation measure or norm. In this setting existence of optimal parameters and differentiability of the solution operator can be proved, and with this an optimality system can be derived: see De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2016). More precisely, for the special case of the TV-denoising model the following theorem holds.

Theorem 4.2 (TV denoising (De los Reyes et al. Reference De los Reyes, Schönlieb and Valkonen2016)).

Consider the denoising problem (2.2) where $g,f_{\text{true}}\in BV(\unicode[STIX]{x1D6FA})\bigcap L^{2}(\unicode[STIX]{x1D6FA})$ , and assume $\operatorname{TV}(g)>\operatorname{TV}(f_{\text{true}})$ . Also, let $\operatorname{TV}^{\unicode[STIX]{x1D6FE}}(f):=\Vert Df\Vert _{{\mathcal{M}}(\unicode[STIX]{x1D6FA};\mathbb{R}^{n})}^{\unicode[STIX]{x1D6FE}}$ . Then there exist $\bar{\unicode[STIX]{x1D707}},\bar{\unicode[STIX]{x1D6FE}}>0$ such that any optimal solution $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D707}}\in [0,\infty ]$ to the problem

satisfies $\unicode[STIX]{x1D6FC}_{\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D707}}>0$ whenever $\unicode[STIX]{x1D707}\in [0,\bar{\unicode[STIX]{x1D707}}]$ and $\unicode[STIX]{x1D6FE}\in [\bar{\unicode[STIX]{x1D6FE}},\infty ]$ .

Theorem 4.2 states that if $g$ is a noisy image which oscillates more than the noise-free image $f_{\text{true}}$ , then the optimal parameter is strictly positive, which is exactly what we would naturally expect. De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2016) proved a similar result for second-order TGV and ICTV regularization for the case when $X=Y$ . The result was not extended to data with a general $Y$ , but it is possible with additional assumptions on the parameter space.

Moreover, in much of the analysis for (4.4) we could allow for spatially dependent parameters $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D706}$ . However, the parameters would then need to lie in a finite-dimensional subspace of $C_{0}(\unicode[STIX]{x1D6FA};\mathbb{R}^{N})$ : see De los Reyes and Schönlieb (Reference De los Reyes and Schönlieb2013) and Van Chung et al. (Reference Van Chung, De los Reyes and Schönlieb2017). Observe that Theorem 4.2 allows for infinite parameters $\unicode[STIX]{x1D6FC}$ . Indeed, for regularization parameter learning it is important not to restrict the parameters to be finite, as this allows us to decide between TGV $^{2}$ , TV and TV $^{2}$ regularization. De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2016) also prove a result on the approximation properties of the bilevel scheme with the smoothed variational model (4.6) as $\unicode[STIX]{x1D6FE}{\nearrow}\infty$ and $\unicode[STIX]{x1D707}{\searrow}0$ . In particular, they prove that as the numerical regularization vanishes, any optimal parameters for the regularized models tend to optimal parameters of the original model (4.4) in the sense of $\unicode[STIX]{x1D6E4}$ -convergence. Moreover, De los Reyes and Schönlieb (Reference De los Reyes and Schönlieb2013) take the limit as $\unicode[STIX]{x1D6FE}{\nearrow}\infty$ but $\unicode[STIX]{x1D707}>0$ fixed for the optimality system and derive an optimality system for the non-smooth problem. Recently Davoli and Liu (Reference Davoli and Liu2018) expanded the analysis for bilevel optimization of total variation regularizers beyond (4.4) to a richer family of anisotropic total variation regularizers in which the parameter of the dual norm and the (fractional) order of derivatives becomes part of the parameter space that (4.4) is optimized over.

Figure 4.1. Parameter optimality for TV denoising in Theorem 4.2. The non-convexity of the loss function, even for this one-parameter optimization problem, is clearly visible. Courtesy of Pan Liu.

Computing optimal solutions to the bilevel learning problems requires a proper characterization of optimal solutions in terms of a first-order optimality system. Since (4.4)–(4.6) constitutes a PDE-constrained optimization problem, suitable techniques from this field may be utilized. For the limit cases, an additional asymptotic analysis needs to be performed in order to get a sharp characterization of the solutions as $\unicode[STIX]{x1D6FE}\rightarrow \infty$ or $\unicode[STIX]{x1D707}\rightarrow 0$ , or both.

Several instances of the abstract problem (4.4) have been considered in the literature. De los Reyes and Schönlieb (Reference De los Reyes and Schönlieb2013) considered the case with TV regularization in the presence of several noise models. They proved the Gâteaux differentiability of the solution operator, which led to the derivation of an optimality system. Thereafter they carried out an asymptotic analysis with respect to $\unicode[STIX]{x1D6FE}\rightarrow \infty$ (with $\unicode[STIX]{x1D707}>0$ ), obtaining an optimality system for the corresponding problem. In that case the optimization problem corresponds to one with variational inequality constraints and the characterization concerns C-stationary points. Also, De los Reyes et al. (Reference De los Reyes, Schönlieb and Valkonen2017) have investigated differentiability properties of higher-order regularization solution operators. They proved a stronger Fréchet-differentiability result for the TGV $^{2}$ case, which also holds for TV. In particular, these stronger results open the door to further necessary and sufficient optimality conditions. Further, using the adjoint optimality condition gradient, formulas for the reduced cost functional can be derived, which in turn feed into the design of numerical algorithms for solving (4.4): see Calatroni et al. (Reference Calatroni, Cao, De los Reyes, Schönlieb and Valkonen2016, Section 3).

4.3.2 From fields of experts to variational networks

The TV-type regularization approaches discussed above used popular hand-crafted variational regularizers in a first attempt to make them more data-driven, staying in the framework of optimizing proposed regularization models in infinite-dimensional function space. In what follows, we change our perspective to discrete variational regularization models (4.2) in which $X=\mathbb{R}^{n}$ , and the functional ${\mathcal{S}}(\cdot )$ is equipped with parametrizations that go beyond TV-type regularizers; however, these are not covered by the theory in the previous section.

MRF parametrizations.

A prominent example of such richer parametrizations of regularizers has its origins in MRF theory. MRFs, first introduced by Besag (Reference Besag1974) and then used by Cross and Jain (Reference Cross and Jain1983) for texture modelling, provide a probabilistic framework for learning image priors: see e.g. Zhu, Wu and Mumford (Reference Zhu, Wu and Mumford1998). Higher-order MRF models such as the celebrated Field of Experts (FoE) model (Roth and Black Reference Roth and Black2005) are the most popular ones, as they are able to express very complex image structures and yield good performance for various image processing tasks. The FoE model is learning a rich MRF image prior by using ideas from sparse coding (Olshausen and Field Reference Olshausen and Field1997), where image patches are approximated by a linear combination of learned filters, variations of which lead to the well-known principal component analysis (PCA) and independent component analysis (ICA), and the Product of Experts (PoE) model (Welling, Osindero and Hinton Reference Welling, Osindero and Hinton2003).

In the context of bilevel learning, the FoE model has been used by Chen, Pock, Ranftl and Bischof (Reference Chen, Pock, Ranftl and Bischof2013) and Chen, Ranftl and Pock (Reference Chen, Ranftl and Pock2014). Here, the regularizer is parametrized in the form

(4.7) $$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f) :=\mathop{\sum }_{i=1}^{N}\left[\unicode[STIX]{x1D6FC}_{i}\mathop{\sum }_{k=1}^{n}\unicode[STIX]{x1D70C}((J_{i}f)_{k})\right]\quad \text{with }\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D6FC}_{1},\ldots ,\unicode[STIX]{x1D6FC}_{N},J_{1},\ldots ,J_{N}),\end{eqnarray}$$

where $N$ is the number of filters $J_{i}$ which are sparse and implemented as a two-dimensional convolution with a kernel $k_{i}$ , that is, $J_{i}f$ is the digitized version of $k_{i}\ast f$ . In the FoE model we use the parametrization $J_{i}=\sum _{j=1}^{m}\unicode[STIX]{x1D6FD}_{i,j}B_{j}$ , for a given set of basis filters $\{B_{1},\ldots ,B_{m}\}$ with $B_{j}\in \mathbb{R}^{n\times n}$ . Moreover, the $\unicode[STIX]{x1D6FC}_{i}$ are non-negative parameters and the non-linear function $\unicode[STIX]{x1D70C}(z)=\log (1+z^{2})$ . The shape of $\unicode[STIX]{x1D70C}$ is motivated by statistics over natural images that were shown to roughly follow a Student $t$ -distribution (Huang and Mumford Reference Huang and Mumford1999). With these parametrizations, $\unicode[STIX]{x1D703}$ in (4.7) can be seen as $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D6FC}_{i},\unicode[STIX]{x1D6FD}_{i,j})$ for $i=1,\ldots ,n$ and $j=1,\ldots ,m$ , and one learns the optimal $\unicode[STIX]{x1D703}$ from supervised training data $(f_{i},g_{i})$ by minimizing the squared $L^{2}$ -loss:

$$\begin{eqnarray}\ell _{X}({\mathcal{R}}_{\unicode[STIX]{x1D703}}(g),f) :=\Vert {\mathcal{R}}_{\unicode[STIX]{x1D703}}(g)-f\Vert _{2}^{2}\quad \text{with }{\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X\text{ given by (4.4)}.\end{eqnarray}$$

Chen et al. (Reference Chen, Ranftl and Pock2014) investigated other non-linearities such as the non-smooth and non-convex $\unicode[STIX]{x1D70C}(z)=\log (1+|z|)$ , and the non-smooth but convex $\ell _{1}$ -norm $\unicode[STIX]{x1D70C}(z)=|z|$ . Their experiments in particular suggested that, while the two $\log$ -type non-linearities gave very similar regularization performance, the convex $\ell _{1}$ -norm regularizer clearly did worse. Moreover, Chen et al. (Reference Chen, Yu and Pock2015) parametrized the non-linearities $\unicode[STIX]{x1D70C}$ with radial basis functions, whose coefficients are learned as well. Earlier MRF-based bilevel learning schemes exist: see e.g. Samuel and Tappen (Reference Samuel and Tappen2009).

A further development is the variational networks introduced by Kobler, Klatzer, Hammernik and Pock (Reference Kobler, Klatzer, Hammernik and Pock2017) and Hammernik et al. (Reference Hammernik, Klatzer, Kobler, Recht, Sodickson, Pock and Knoll2018). The idea was to replace the above bilevel scheme with learning to optimize (Section 4.9) using a supervised loss, thereby leading to a learned iterative method. This will be discussed in more detail in Section 5.1.4, where the variational networks re-emerge in the framework of learned iterative schemes.

Patch-based local models.

A particular class of sparse coding models are those which impose sparsity on patches (segmented) that are extracted from the signal. Let $g={\mathcal{A}}(f_{\text{true}})+e$ and $\hat{f}\in X$ . Further, let $\mathbb{D}=\{\unicode[STIX]{x1D719}_{k}\}\subset X$ be a dictionary, and $P_{j}:X\rightarrow X$ with $P_{j}(f)=f|_{\unicode[STIX]{x1D6FA}_{j}}$ for $j=1,\ldots N$ patches of $f$ . Assume that

(4.8) $$\begin{eqnarray}\hat{f}=\mathop{\sum }_{j=1}^{N}P_{j}(f_{\text{true}}),\end{eqnarray}$$

where $P_{j}(f_{\text{true}})\in X$ is compressible with respect to $\mathbb{D}$ .

Some of the currently leading denoising methods are based on patch-based local models. Examples include K-SVD, which is a sparse coding approach on image patches (Elad and Aharon Reference Elad and Aharon2006), BM3D, which combines sparsity and self-similarity (Dabov, Foi, Katkovnik and Egiazarian Reference Dabov, Foi, Katkovnik and Egiazarian2007), and EPLL, which is a Gaussian mixture model of image patches (Zoran and Weiss Reference Zoran and Weiss2011). Other examples are NCSR (Dong, Zhang, Shi and Li Reference Dong, Zhang, Shi and Li2013), weighted nuclear norm regularization of image patches (WNNM) (Gu, Zhang, Zuo and Feng Reference Gu, Zhang, Zuo and Feng2014) and SSC-GSM, a non-local sparsity approach with a Gaussian scale mixture (Dong, Shi, Ma and Li Reference Dong, Shi, Ma and Li2015).

Sparse-land models are a particular subclass of patch-based models. Here, each patch is sparse with respect to some global dictionary. Examples include sparse coding approaches that are applied patch-wise (Elad and Aharon Reference Elad and Aharon2006, Dong, Zhang, Shi and Wu Reference Dong, Zhang, Shi and Wu2011, Mairal and Ponce Reference Mairal and Ponce2014, Romano and Elad Reference Romano and Elad2015, Sulam and Elad Reference Sulam and Elad2015). In sparse-land models $\hat{f}\in X$ is reconstructed by solving

(4.9) $$\begin{eqnarray}\min _{f\in X,\unicode[STIX]{x1D709}_{i}\in \unicode[STIX]{x1D6EF}}\{{\mathcal{L}}({\mathcal{A}}(f),g)+{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f,\unicode[STIX]{x1D709}_{1},\ldots ,\unicode[STIX]{x1D709}_{N})\},\end{eqnarray}$$

where

(4.10) $$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f,\unicode[STIX]{x1D709}_{1},\ldots ,\unicode[STIX]{x1D709}_{N}) :=\mathop{\sum }_{j=1}^{N}[\unicode[STIX]{x1D706}_{j}\Vert P_{j}(f)-{\mathcal{E}}_{\mathbb{D}}^{\ast }(\unicode[STIX]{x1D709}_{j})\Vert _{2}^{2}+\unicode[STIX]{x1D707}_{j}\Vert \unicode[STIX]{x1D709}_{j}\Vert _{p}^{p}],\end{eqnarray}$$

with $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D706}_{j},\unicode[STIX]{x1D707}_{j})_{j=1}^{N}\in (\mathbb{R}^{2})^{N}$ and ${\mathcal{E}}_{\mathbb{D}}^{\ast }:\unicode[STIX]{x1D6EF}\rightarrow X$ denoting the synthesis operator associated with the given dictionary $\mathbb{D}$ . Bai et al. (Reference Bai, Yan, Jia, Jiang, Wang and Mou2017) propose the following alternating scheme to solve the sparse-land reconstruction problem:

(4.11)

The advantage of these approaches over plain-vanilla dictionary learning is that sparse-land models are computationally more feasible. Sparse-land models are one example of a dictionary learning approach, which will be discussed in the next section.

4.4.1 Reconstruction and dictionary learning

In sparse models, there are three ways to determine the dictionary. First, the dictionary is specified analytically. Second, the dictionary is determined from example data (dictionary learning). Third (as in the sparse-land models discussed above) the dictionary is determined jointly while performing reconstruction (joint reconstruction and dictionary learning). In what follows, we particularly focus on the second and third options for determining the dictionary.

Joint reconstruction and dictionary learning.

The recent paper by Chambolle, Holler and Pock (Reference Chambolle, Holler and Pock2018) proposes a convex variational model for joint reconstruction and (convolutional) dictionary learning that applies to inverse problems where data are corrupted and/or incomplete. Chambolle et al. (Reference Chambolle, Holler and Pock2018) prove rigorous results on well-posedness and stability stated in the infinite-dimensional (non-digitized) setting.

Earlier approaches were less rigorous and considered the finite-dimensional setting with a patch-based dictionary (sparse-land model). One such example is adaptive dictionary-based statistical iterative reconstruction (ADSIR) (Zhang et al. Reference Zhang, Zhang, Li, Peng, Liu and Zheng2016), which adds the dictionary $\mathbb{D}=\{\unicode[STIX]{x1D719}_{i}\}\subset X$ as an unknown to the sparse-land model as in global dictionary-based statistical iterative reconstruction (GDSIR) (Xu et al. Reference Xu, Yu, Mou, Zhang, Hsieh and Wang2012; see also Chun, Zheng, Long and Fessler Reference Chun, Zheng, Long and Fessler2017). This results in the following problem:

(4.12) $$\begin{eqnarray}\min _{f\in X,\unicode[STIX]{x1D709}_{i}\in \unicode[STIX]{x1D6EF},\mathbb{D}}\{{\mathcal{L}}({\mathcal{A}}(f),g)+{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f,\unicode[STIX]{x1D709}_{1},\ldots ,\unicode[STIX]{x1D709}_{N},\mathbb{D})\},\end{eqnarray}$$

where

(4.13) $$\begin{eqnarray}{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f,\unicode[STIX]{x1D709}_{1},\ldots ,\unicode[STIX]{x1D709}_{N},\mathbb{D}):=\mathop{\sum }_{j=1}^{N}[\unicode[STIX]{x1D706}_{j}\Vert P_{j}(f)-{\mathcal{E}}_{\mathbb{D}}^{\ast }(\unicode[STIX]{x1D709}_{j})\Vert _{2}^{2}+\unicode[STIX]{x1D707}_{j}\Vert \unicode[STIX]{x1D709}_{j}\Vert _{p}^{p}],\end{eqnarray}$$

with $\unicode[STIX]{x1D703}=((\unicode[STIX]{x1D706}_{j},\unicode[STIX]{x1D707}_{j})_{j=1}^{N})\in (\mathbb{R}^{2})^{N}$ , and ${\mathcal{E}}_{\mathbb{D}}^{\ast }:\unicode[STIX]{x1D6EF}\rightarrow X$ denotes the synthesis operator associated with the dictionary $\mathbb{D}$ . Usually an alternating minimization scheme is used to optimize over the three variables in (4.12).

Dictionary learning.

Let $\ell _{X}:X\times X\rightarrow \mathbb{R}$ be a given loss function (e.g. the $\ell _{2}$ - or $\ell _{1}$ -norm). Further, let $f_{1}\ldots f_{N}\in X$ be a given unsupervised training data, $\mathbb{D}=\{\unicode[STIX]{x1D719}_{i}\}\subset X$ a dictionary, and ${\mathcal{E}}_{\mathbb{D}}^{\ast }:\unicode[STIX]{x1D6EF}\rightarrow X$ the synthesis operator given as ${\mathcal{E}}_{\mathbb{D}}^{\ast }(\unicode[STIX]{x1D709})=\sum _{i}\unicode[STIX]{x1D709}_{i}\unicode[STIX]{x1D719}_{i}$ for $\unicode[STIX]{x1D709}\in \unicode[STIX]{x1D6EF}$ . One approach in dictionary learning is based on a sparsity assumption and solves

(4.14)

such that $\Vert \unicode[STIX]{x1D709}_{i}\Vert _{0}\leq s$ for $i=1,\ldots ,N$ , for a given sparsity level $s$ . Alternatively, one looks for a dictionary that minimizes the total cost for representing signals in the training data while enforcing a constraint on the precision in the following way:

(4.15)

such that $\ell _{X}(f_{i},{\mathcal{E}}_{\mathbb{D}}^{\ast }(\unicode[STIX]{x1D709}_{i}))\leq \unicode[STIX]{x1D716}$ for $i=1,\ldots ,N$ . A unified formulation is given by the unconstrained problem

(4.16)

All three formulations are posed in terms of the $\ell _{0}$ -norm and are NP-hard to compute. This suggests the use of convex relaxation, by which (4.16) becomes

(4.17)

If $\mathbb{D}$ is fixed then the sum in (4.17) decouples and leads to the convex relaxation of the sparse coding problem in (2.20) for ${\mathcal{A}}=\text{id}$ .

In the finite-dimensional setting $X=\mathbb{R}^{m}$ and $\unicode[STIX]{x1D6EF}$ is replaced by $\mathbb{R}^{n}$ for some $n$ . Then the dictionary is $\mathbb{D}:=\{{\unicode[STIX]{x1D719}_{k}\}}_{k=1\ldots n}\subset \mathbb{R}^{m}$ represented by an $n\times m$ matrix $\mathbf{\mathsf{D}}$ leading to the synthesis operator becoming a mapping ${\mathcal{E}}_{\mathbf{\mathsf{D}}}^{\ast }:\mathbb{R}^{m}\rightarrow \mathbb{R}^{n}$ with ${\mathcal{E}}_{\mathbb{D}}^{\ast }(\unicode[STIX]{x1D709})=\mathbf{\mathsf{D}}\cdot \unicode[STIX]{x1D709}$ . In this setting (4.17) becomes

(4.18)

Here, if $\mathbf{\mathsf{D}}$ satisfies the RIP then the convex relaxation preserves the sparse solution (Candès et al. Reference Candès, Romberg and Tao2006). State-of-the-art dictionary learning algorithms are K-SVD (Aharon, Elad and Bruckstein Reference Aharon, Elad and Bruckstein2006), geometric multi-resolution analysis (GRMA) (Allard, Chen and Maggioni Reference Allard, Chen and Maggioni2012) and online dictionary learning (Mairal, Bach, Ponce and Sapiro Reference Mairal, Bach, Ponce and Sapiro2010). Most work on dictionary learning to date has been done in the context of denoising, i.e.  ${\mathcal{A}}=\text{id}$ ; see also Rubinstein et al. (Reference Rubinstein, Bruckstein and Elad.2010).

4.4.2 Convolutional sparse coding and convolutional dictionary learning

Dictionary learning in the context of sparse coding has been very popular and successful, with several seminal approaches arising in this field, as outlined in the previous section. However, there are still several issues with sparse coding related to the locality of learned structures and the computational effort needed. Sparse-land models (4.9), for instance, perform sparse coding over all the patches, and this tends to be a slow process.

The computational performance can be addressed by using learning to improve upon an optimization scheme (Section 4.9), for example the Learned Iterative Soft-Thresholding Algorithm (LISTA) (see Section 4.9.2) learns a finite number of unrolled Iterative Soft-Thresholding Algorithm (ISTA) iterates using unsupervised training data to match ISTA solutions (Gregor and LeCun Reference Gregor and LeCun2010). Moreover, learning a dictionary over each patch independently as in (4.12) cannot account for global information, e.g. shift-invariance in images. What is needed is a computational feasible approach that introduces further structure and invariances in the dictionary, e.g. shift-invariance, and that makes each atom dependent on the whole signal instead of just individual patches. In this realm convolutional dictionaries have been introduced. Here atoms are given by convolution kernels and act on signal features via convolution, that is, $\mathbf{\mathsf{D}}$ is a concatenation of Toeplitz matrices (a union of banded and circulant matrices). Set up like this, convolutional dictionaries render computationally feasible shift-invariant dictionaries, where atoms depend on the entire signal.

Convolutional sparse coding.

Consider now the inverse problem of recovering $f_{\text{true}}\in X$ from (2.1) with the assumption that $f_{\text{true}}$ is compressible with respect to convolution dictionary $\mathbb{D}:=\{\unicode[STIX]{x1D719}_{i}\}\subset X$ .

In convolutional sparse coding (CSC), this is done by performing a synthesis using convolutional dictionaries, that is, atoms act by convolutions. More precisely, the reconstruction operator ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ is given as

(4.19) $$\begin{eqnarray}{\mathcal{R}}(g):=\mathop{\sum }_{i}\hat{\unicode[STIX]{x1D709}}_{i}\ast \unicode[STIX]{x1D719}_{i},\end{eqnarray}$$

where

Computational methods for solving (4.19) for denoising use convex relaxation followed by the alternating direction method of multipliers (ADMM) in frequency space (Bristow, Eriksson and Lucey Reference Bristow, Eriksson and Lucey2013) and its variants. See also Sreter and Giryes (Reference Sreter and Giryes2017) on using Learned Iterative Soft-Thresholding Algorithm (LISTA) in this context. So far CSC has only been analysed in the context of denoising (Bristow et al. Reference Bristow, Eriksson and Lucey2013, Wohlberg Reference Wohlberg2014, Gu et al. Reference Gu, Zuo, Xie, Meng, Feng and Zhang2015, Garcia-Cardona and Wohlberg Reference Garcia-Cardona and Wohlberg2017) with theoretical properties given in Papyan, Sulam and Elad (Reference Papyan, Sulam and Elad2016a , Reference Papyan, Sulam and Elad2016b ).

Convolutional dictionary learning.

Learning a dictionary in the context of CSC is called convolutional dictionary learning. Here, given unsupervised training data $f_{1},\ldots ,f_{m}\in X$ and a loss function $\ell _{X}:X:X\rightarrow X$ , one solves

(4.20)

where $\Vert \unicode[STIX]{x1D719}_{i}\Vert _{2}=1$ . For instance, Garcia-Cardona and Wohlberg (Reference Garcia-Cardona and Wohlberg2017) solved (4.20) with a squared $L_{2}$ -loss using an ADMM-type scheme. Extension of convolutional dictionary learning to a supervised data setting has been considered by Affara, Ghanem and Wonka (Reference Affara, Ghanem and Wonka2018), for instance. Here, discriminative dictionaries instead of purely reconstructive ones are learned by introducing a supervised regularization term in the usual CSC objective that encourages the final dictionary elements to be discriminative.

Multi-layer convolutional sparse coding.

A multi-layer extension of CSC, referred to as multi-layer convolutional sparse coding (ML-CSC), is proposed by Sulam, Papyan, Romano and Elad (Reference Sulam, Papyan, Romano and Elad2017). Given $L$ convolutional dictionaries $\mathbb{D}_{1},\ldots ,\mathbb{D}_{L}\subset X$ with atoms $\mathbb{D}_{k}:=\{\unicode[STIX]{x1D719}_{j}^{k}\}_{j}$ , a model parameter $f\in X$ admits a representation in terms of the corresponding multi-layer convolutional sparse coding (ML-CSC) model if there are $s_{1},\ldots ,s_{L}\in \mathbb{R}$ such that

$$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}f=\mathop{\sum }_{j}(\unicode[STIX]{x1D709}_{j}^{1}\ast \unicode[STIX]{x1D719}_{j}^{1}),\quad & \\[4.0pt] \boldsymbol{\unicode[STIX]{x1D709}}^{k}=\mathop{\sum }_{j}(\unicode[STIX]{x1D709}_{j}^{k+1}\ast \unicode[STIX]{x1D719}_{j}^{k+1})\quad \text{for }k=1,\ldots ,L-1,\quad & \end{array}\right.\end{eqnarray}$$

and $\Vert \boldsymbol{\unicode[STIX]{x1D709}}^{k}\Vert _{0,\infty }\leq s_{k}$ for $k=1,\ldots ,L$ . Hence, atoms $\unicode[STIX]{x1D719}_{k,i}\in \mathbb{D}_{k}$ in the $k$ th convolution dictionary are compressible in the $(k+1)$ th dictionary $\mathbb{D}_{k+1}$ for $k=1,\ldots ,L-1$ .

The ML-CSC model is a special case of CSC where intermediate representations have a specific structure (Sulam et al. Reference Sulam, Papyan, Romano and Elad2017, Lemma 1). Building on the theory for CSC, Sulam et al. (Reference Sulam, Papyan, Romano and Elad2017) provide a theoretical study of this novel model and its associated pursuits for dictionary learning and sparse coding in the context of denoising. Further, consequences for the theoretical analysis of partial differential equations (CNNs) can be extracted from ML-CSC using the fact that the resulting layered thresholding algorithm and the layered basis pursuit share many similarities with a forward pass of a deep CNN.

Indeed, Papyan, Romano and Elad (Reference Papyan, Romano and Elad2017) show that ML-CSC yields a Bayesian model that is implicitly imposed on $\hat{f}$ when deploying a CNN, and that consequently characterizes signals belonging to the model behind a deep CNN. Among other properties, one can show that the CNN is guaranteed to recover an estimate of the underlying representations of an input signal, assuming these are sparse in a local sense (Papyan et al. Reference Papyan, Romano and Elad2017, Theorem 4) and the recovery is stable (Papyan et al. Reference Papyan, Romano and Elad2017, Theorems 8 and 10). Many of these results also hold for fully connected networks, and they can be used to formulate new algorithms for partial differential equations (CNNs), for example to propose an alternative to the commonly used forward pass algorithm in CNN. This is related to both deconvolutional (Zeiler, Krishnan, Taylor and Fergus Reference Zeiler, Krishnan, Taylor and Fergus2010, Pu et al. Reference Pu, Yuan, Stevens, Li and Carin2016) and recurrent networks (Bengio, Simard and Frasconi Reference Bengio, Simard and Frasconi1994).

An essential technique for proving the key results in the cited references for ML-CSC is based on unrolling, which establishes a link between sparsity-promoting regularization (compressed sensing) and deep neural networks. More precisely, one starts with a variational formulation like that in (4.20) and specifies a suitable iterative optimization scheme. In this setting one can prove several theoretical results, such as convergence, stability and error estimates. Next, one unrolls the truncated optimization iterates and identifies the updating between iterates as layers in a deep neural network (Section 4.9.1). The properties of this network, such as stability and convergence, can now be analysed using methods from compressed sensing.

4.4.3 TV-type regularizer with learned sparsifying transforms.

The idea here is to learn the underlying sparsifying transform in a TV-type regularizer. One such approach is presented by Chun, Zheng, Long and Fessler (Reference Chun, Zheng, Long and Fessler2017), who consider the case of a single sparsifying transform. This idea is further developed by Zheng, Ravishankar, Long and Fessler (Reference Zheng, Ravishankar, Long and Fessler2018), who consider a regularizer that is given by a union of sparsifying transforms (ULTRA), which act on image patches and quantify the sparsification error of each patch using its best-matched sparsifying transform: see Zheng et al. (Reference Zheng, Ravishankar, Long and Fessler2018, equation (3)). The resulting optimization problem is solved by an intertwined scheme that alternates between a CT image reconstruction step, calculated via a relaxed linearized augmented Lagrangian method with ordered subsets, and a sparse coding and clustering step, which simultaneously groups the training patches into a fixed number of clusters and learns a transform in each cluster along with the sparse coefficients (in the transform domain) of the patches.

A variant of Zheng et al. (Reference Zheng, Ravishankar, Long and Fessler2018) is considered by Ye, Ravishankar, Long and Fessler (Reference Ye, Ravishankar, Long and Fessler2018b ), who use the same pre-learned union of sparsifying transforms as in Zheng et al. (Reference Zheng, Ravishankar, Long and Fessler2018), but the alternating scheme also includes iteratively updating a quadratic surrogate functions for the data discrepancy. These variants of learned sparsifying transforms are all applied to low-dose CT image reconstruction. Finally, Chen et al. (Reference Chen, Xiang, Gong, Wang and Tan2018) use essentially the method of Zheng et al. (Reference Zheng, Ravishankar, Long and Fessler2018) applied to cone-beam CT.

The Plug-and-Play Prior ( $P^{3}$ ) method.

This approach, which was introduced in Venkatakrishnan, Bouman and Wohlberg (Reference Venkatakrishnan, Bouman and Wohlberg2013), is based on the observation that a split operator method for minimizing the objective in (2.7) can be posed as an equality constrained problem, that is,

is equivalent to

The latter can be solved using ADMM (Section A.4), where the update in $h$ is computed using a proximal operator,

(4.21) $$\begin{eqnarray}h^{k+1 )}=\operatorname{prox}_{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D706}{\mathcal{S}}}(h^{(k)}-f+u),\end{eqnarray}$$

where $u$ is a Lagrange (dual) variable.

The idea is now to replace the proximal operator with a generic denoising operator, which implies that the regularization functional ${\mathcal{S}}$ is not necessarily explicitly defined. This opens the door to switching in any demonstrably successful denoising algorithms without redesigning a reconstruction algorithm – hence the name ‘Plug-and-Play’. However, the method comes with some disadvantages. The lack of an explicit representation of the regularizer detracts from a strict Bayesian interpretation of the regularizer as a prior probability distribution in MAP estimation and prevents explicit monitoring of the change in posterior probability of the iterative estimates of the solution. Next, the method is by design tied to the ADMM iterative scheme, which may not be optimal and which requires a non-trivial tuning of parameters of the ADMM algorithm itself (e.g. the Lagrangian penalty weighting term). Finally, it is not provably convergent for arbitrary denoising ‘engines’.

Regularization by denoising (RED).

The RED method (Romano, Elad and Milanfar Reference Romano, Elad and Milanfar2017a ) is motivated by the $P^{3}$ method. It is a variational method where the reconstruction operator is given as in (2.7), where the regularization functional is explicitly given as

(4.22) $$\begin{eqnarray}{\mathcal{S}}(f):=\langle f,f-\unicode[STIX]{x1D6EC}(f)\rangle \quad \text{for some }\unicode[STIX]{x1D6EC}:X\rightarrow X.\end{eqnarray}$$

The $\unicode[STIX]{x1D6EC}$ operator above is a general (non-linear) denoising operator: it can for example be a trained deep neural network. It does, however, need to satisfy two two key properties, which are justifiable in terms of the desirable features of an image denoiser.

Local homogeneity. The denoising operator should locally commute with scaling, that is,

$$\begin{eqnarray}\unicode[STIX]{x1D6EC}(cf)=c\unicode[STIX]{x1D6EC}(f)\end{eqnarray}$$

for all $f\in X$ and $|c-1|\leq \unicode[STIX]{x1D716}$ with $\unicode[STIX]{x1D716}$ small.Footnote ⁵

Strong passivity. The derivative of $\unicode[STIX]{x1D6EC}$ should have a spectral radius less than unity:

$$\begin{eqnarray}\unicode[STIX]{x1D70C}(\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}(f))\leq 1.\end{eqnarray}$$

This is justified by imposing a condition that the effect of the denoiser should not increase the norm of the model parameter:

$$\begin{eqnarray}\Vert \unicode[STIX]{x1D6EC}(f)\Vert =\Vert \unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}(f)f\Vert \leq \Vert \unicode[STIX]{x1D70C}(\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}(f))\Vert \,\Vert f\Vert \leq \Vert f\Vert .\end{eqnarray}$$

The key implication from local homogeneity is that the directional derivative of $\unicode[STIX]{x1D6EC}$ along $f$ is just the application of the denoiser to the $f$ itself:

$$\begin{eqnarray}\unicode[STIX]{x2202}\unicode[STIX]{x1D6EC}(f)f=\unicode[STIX]{x1D6EC}(f).\end{eqnarray}$$

The key implication from strong passivity is that it allows for convergence of the proposed RED methods by ensuring convexity of their associated regularization functionals.

Defining $W:X\rightarrow \mathbb{R}$ implicitly through the relation $W(f)f=\unicode[STIX]{x1D6EC}(f)$ and assuming local homogeneity and strong passivity yields the following computationally feasible expression for the gradient of the regularizer:

(4.23) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}{\mathcal{S}}(f)=f-\unicode[STIX]{x1D6EC}(f)=(\text{id}-W(f))(f).\end{eqnarray}$$

Here, the operator $(\text{id}-W(f))$ is further interpreted as an ‘image adaptive Laplacian-based’ regularizer. The above allows one to implement the RED framework in any optimization such as gradient-descent, fixed-point or ADMM in contrast to the $P^{3}$ approach, which is coupled to the ADMM scheme. See also Reehorst and Schniter (Reference Reehorst and Schniter2018) for further clarifications and new interpretations of the regularizing properties of the RED method.

4.7.1 Adversarial regularizer

A recent proposal by Lunz et al. (Reference Lunz, Öktem, Schönlieb and Bengio2018) for the construction of data-driven regularizers is inspired by how discriminative networks are trained using modern generative adversarial network (GAN) architectures. Our aim is to learn a regularizer ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ , which in some cases (Section 3.3.2) can be interpreted as being proportional to the negative log-likelihood of the prior $\unicode[STIX]{x1D6F1}_{\text{prior}}$ in a MAP estimator.

Consider the statistical setting in (3.4) and let $f_{i}\in X$ be samples of the $X$ -valued random variable $\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D6F1}_{\text{prior}}$ . Likewise, let $g_{i}\in X$ be samples of the $Y$ -valued random variable $\unicode[STIX]{x1D558}\sim \unicode[STIX]{x1D70E}$ that are independent of the samples $f_{i}$ , that is, we have unmatched training data. We also assume there exists a (potentially regularizing) pseudo-inverse $\{^{\dagger }:Y\rightarrow X$ to the forward operator ${\mathcal{A}}$ in (3.2) and define the measure $\unicode[STIX]{x1D70C}\in \mathscr{P}_{X}$ as $\unicode[STIX]{x1D70C}:=\{_{\#}^{\dagger }(\unicode[STIX]{x1D70E})$ for $\unicode[STIX]{x1D70E}\in \mathscr{P}(Y)$ . Note that both $\unicode[STIX]{x1D6F1}_{\text{prior}}$ and $\unicode[STIX]{x1D70E}$ are replaced by their empirical counterparts given by the training data $f_{i}\in X$ and $g_{i}\in Y$ , respectively.

The idea is to train the regularizer ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ parametrized by a neural network (see Section 7.4 for an exemplar architecture and application of this approach) in order to discriminate between the distributions $\unicode[STIX]{x1D6F1}_{\text{prior}}$ and $\unicode[STIX]{x1D70C}$ , the latter representing the distribution of imperfect solutions $\{^{\dagger }(g_{i})$ . More specifically, we compute

(4.24) $$\begin{eqnarray}{\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}:X\rightarrow \mathbb{R},\end{eqnarray}$$

where $\widehat{\unicode[STIX]{x1D703}}\in \text{arg}\,\text{min}_{\unicode[STIX]{x1D703}}L(\unicode[STIX]{x1D703})$ , with the loss function $\unicode[STIX]{x1D703}\mapsto L(\unicode[STIX]{x1D703})$ defined as

(4.25) $$\begin{eqnarray}L(\unicode[STIX]{x1D703}):=\mathbb{E}_{\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D6F1}_{\text{prior}}}[G_{1}({\mathcal{S}}_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D557}))]-\mathbb{E}_{\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D70C}}[G_{2}({\mathcal{S}}_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D557}))].\end{eqnarray}$$

Here, $G_{1},G_{2}:\mathbb{R}\rightarrow \mathbb{R}$ are monotone functions that have to be chosen. Popular choices for $G_{i}$ are logarithms, as in the original GAN paper (Goodfellow et al. Reference Goodfellow, Pouget-Abadie, Mirza, Xu, Warde-Farley, Ozair, Courville, Bengio and Ghahramani2014), and the $G_{i}$ associated with the Wasserstein loss, as in Arjovsky, Chintala and Bottou (Reference Arjovsky, Chintala and Bottou2017) and Gulrajani et al. (Reference Gulrajani, Ahmed, Arjovsky, Dumoulin, Courville and Guyon2017). The heuristic behind this choice for the loss function for training a regularizer is that a network trained with (4.25) will penalize noise and artefacts generated by the pseudo-inverse (and contained in $\unicode[STIX]{x1D70C}$ ). When used as a regularizer, it will hence prevent these undesirable features from occurring. Note also that in practical applications, the measures $\unicode[STIX]{x1D6F1}_{\text{prior}},\unicode[STIX]{x1D70C}\in \mathscr{P}_{X}$ are replaced with their empirical counterparts $\widehat{\unicode[STIX]{x1D6F1}}_{\text{prior}}$ and $\hat{\unicode[STIX]{x1D70C}}$ , given by training data $f_{i}$ and $\{^{\dagger }(g_{i})$ , respectively. The training problem in (4.24) for computing $\widehat{\unicode[STIX]{x1D703}}$ then reads as

(4.26)

We also point out that, unlike other data-driven approaches for inverse problems, the above method can be adapted to work with only unsupervised training data. A special case is to have $g_{i}\approx {\mathcal{A}}(f_{i})$ , which gives a unsupervised formulation of (4.26).

Lunz et al. (Reference Lunz, Öktem, Schönlieb and Bengio2018) chose a Wasserstein-flavoured loss functional (Gulrajani et al. Reference Gulrajani, Ahmed, Arjovsky, Dumoulin, Courville and Guyon2017) to train the regularizer, that is, one solves (4.24) with the loss function

(4.27) $$\begin{eqnarray}L(\unicode[STIX]{x1D703}):=\mathbb{E}_{\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D6F1}_{\text{prior}}}[{\mathcal{S}}_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D557})]-\mathbb{E}_{\unicode[STIX]{x1D557}\sim \unicode[STIX]{x1D70C}}[{\mathcal{S}}_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D557})]+\unicode[STIX]{x1D706}\mathbb{E}[(\Vert \unicode[STIX]{x1D6FB}{\mathcal{S}}_{\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D557})\Vert -1)_{+}^{2}].\end{eqnarray}$$

The last term in the loss function serves to enforce the trained regularizer ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ to be Lipschitz-continuous with constant one (Gulrajani et al. Reference Gulrajani, Ahmed, Arjovsky, Dumoulin, Courville and Guyon2017). Under appropriate assumptions on $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D70B}$ (see Assumptions 4.4 and 4.5) and for the asymptotic case of ${\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}$ having been trained to perfection, the loss (4.27) coincides with the 1-Wasserstein distance defined in (B.2).

A list of qualitative properties of ${\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}$ can be proved: for example, Theorem 1 of Lunz et al. (Reference Lunz, Öktem, Schönlieb and Bengio2018) shows that under appropriate regularity assumptions on the Wasserstein distance between $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D70B}$ , starting from elements in $\unicode[STIX]{x1D70C}$ and taking a gradient descent step of ${\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}$ (which results in a new distribution $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D702}}$ ) strictly decreases the Wasserstein distance between the new distribution $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D702}}$ and $\unicode[STIX]{x1D70B}$ . This is a good indicator that using ${\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}$ as a variational regularization term, and consequently penalizing it, indeed introduces the highly desirable incentive to align the distribution of regularized solutions with the distribution of ground truth samples $\unicode[STIX]{x1D6F1}_{\text{prior}}$ . Another characterization of such a trained regularizer ${\mathcal{S}}_{\widehat{\unicode[STIX]{x1D703}}}$ using the Wasserstein loss in (B.2) is in terms of a distance function to a manifold of desirable solutions ${\mathcal{M}}$ . To do so, we make the following assumptions.

This assumption captures the intuition that real data lie in a lower-dimensional submanifold of $X$ , which is a common assumption for analysing adversarial networks (Goodfellow et al. Reference Goodfellow, Pouget-Abadie, Mirza, Xu, Warde-Farley, Ozair, Courville, Bengio and Ghahramani2014). Moreover, it is assumed that the distributions $\unicode[STIX]{x1D70B}$ can be recovered from the distribution $\unicode[STIX]{x1D70B}$ through an appropriate projection onto the data manifold ${\mathcal{M}}$ .

Note that Assumption 4.5 is weaker than assuming that any given $f$ can be recovered by projecting the pseudo-inverse of the corresponding $g$ back onto the data manifold. Assumption 4.5 can instead be considered as a low-noise assumption. These assumptions yield the following theorem.

Theorem 4.6 shows that, if ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ were trained to perfection, i.e. trained so that ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ solves (B.2), then ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ would be given by the $L^{2}$ -distance function to ${\mathcal{M}}$ . This is implicitly also done in the RED approach described in Section 4.6. Similarly, Wu, Kim, Fakhri and Li (Reference Wu, Kim, Fakhri and Li2017) learn a regularizer in a variational model given as in Wu et al. (Reference Wu, Kim, Fakhri and Li2017, equation (3)) from unsupervised data by means of a K-sparse auto-encoder. This yields a regularizer that minimizes the distance of the image to the data manifold.

Finally, a weak stability result for ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ is proved in the spirit of the classical theory provided in Engl et al. (Reference Engl, Hanke and Neubauer2000). Since ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ is not necessarily bounded from below, the $1$ -Lipschitz property of ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ is used instead to prove this stability result.

Theorem 4.7 (Lunz, Öktem and Schönlieb Reference Lunz, Öktem, Schönlieb and Bengio2018).

Under appropriate assumptions on ${\mathcal{A}}:X\rightarrow Y$ given in Lunz et al. (Reference Lunz, Öktem, Schönlieb and Bengio2018, Appendix A), the following holds. Consider a sequence $\{{g_{n}\}}_{n}\subset Y$ with $g_{n}\rightarrow g$ in the norm topology in $Y$ and let $\{f_{n}\}\subset X$ denote a sequence of corresponding minimizers of (2.7), that is,

Then $f_{n}$ has a subsequence that converges weakly to a minimizer of $f\mapsto \Vert {\mathcal{A}}(f)-g\Vert ^{2}+\unicode[STIX]{x1D706}{\mathcal{S}}_{\unicode[STIX]{x1D703}}(f)$ .

Theorem 4.7 constitutes a starting point for further investigation into the regularizing properties of adversarial regularizers ${\mathcal{S}}_{\unicode[STIX]{x1D703}}$ . Section 7.4 shows the application of the adversarial regularizer to CT reconstruction.

4.7.2 The neural network Tikhonov (NETT) approach

Another proposal for learning a regularizer in (2.7) by a neural network is given in Li et al. (Reference Li, Schwab, Antholzer and Haltmeier2018b ) and called the NETT approach. NETT is based on composing a pre-trained network $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}:X\rightarrow \unicode[STIX]{x1D6EF}$ with a regularization functional ${\mathcal{S}}:\unicode[STIX]{x1D6EF}\rightarrow [0,\infty ]$ , such that ${\mathcal{S}}\circ \unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}:X\rightarrow [0,\infty ]$ takes small values for desired model parameters and penalizes (larger values) model parameters with artefacts or other unwanted structures.

Here, $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}:X\rightarrow \unicode[STIX]{x1D6EF}$ and ${\mathcal{S}}:\unicode[STIX]{x1D6EF}\rightarrow [0,\infty ]$ and the deep neural network $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ is allowed to be a rather general network model, a typical choice for $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ being an auto-encoder network. Once trained, the reconstruction operator is given by

(4.28)

where ${\mathcal{J}}_{\unicode[STIX]{x1D703}}(f):={\mathcal{L}}({\mathcal{A}}(f),g)+\unicode[STIX]{x1D706}{\mathcal{S}}(\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}(f))$ .

The main focus of Li et al. (Reference Li, Schwab, Antholzer and Haltmeier2018b ) is a discussion on analytic conditions, which guarantees that the NETT approach is indeed a regularization method in the sense of functional analytic regularization. In particular, Li et al. discuss assumptions on ${\mathcal{S}}$ and $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ such that the functional analytic regularization theory of Grasmair et al. (Reference Grasmair, Haltmeier and Scherzer2008) can be applied. This theory requires a weakly lower semicontinuous and coercive regularization term, which obviously holds for many deep neural networks $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ with coercive activation functions. Accordingly, Li et al. (Reference Li, Schwab, Antholzer and Haltmeier2018b ) discuss replacing the usual ReLU activation function with leaky ReLU, defined with a small $\unicode[STIX]{x1D70F}>0$ as

$$\begin{eqnarray}\ell \operatorname{ReLU}_{\unicode[STIX]{x1D70F}}(s):=\max (\unicode[STIX]{x1D70F}s,s).\end{eqnarray}$$

For $s\rightarrow -\infty$ , leaky ReLU also tends to $-\infty$ , which in combination with the affine linear maps ${\mathcal{W}}$ in $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ yields a coercive and weakly lower semicontinuous regularization function ${\mathcal{S}}\circ \unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ for standard choices of ${\mathcal{S}}$ , such as weighted $\ell _{p}$ -norms ${\mathcal{S}}(\unicode[STIX]{x1D709})=\sum _{i}v_{i}|\unicode[STIX]{x1D709}_{i}|^{p}$ with uniformly positive weights $v_{i}$ and $p\geq 1$ . They even go beyond these classical results and, by introducing the novel concept of absolute Bregman distances (Li et al. Reference Li, Schwab, Antholzer and Haltmeier2018b ), they obtain convergence results and convergence rates in the underlying function space norm.

Li et al. (Reference Li, Schwab, Antholzer and Haltmeier2018b ) discuss the application of NETT to PAT reconstruction from limited data. The neural network $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ is trained against supervised data $(f_{i},h_{i})\in X\times X$ , where $f_{i}$ serves as ground truth model parameter and $h_{i}=(\{^{\dagger }\circ {\mathcal{A}})(f_{i})$ , where $\{^{\dagger }:Y\rightarrow X$ is some (regularized) pseudo-inverse of the forward operator ${\mathcal{A}}$ . In tomographic applications it is taken as FBP, so $h_{i}$ will typically contain sampling artefacts. The network $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ is then modelled as an auto-encoder, more specifically as the encoder part of an encoder–decoder neural network. The parameters $\unicode[STIX]{x1D703}$ are trained by minimizing a loss function that evaluates the capability of the encoder–decoder pair to reproduce desirable $f$ as well as $h=(\{^{\dagger }\circ {\mathcal{A}})(f)$ with artefacts using an appropriate distance function. After training $\unicode[STIX]{x1D6F9}_{\unicode[STIX]{x1D703}}$ as described above, the NETT functional (4.28) is minimized by a generalized gradient descent method. The numerical results obtained with artificial data confirm the theoretical findings but lack comparison or seem to be slightly less favourable when compared with results obtained with other neural network approaches: see e.g. Hauptmann et al. (Reference Hauptmann, Lucka, Betcke, Huynh, Adler, Cox, Beard, Ourselin and Arridge2018).

4.9.1 General principle of unrolling

The aim of unrolling is to find a deep neural network architecture ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ that is especially suited to approximating an operator ${\mathcal{R}}:Y\rightarrow X$ that is implicitly defined via an iterative scheme. A typical example is when ${\mathcal{R}}(g)$ is implicitly defined via an iterative scheme that is designed to converge to a (local) minimum of $f\mapsto {\mathcal{J}}(f,g)$ .

We start by giving an illustrative example of how to unroll an iterative gradient descent scheme. This is followed by an outline of the general principle of unrolling. Two specific cases are considered in Sections 4.9.2 and 4.9.3.

Unrolled gradient descent as a neural network.

Consider the case when ${\mathcal{J}}(\,\cdot \,,g):X\rightarrow \mathbb{R}$ is smooth for any $g\in Y$ and assume ${\mathcal{R}}(g)$ is given by a standard gradient descent algorithm, that is,

$$\begin{eqnarray}{\mathcal{R}}(g)=\lim _{k\rightarrow \infty }f^{k}\quad \text{where }\left\{\begin{array}{@{}ll@{}}f^{0}=f_{0}\quad \text{is given,}\quad & \\ f^{k}=f^{k-1}-\unicode[STIX]{x1D714}_{k}\unicode[STIX]{x1D6FB}_{f}{\mathcal{J}}(f^{k},g)\quad \text{for }k=1,2,\ldots .\quad & \end{array}\right.\end{eqnarray}$$

Each parameter $\unicode[STIX]{x1D714}_{k}$ is a step length for the $k$ th iteration, and normally these are chosen via the Goldstein rule or backtracking line search (Armijo rule) (Bertsekas Reference Bertsekas1999), which under suitable conditions ensures convergence to a minimum.

Computational feasibility of the above scheme is directly tied to the number of iterates necessary to get sufficiently close to the desired local minima. Time limitations imposed by many applications limit the maximum number of evaluations of $f\mapsto {\mathcal{J}}(f,g)$ , that is, we have an a priori bound $N$ on the number of iterates. Unrolling the gradient descent scheme defining ${\mathcal{R}}$ and stopping the iterates after $N$ steps allows us to express the $N$ th iterate as

(4.32) $$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D703}}(g):=(\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D714}_{N}}\circ \ldots \circ \unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D714}_{1}})(f_{0})\quad \text{with }\unicode[STIX]{x1D703}:=(\unicode[STIX]{x1D714}_{1},\ldots ,\unicode[STIX]{x1D714}_{N}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D714}_{k}}:X\rightarrow X$ are updating operators given by

(4.33) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D714}_{k}}:=\text{id}-\unicode[STIX]{x1D714}_{k}\unicode[STIX]{x1D6FB}\!{\mathcal{J}}(\,\cdot \,,g)\quad \text{for }k=1,\ldots ,N.\end{eqnarray}$$

Now, note that ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ can be seen as a feed-forward neural network where each layer in the network evaluates $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D714}_{k}}$ and the parameters of the network are $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D714}_{1},\ldots ,\unicode[STIX]{x1D714}_{N})$ . Moreover, if the step length is fixed, i.e.  $\unicode[STIX]{x1D714}_{1}=\cdots =\unicode[STIX]{x1D714}_{N}=\unicode[STIX]{x1D714}$ for some $\unicode[STIX]{x1D714}$ , the gradient descent algorithm can in fact be interpreted as a recurrent neural network. For both cases, the best choice of step lengths $\unicode[STIX]{x1D714}_{1},\ldots ,\unicode[STIX]{x1D714}_{N}$ for approximating ${\mathcal{R}}:Y\rightarrow X$ with $N$ gradient descent iterates can be obtained by unsupervised learning: simply solve (4.31) with ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ as in (4.32).

A further option is to replace the explicit updating in (4.33), used to define ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ in (4.32), with generic deep neural networks:

(4.34) $$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}:=\text{id}+\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}(\,\cdot \,,\unicode[STIX]{x1D6FB}\!{\mathcal{J}}(\,\cdot \,,g))\quad \text{where }\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}:X\times X\rightarrow X.\end{eqnarray}$$

Each $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}$ is now a deep neural network, for example a CNN, that is trained against unsupervised data by solving (4.31) with ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ as in (4.34).

Abstract unrolled schemes.

The above example of unrolling a gradient descent and replacing its explicit updates with a neural network trained by unsupervised learning applies to many other iterative schemes. It essentially amounts to optimizing over optimization solvers. Since we seek $g\mapsto {\mathcal{R}}(g)$ given by (4.29), we are not only interested in optimizing a single objective, but rather an infinite family $\{{\mathcal{J}}(\,\cdot \,,g)\}_{g}$ of objectives parametrized by data $g\in Y$ .

More precisely, if $f\mapsto {\mathcal{J}}(f,g)$ in (4.29) is smooth, then it is natural to consider an iterative scheme that makes use of the gradient of the objective. Given an initial model parameter $f^{0}\in X$ (usually set to zero), such schemes can be written abstractly as

$$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}f^{0}=f_{0}\in X\text{ chosen},\quad & \\ f^{k+1}:=\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}(f^{k},\mathbf{f}_{m}^{k},\unicode[STIX]{x1D6FB}_{\!f}\!{\mathcal{J}}(f^{k},g)),\quad \end{array}\right.\end{eqnarray}$$

for some updating operator $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}:X\times X^{m}\times X\rightarrow X$ . The above formulation includes accelerated schemes, like fast gradient methods (Nesterov Reference Nesterov2004) and quasi-Newton methods (Nocedal and Wright Reference Nocedal and Wright2006). Now, stopping and unrolling the above scheme after $N$ iterates amounts to defining ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ with $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N})$ as

$$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D703}}(g)=(\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{N}}\circ \ldots \circ \unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{1}})(f^{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}:X\times X^{m}\rightarrow X\times X^{m}$ is defined as

$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}(f,\mathbf{h}):=(\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}(f,\mathbf{h},\unicode[STIX]{x1D6FB}_{\!f}\!{\mathcal{J}}(f,g)),\mathbf{h}^{\prime })\quad \text{with }\mathbf{h}^{\prime }:=(h_{2},\ldots ,h_{m})\in X^{m-1}\end{eqnarray}$$

and ${\mathcal{P}}_{X}:X\times X^{m}\rightarrow X$ is the usual projection. The final step is to replace the hand-crafted updating operators $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}:X\times X^{m}\times X\rightarrow X$ with deep neural networks that then train the resulting ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ against unsupervised data as in (4.31).

In inverse problem applications, the objective in (4.29) often has further structure that can be utilized. A typical example is (2.7), where the objective is of the form

(4.35) $$\begin{eqnarray}{\mathcal{J}}(f,g):={\mathcal{L}}({\mathcal{A}}(f),g)+{\mathcal{S}}_{\unicode[STIX]{x1D706}}(f).\end{eqnarray}$$

A wide range of iterative schemes that better utilize such a structure include updating in both $X$ (primal) and $Y$ (dual) spaces. These can be written abstractly as

(4.36) $$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}(f^{0},g^{0})=(f_{0},g)\in X\times Y\text{ with }f_{0}\text{ chosen},\quad & \\ (f^{k+1},g^{k+1}):=\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}(f^{k},\mathbf{f}_{m_{1}}^{k},g^{k},\mathbf{g}_{m_{2}}^{k},[\unicode[STIX]{x2202}\!{\mathcal{A}}(f^{k})]^{\ast }(g^{k}),{\mathcal{A}}(f^{k}),\unicode[STIX]{x1D6FB}\!{\mathcal{S}}_{\unicode[STIX]{x1D706}}(f^{k}),)\quad \end{array}\right.\end{eqnarray}$$

for some updating operator

(4.37) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}:X\times X^{m_{1}}\times Y\times Y^{m_{2}}\times X\times Y\times X\rightarrow X\times Y.\end{eqnarray}$$

Here, $m_{1},m_{2}$ is the memory in the $X$ - and $Y$ -iterates, so

$$\begin{eqnarray}\mathbf{f}_{m_{1}}^{k}=(f^{k-1},\ldots ,f^{k-m_{1}})\in X^{m_{1}}\quad \text{and}\quad \mathbf{g}_{m_{2}}^{k}=(g^{k-1},\ldots ,g^{k-m_{2}})\in Y^{m_{2}},\end{eqnarray}$$

with the convention that $f^{k-l}=f^{0}$ and $g^{k-l}=g$ whenever $k-l<0$ . Stopping and unrolling such a scheme after $N$ iterates amounts to defining ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ with $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D703}_{1},\ldots ,\unicode[STIX]{x1D703}_{N})$ as

(4.38) $$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D703}}(g):=({\mathcal{P}}_{X}\circ \unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{N}}\circ \ldots \circ \unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{1}})(f^{0},\mathbf{f}_{m_{1}}^{0},g,\mathbf{g}_{m_{2}}^{0}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}:X\times X^{m_{1}}\times Y\times Y^{m_{2}}\rightarrow X\times X^{m_{1}}\times Y\times Y^{m_{2}}$ for $k=1,\ldots ,N$ is defined as

$$\begin{eqnarray}\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}(f,\mathbf{h},v,\mathbf{v}):=(\hat{f}_{\unicode[STIX]{x1D703}_{k}},\mathbf{h}^{\prime },\hat{v}_{\unicode[STIX]{x1D703}_{k}},\mathbf{v}^{\prime }),\end{eqnarray}$$

with $\mathbf{h}^{\prime }:=(h_{2},\ldots ,h_{m_{1}})\in X^{m_{1}-1}$ , $\mathbf{v}^{\prime }:=(v_{2},\ldots ,v_{m_{2}})\in X^{m_{2}-1}$ , and

$$\begin{eqnarray}(\hat{f}_{\unicode[STIX]{x1D703}_{k}},\hat{v}_{\unicode[STIX]{x1D703}_{k}}):=\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D703}_{k}}(f,\mathbf{f},v,\mathbf{v},[\unicode[STIX]{x2202}\!{\mathcal{A}}(f)]^{\ast }(v),{\mathcal{A}}(f),\unicode[STIX]{x1D6FB}\!{\mathcal{S}}_{\unicode[STIX]{x1D706}}(f))\in X\times Y.\end{eqnarray}$$

In the above, ${\mathcal{P}}_{X}:X\times X^{m_{1}}\times Y\times Y^{m_{2}}\rightarrow X$ is the usual projection and, just as before, these hand-crafted updating operators are replaced with deep neural networks and the resulting deep neural network ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ is trained against training data, for example as in (4.31) in the case with unsupervised data.

Likewise, in the setting where ${\mathcal{S}}_{\unicode[STIX]{x1D706}}:X\rightarrow \mathbb{R}$ in (4.35) is non-smooth, the prototype proximal-gradient scheme without memory reads as

(4.39) $$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}f^{k+1/2}:=\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}^{1}}(f^{k},\unicode[STIX]{x1D6FB}_{\!f}{\mathcal{L}}({\mathcal{A}}(f^{k}),g)),\quad & \\ f^{k+1}:=\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}^{2}}(f^{k},f^{k+1/2},\operatorname{prox}_{\unicode[STIX]{x1D6FE}{\mathcal{S}}_{\unicode[STIX]{x1D706}}}(f^{k+1/2})),\quad \end{array}\right.\end{eqnarray}$$

for updating operators $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}^{1}}:X\times X\rightarrow X$ and $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}^{2}}:X\times X\times X\rightarrow X$ . Now, just as before, one can terminate and unroll the above scheme after $N$ iterates and replace the updating operators with deep neural networks. The resulting deep neural network ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow X$ , where $\unicode[STIX]{x1D703}=(\unicode[STIX]{x1D703}_{1}^{1},\ldots ,\unicode[STIX]{x1D703}_{N}^{1},\unicode[STIX]{x1D703}_{1}^{2},\ldots ,\unicode[STIX]{x1D703}_{N}^{2})$ , is then trained against unsupervised data as in (4.31). Just as in the smooth case, it is furthermore possible to add memory and/or ‘break up’ the objective even further, the latter in order to better account for the structure in the problem. For example, variational methods for linear inverse problems often result in minimizing an objective of the form

$$\begin{eqnarray}{\mathcal{J}}(f,g):={\mathcal{S}}_{\unicode[STIX]{x1D706}}(f)+\mathop{\sum }_{i=1}^{m}{\mathcal{L}}_{i}({\mathcal{A}}_{i}(f),g_{i})\end{eqnarray}$$

where ${\mathcal{A}}_{i}:X\rightarrow Y_{i}$ are linear and ${\mathcal{S}}_{\unicode[STIX]{x1D706}}:X\rightarrow [-\infty ,\infty ]$ and ${\mathcal{L}}_{i}:Y_{i}\times Y_{i}\rightarrow [-\infty ,\infty ]$ are proper, convex and lower semicontinuous. One can then consider unrolling iterative schemes that utilize this structure, such as operator splitting techniques: see Eckstein and Bertsekas (Reference Eckstein and Bertsekas1992), Beck and Teboulle (Reference Beck and Teboulle2009), Chambolle and Pock (Reference Chambolle and Pock2011), Boyd et al. (Reference Boyd, Parikh, Chu, Peleato and Eckstein2011), Combettes and Pesquet (Reference Combettes, Pesquet and Bauschke2011, Reference Combettes and Pesquet2012), He and Yuan (Reference He and Yuan2012), Boţ and Hendrich (Reference Boţ and Hendrich2013), Boţ and Csetnek (Reference Boţ and Csetnek2015), Ko, Yu and Won (Reference Ko, Yu and Won2017), Latafat and Patrinos (Reference Latafat and Patrinos2017) and Bauschke and Combettes (Reference Bauschke and Combettes2017).

A natural question is to investigate the error in using ${\mathcal{R}}_{\unicode[STIX]{x1D703}}$ to approximate ${\mathcal{R}}$ as a function of $N$ and properties of the unsupervised training data. Such estimates, which are referred to as time–accuracy trade-offs, have been proved for LISTA (Section 4.9.2) by Giryes, Eldar, Bronstein and Sapiro (Reference Giryes, Eldar, Bronstein and Sapiro2017); see also Oymak and Soltanolkotabi (Reference Oymak and Soltanolkotabi2017) for further development along these lines. Another related type of investigation has been pursued by Banert et al. (Reference Banert, Ringh, Adler, Karlsson and Öktem2018), who derive conditions on the unrolled iterative scheme that is terminated after $N$ iterates, so that training this yields a scheme that is convergent in the limit. It turns out that imposing such convergence constraints has only a minor impact on the performance at $N$ iterates. Furthermore, it improves the generalization properties, that is, one can use the same trained solver even when the objective is slightly changed, for example by considering a different forward operator.

To summarize, the iterative update in most schemes for solving optimization problems has the structure of a composition of (possibly multiple) affine operations followed by a non-linear operation. If each iterate is identified with a layer and the non-linear operation plays the role of an activation function, e.g.  ${\mathcal{S}}_{\unicode[STIX]{x1D706}}(f)=\unicode[STIX]{x1D706}\Vert f\Vert$ , then $\operatorname{prox}_{\unicode[STIX]{x1D706}{\mathcal{S}}}$ corresponds to soft thresholding, which in turn is very close to the well-known ReLU activation function. In this way the iterative scheme stopped after $N$ steps can be represented via a neural network with an architecture that is adapted to that iterative scheme.

4.9.2 Learned Iterative Soft-Thresholding Algorithm (LISTA)

LISTA is the earliest example of unrolling an optimization scheme and it was first introduced in Gregor and LeCun (Reference Gregor and LeCun2010). It is the abstract unrolling scheme in Section 4.9.1 adapted to the specific case of ISTA iterates (Section A.3). This results in a fully connected network with $N$ -internal layers of identical size that is adapted to evaluating the solution operator for the following convex non-smooth optimization problem:

(4.40) $$\begin{eqnarray}{\mathcal{R}}(g):={\mathcal{E}}^{\ast }(\hat{\unicode[STIX]{x1D709}}),\end{eqnarray}$$

where

Such optimization problems arise as the convex relaxation of (2.20), which is sparsity-promoting regularization of an ill-posed linear inverse problem (Section 2.7).

The ISTA iterative scheme for evaluating ${\mathcal{R}}$ reads as

(4.41) $$\begin{eqnarray}\unicode[STIX]{x1D709}^{n+1}=\operatorname{prox}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}\Vert \cdot \Vert _{1}}=S_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}}[\unicode[STIX]{x1D70F}{\mathcal{A}}_{{\mathcal{E}}^{\ast }}^{\ast }g+(\text{id}-\unicode[STIX]{x1D70F}{\mathcal{A}}_{{\mathcal{E}}^{\ast }}^{\ast }{\mathcal{A}}_{{\mathcal{E}}^{\ast }})\unicode[STIX]{x1D709}^{n}],\end{eqnarray}$$

where $0<\unicode[STIX]{x1D70F}<2/L$ is the step length approximated by the reciprocal of the Lipschitz constant of ${\mathcal{A}}_{{\mathcal{E}}^{\ast }}$ , and $S$ is the shrinkage operator (see Section A.3). Now the insight is to recognize $b_{k}(g):=\unicode[STIX]{x1D70F}_{k}{\mathcal{A}}_{{\mathcal{E}}^{\ast }}^{\ast }g$ as a bias, $W_{k}:=(\text{id}-\unicode[STIX]{x1D70F}_{k}{\mathcal{A}}_{{\mathcal{E}}^{\ast }}^{\ast }{\mathcal{A}}_{{\mathcal{E}}^{\ast }})$ as a fixed linear operator and $\unicode[STIX]{x1D713}_{k}=S_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}_{k}}$ as a pointwise non-linearity. Note that $W_{k}$ is symmetric positive definite by construction. Assuming a fixed iteration number $N$ , we may write the approximation to (4.40) as

$$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D709}^{(0)}\in \unicode[STIX]{x1D6EF}\text{given}\quad & \\ \unicode[STIX]{x1D709}^{(k)}=\unicode[STIX]{x1D713}(W_{k}\unicode[STIX]{x1D709}^{(k-1)}+b_{k}(g))\quad & \end{array}\right.\quad \text{for }k=1,\ldots ,N.\end{eqnarray}$$

In unrolled form, the above defines ${\mathcal{R}}_{\unicode[STIX]{x1D703}}:Y\rightarrow \unicode[STIX]{x1D6EF}$ with $\unicode[STIX]{x1D703}_{k}:=(\unicode[STIX]{x1D713},\unicode[STIX]{x1D70F}_{k},W_{k},b_{k})$ as

(4.42) $$\begin{eqnarray}{\mathcal{R}}_{\unicode[STIX]{x1D703}}:=(\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{N}}\circ \ldots \circ \unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{1}})(\unicode[STIX]{x1D709}^{(0)})\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}:\unicode[STIX]{x1D6EF}\rightarrow \unicode[STIX]{x1D6EF}$ is $\unicode[STIX]{x1D6EC}_{\unicode[STIX]{x1D703}_{k}}:=\unicode[STIX]{x1D713}_{k}\circ (W_{k}+b_{k}(g)\text{id})$ . This is precisely the form of a feed-forward network with $N$ layers, each with identical fixed affine map ${\mathcal{W}}_{k}(\unicode[STIX]{x1D709}):=W_{k}\unicode[STIX]{x1D709}+b_{k}$ ; the activation function $\unicode[STIX]{x1D713}_{k}:=S_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D70F}_{k}}$ is given by the shrinkage operator (see Section A.3).

The affine maps between the layers are assumed to be identical. LISTA then trains $W$ and potentially $\unicode[STIX]{x1D713}$ on some given unsupervised training data by solving (4.31). The derivation of Gregor and LeCun (Reference Gregor and LeCun2010) is framed in terms of learning a sparse encoder. Here we rephrase the idea in terms of a reconstruction problem and a learned decoder as follows.