3.1 Total variation and related regularizations

As mentioned above, total variation regularization has been one of the driving examples for development of regularization methods in Banach spaces starting from Rudin et al. (Reference Rudin, Osher and Fatemi1992) and Acar and Vogel (Reference Acar and Vogel1994). Since then it has been a constant source of motivation for further development of mathematical analysis (e.g. Chambolle and Lions Reference Chambolle and Lions1997, Strong et al. Reference Strong and Chan1996, Scherzer Reference Scherzer1998, Chavent and Kunisch Reference Chavent and Kunisch1997, Ring Reference Ring2000, Strong and Chan Reference Strong and Chan2003, Burger and Osher Reference Burger and Osher2004, Caselles, Chambolle and Novaga Reference Caselles, Chambolle and Novaga2007, Allard Reference Allard2007), computational optimization techniques for non-smooth problems (e.g. Chan, Golub and Mulet Reference Chan, Golub and Mulet1999, Vogel Reference Vogel2002, Chambolle Reference Chambolle2004, Kunisch and Hintermüller Reference Kunisch and Hintermüller2004, Chambolle and Pock Reference Chambolle and Pock2011), and development of advanced models (e.g. Scherzer Reference Scherzer1998, Osher et al. Reference Osher, Burger, Goldfarb, Xu and Yin2005, Burger, Osher, Xu and Gilboa Reference Burger, Osher, Xu, Gilboa and Paragios2005, Burger et al. Reference Burger, Gilboa, Osher and Xu2006, Bredies and Holler Reference Bredies and Holler2015a , Hu and Jacob Reference Hu and Jacob2012, Lenzen, Becker and Lellmann Reference Lenzen, Becker, Lellmann and Aujol2013, Benning, Brune, Burger and Müller Reference Benning, Brune, Burger and Müller2013).

The key step for modern analysis and computational methods is the (pre-) dual formulation of total variation,

(3.12) $$\begin{eqnarray}TV(u)=|u|_{BV}:=\sup _{g\in C_{0}^{\infty }(\unicode[STIX]{x1D6FA})^{d},g\in {\mathcal{C}}}\int _{\unicode[STIX]{x1D6FA}}u\unicode[STIX]{x1D6FB}\cdot g\,\text{d}x,\end{eqnarray}$$

with the convex set

$$\begin{eqnarray}{\mathcal{C}}=\{g\in L^{\infty }(\unicode[STIX]{x1D6FA})\mid |g(x)|\leq 1\text{ a.e. in }\unicode[STIX]{x1D6FA}\}.\end{eqnarray}$$

This characterization allows us to understand the structure of subgradients as elements of ${\mathcal{C}}$ absolutely continuous with respect to the gradient measure $D$ such that

$$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}g\cdot \,\text{d}Du=|u|_{BV}.\end{eqnarray}$$

The optimality condition (3.3) provides

(3.13) $$\begin{eqnarray}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku,f)+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FB}\cdot g=0,\end{eqnarray}$$

where $g$ is a vector field such that $g|Du|$ is a polar decomposition of the vector measure (Ambrosio et al. Reference Ambrosio, Fusco and Pallara2000).

In one spatial dimension the structure of solutions can be understood directly from the optimality condition. If there is an open set where $u$ is not constant, either with positive or negative derivative, then $g$ equals $+1$ or $-1$ , hence its derivative vanishes. Thus, in such regions the generalized residual $K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku,f)$ vanishes. In the case of noisy data this is usually not happening for larger sets, thus $u$ is typically piecewise constant. In higher spatial dimensions this is not completely true, but still the case $|g(x)|<1$ is the canonical one, so in many cases solutions are piecewise constant. On the other hand, piecewise constant structures are not optimal in all instances. In particular, total variation methods are well known to exhibit staircasing phenomena, that is, smoothly varying parts in the solution are often approximated by piecewise constant structures with many jumps resembling a stair structure. For this reason many modifications and variants of total variation regularization have been investigated in recent decades. An immediate option is that of higher-order total variation approaches, which formally replace the one-norm of the gradient with the one-norm of a higher-order derivative such as the Laplacian, the Hessian or the symmetric part of the Hessian (e.g. Scherzer Reference Scherzer1998, Chan, Esedoglu and Park Reference Chan, Esedoglu and Park2010, Hinterberger and Scherzer Reference Hinterberger and Scherzer2006, Papafitsoros and Schönlieb Reference Papafitsoros and Schönlieb2014). The disadvantage of such an approach is that solutions of the regularization model will be too regular and discontinuity sets (edges) are lost. In view of (3.12) such approaches can be characterized by ${\mathcal{C}}$ not being a bounded set in $L^{\infty }(\unicode[STIX]{x1D6FA})$ , but rather being derivatives of bounded measurable functions.

An alternative model trying to take advantage of total variation and higher-order total variation is a decomposition into two or more parts, that is, $u=u_{1}+u_{2}$ with $u_{1}$ and $u_{2}$ being regularized differently. This was proposed for the first time by Chambolle and Lions (Reference Chambolle and Lions1997) as an infimal convolution of first- and second-order total variations; the effective regularization functional is given by

$$\begin{eqnarray}J(u)=\inf _{u_{1}+u_{2}=u}(|u_{1}|_{BV}+|\unicode[STIX]{x1D6FB}u_{2}|_{BV}).\end{eqnarray}$$

The TGV-type models, proposed by Bredies, Kunisch and Pock (Reference Bredies, Kunisch and Pock2010), became a popular alternative. Instead of decomposing $u$ , these instead decompose the gradient measure $Du$ into $Du_{1}$ and some vector field $u_{2}$ . One version of the regularization functional is then given by

$$\begin{eqnarray}J(u)=\inf _{Du_{1}+u_{2}=Du}(|u_{1}|_{BV}+|u_{2}|_{BV}).\end{eqnarray}$$

The fact that the higher-order part is an arbitrary vector field provides additional freedom that can be beneficial compared to the infimal convolution model (Bredies and Holler Reference Bredies and Holler2015a , Benning, Brune, Burger and Müller Reference Benning, Brune, Burger and Müller2013, Grah Reference Grah2017). We also mention that the original TGV model in Bredies et al. (Reference Bredies, Kunisch and Pock2010) does not use a bounded variation model for $u_{2}$ , but only bounded deformations, that is, the symmetric part of the gradient. Moreover, the approach can be formulated for arbitrary order of regularization. In the dual formulation (3.12) approaches such as infimal convolution or TGV still lead to ${\mathcal{C}}$ being a subset of the unit ball in $L^{\infty }(\unicode[STIX]{x1D6FA})$ , which implies

$$\begin{eqnarray}J(u)\leq |u|_{BV}\quad \text{for all}~u\in BV(\unicode[STIX]{x1D6FA}).\end{eqnarray}$$

On the other hand, for many of them a lower bound inequality can be shown at least when excluding a low- (finite-) dimensional null space (Benning et al. Reference Benning, Brune, Burger and Müller2013), that is, there exists a positive constant $c$ and some linear functionals $\ell _{i}$ such that

$$\begin{eqnarray}J(u)\geq c|u|_{BV}\quad \text{for all}~u\in BV(\unicode[STIX]{x1D6FA}),\quad \text{such that }\ell _{i}(u)=0,i=1,\ldots M.\end{eqnarray}$$

Hence, $J$ is an equivalent norm on the subspace of $BV$ excluding the null space. For the combination of first- and second-order derivatives the null space naturally consists of piecewise affine functions (thus $M=d+1$ ). For further discussion and advanced aspects we refer to Bredies et al. (Reference Bredies, Kunisch and Pock2010), Benning et al. (Reference Benning, Brune, Burger and Müller2013), Ranftl, Pock and Bischof (Reference Ranftl, Pock, Bischof and Kuijper2013), Bredies and Holler (Reference Bredies and Holler2014, Reference Bredies and Holler2015a , Reference Bredies and Holler2015b ), Burger, Papafitsoros, Papoutsellis and Schönlieb (Reference Burger, Papafitsoros, Papoutsellis, Schönlieb and Bociu2015b , Reference Burger, Papafitsoros, Papoutsellis and Schönlieb2016c ), Bergounioux (Reference Bergounioux2016), Setzer, Steidl and Teuber (Reference Setzer, Steidl and Teuber2011), Holler and Kunisch (Reference Holler and Kunisch2014), Gao and Bredies (Reference Gao and Bredies2017) and Bergounioux and Papoutsellis (Reference Bergounioux and Papoutsellis2018).

In certain cases it is also interesting to use total variation regularization on some transform of the image. Motivated by research in image analysis taking into account orientations via local Radon transforms (Krause, Alles, Burgeth and Weickert Reference Krause, Alles, Burgeth and Weickert2016), Burger, Müller, Papoutsellis and Schönlieb (Reference Burger, Müller, Papoutsellis and Schönlieb2014) investigated total variation regularization on the Radon transform, combined with total variation on the image itself, to promote piecewise constant images with very thin structures resembling lines. Grah (Reference Grah2017) investigated total variation on the spherical Radon transform (equivalent to circular Hough transform in computer vision) in order to reconstruct small circular structures.

Another variant is total variation regularization for vector fields, for example arising for colour images (Bresson and Chan Reference Bresson and Chan2008, Blomgren and Chan Reference Blomgren and Chan1998), flow fields (Hinterberger, Scherzer, Schnörr and Weickert Reference Hinterberger, Scherzer, Schnörr and Weickert2002, Zach, Pock and Bischof Reference Zach, Pock, Bischof and Hamprecht2007) or joint reconstruction problems (Knoll et al. Reference Knoll, Holler, Koesters, Otazo, Bredies and Sodickson2017). While many aspects remain the same as in the scalar case, it is particularly interesting to consider which matrix norm should be used for $Du$ , and which dual norm for $g$ (observing that this becomes a matrix in (3.12)).

3.2 Sparsity regularization

Total variation regularization, in particular its discrete version, can be interpreted as a functional favouring sparsity, in this case of the gradient. The paradigm of sparsity has developed in parallel to the total variation regularization (Donoho Reference Donoho1992). A key insight driving sparsity priors was the (approximate) sparsity of signals and natural images in wavelet bases (Mallat and Zhang Reference Mallat and Zhang1993, Huang and Mumford Reference Huang and Mumford1999, Mallat Reference Mallat2008, Starck, Murtagh and Fadili Reference Starck, Murtagh and Fadili2010). Further improvements were made by replacing the orthonormal bases with frames (Christensen Reference Christensen2003) such as curvelets (Candès and Donoho Reference Candès and Donoho2000a , Reference Candès and Donoho2000b ) or shearlets (Labate, Lim, Kutyniok and Weiss Reference Labate, Lim, Kutyniok and Weiss2005, Guo and Labate Reference Guo and Labate2007, Kutyniok and Labate Reference Kutyniok and Labate2012).

Sparsity is naturally measured by the $\ell ^{0}$ -norm,Footnote ¹ the number of non-zero entries. Since the minimization of the $\ell ^{0}$ -norm is highly non-convex and even NP-complete, it is usually relaxed to the convex $\ell ^{1}$ -norm. In the analysis formulation a frame system $\unicode[STIX]{x1D719}_{i}$ is used to test sparsity of $\langle u,\unicode[STIX]{x1D719}_{i}\rangle$ , and the corresponding regularization functional is given by

(3.14) $$\begin{eqnarray}J(u)=\mathop{\sum }_{i}|\langle u,\unicode[STIX]{x1D719}_{i}\rangle |.\end{eqnarray}$$

If $(\unicode[STIX]{x1D719}_{i})$ is an orthonormal system, this is equivalent to the synthesis formulation, which is based on writing

(3.15) $$\begin{eqnarray}J(u)=\mathop{\sum }_{i}|c_{i}|\quad \text{where }u=\mathop{\sum }_{i}c_{i}\unicode[STIX]{x1D719}_{i}.\end{eqnarray}$$

Note that in general the two formulations may differ for frames (Elad, Milanfar and Rubinstein Reference Elad, Milanfar and Rubinstein2007).

In the synthesis formulation we can effectively define the variational problem on the coefficient vector $c$ , that is,

$$\begin{eqnarray}\tilde{K}:\ell ^{2}(\mathbb{N})\rightarrow {\mathcal{V}},\quad c\mapsto \mathop{\sum }_{i}c_{i}K\unicode[STIX]{x1D719}_{i},\end{eqnarray}$$

and compute

$$\begin{eqnarray}\hat{u} =\mathop{\sum }_{i}{\hat{c}}_{i}\unicode[STIX]{x1D719}_{i},\quad c_{i}\in \underset{c}{\operatorname{argmin}}\,F(\tilde{K}c,f)+\unicode[STIX]{x1D6FC}|c|_{1}.\end{eqnarray}$$

The corresponding optimality condition is given by

$$\begin{eqnarray}(\tilde{K}^{\ast }\unicode[STIX]{x2202}_{x}F(\tilde{K}c,f))_{i}+\unicode[STIX]{x1D6FC}s_{i}=0,\end{eqnarray}$$

where $s_{i}$ is a multivalued sign of $c_{i}$ , that is, an element of $[-1,1]$ for $c_{i}=0$ . If $\tilde{K}$ is a bounded linear operator on $\ell ^{2}(\mathbb{N})$ , then its adjoint maps into the same space, and hence $(s_{i})\in \ell ^{2}(\mathbb{N})$ . This implies in particular that $|s_{i}|<1$ , hence $c_{i}=0$ , for $i$ sufficiently large. Thus, we always obtain some sparsity with this model.

In the analysis formulation the optimality condition is instead given by

$$\begin{eqnarray}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku,f)+\unicode[STIX]{x1D6FC}s_{i}\unicode[STIX]{x1D719}_{i}=0,\end{eqnarray}$$

where $s_{i}$ is a multivalued sign for $\langle u,\unicode[STIX]{x1D719}_{i}\rangle$ . Here the understanding of the sparsity property is more complicated: the $s_{i}$ are actually related to the residual via the linear system

$$\begin{eqnarray}\mathop{\sum }_{j}\langle \unicode[STIX]{x1D719}_{i},\unicode[STIX]{x1D719}_{j}\rangle s_{j}=-\displaystyle \frac{1}{\unicode[STIX]{x1D6FC}}\langle K\unicode[STIX]{x1D719}_{i},\unicode[STIX]{x2202}F(Ku,f)\rangle .\end{eqnarray}$$

We refer to Vaiter et al. (Reference Vaiter, Deledalle, Peyré, Dossal and Fadili2013a , Reference Vaiter, Peyré, Dossal and Fadili2013b ) for a detailed analysis in this case. Sparsity models for inverse problems have been studied with different frames and applications extensively in the past decade (e.g. Cotter, Rao, Engan and Kreutz-Delgado Reference Cotter, Rao, Engan and Kreutz-Delgado2005, Chaux, Combettes, Pesquet and Wajs Reference Chaux, Combettes, Pesquet and Wajs2007, Colonna, Easley, Guo and Labate Reference Colonna, Easley, Guo and Labate2010, Recht, Fazel and Parrilo Reference Recht, Fazel and Parrilo2010).

There are several relevant extensions of sparsity priors to multidimensional systems, in particular via a synthesis-type formulation:

$$\begin{eqnarray}u=\mathop{\sum }_{i,j}c_{ij}\unicode[STIX]{x1D719}_{i}\otimes \unicode[STIX]{x1D713}_{j}.\end{eqnarray}$$

The different dimensions are often space (characterized by basis functions $\unicode[STIX]{x1D719}_{i}$ ) and time or frequency (characterized by basis functions $\unicode[STIX]{x1D713}_{j}$ ). Instead of overall sparsity, more detailed prior knowledge can be introduced. The most popular example is joint or collaborative sparsity, which means that only a few of the basis functions, for example in the second dimension, can be used to explain the solution. This means that $c_{\cdot j}$ vanishes for most $j$ , as does the norm of $c_{\cdot j}$ . A common regularization for this case is the joint or collaborative sparsity prior

$$\begin{eqnarray}J(u)=\mathop{\sum }_{j}\Vert c_{\cdot j}\Vert _{\ell ^{r}},\end{eqnarray}$$

usually with $r=2$ or $r=\infty$ (Duarte et al. Reference Duarte, Sarvotham, Wakin, Baron and Baraniuk2005, Teschke and Ramlau Reference Teschke and Ramlau2007, Fornasier and Rauhut Reference Fornasier and Rauhut2008, Gholami and Siahkoohi Reference Gholami and Siahkoohi2010, Lee, Kim, Bresler and Ye Reference Lee, Kim, Bresler and Ye2011). An alternative type of prior knowledge is local sparsity, which means that for each $i$ only a few basis functions $\unicode[STIX]{x1D713}_{j}$ are used. The term local is due to an imaging interpretation of the $\unicode[STIX]{x1D719}_{i}$ as basis functions local in space (e.g. for each pixel). This is a common issue in dynamic or spectral imaging, where one can assume that only a few materials and their characteristic evolutions or spectral curves can be found in each pixel. A regularization functional proposed for this issue (Heins, Moeller and Burger Reference Heins, Moeller and Burger2015) is

$$\begin{eqnarray}J(u)=\max _{i}\Vert c_{i\cdot }\Vert _{\ell ^{1}}+\unicode[STIX]{x1D6FD}\mathop{\sum }_{i}\Vert c_{i\cdot }\Vert _{\ell ^{1}}.\end{eqnarray}$$

An infinite-dimensional extension of the above sparsity models is sparsity in a space of Radon measures, that is, the regularization functional is given as the total variation norm of the measure $u$ :

$$\begin{eqnarray}J(u)=\int _{\unicode[STIX]{x1D6FA}}\,\text{d}|u|=\sup _{g\in C_{0}(\unicode[STIX]{x1D6FA}),\Vert g\Vert _{\infty }\leq 1}\int _{\unicode[STIX]{x1D6FA}}g\,\text{d}u.\end{eqnarray}$$

This yields a convex regularization functional for reconstructing multiple peaks at unknown locations, and was proposed for inverse problems in Bredies and Pikkarainen (Reference Bredies and Pikkarainen2013) and for super-resolution problems in Candès and Fernandez-Granda (Reference Candès and Fernandez-Granda2013, Reference Candès and Fernandez-Granda2014) and Aja-Fernandez, Alberola-Lopez and Westin (Reference Aja-Fernandez, Alberola-Lopez and Westin2008). The reconstruction properties in deconvolution problems were analysed in Duval and Peyré (Reference Duval and Peyré2017a , Reference Duval and Peyré2017b ) and Denoyelle, Duval and Peyré (Reference Denoyelle, Duval and Peyré2017), and asymptotics from finite-dimensional problems with sparsity priors are found in Heins (Reference Heins2014) and Duval and Peyré (Reference Duval and Peyré2017a , Reference Duval and Peyré2017b ).

3.3 Low-rank regularization

In many applications one seeks a decomposition of the form

(3.16) $$\begin{eqnarray}U=\mathop{\sum }_{i}\unicode[STIX]{x1D6F7}_{i}\otimes \unicode[STIX]{x1D6F9}_{i}\end{eqnarray}$$

with unknown $\unicode[STIX]{x1D6F7}_{i},\unicode[STIX]{x1D6F9}_{i}$ and the additional prior knowledge that there are as few elements as possible in the sum. In a finite-dimensional setting this means that the matrix $U$ has low rank, that is, the rank of $U$ would be the obvious regularization functional. However, since the rank is very far from being convex, several relaxations have been proposed instead. The most popular one, originally proposed for matrix completion problems, is the nuclear norm (Candès and Recht Reference Candès and Recht2009, Recht, Fazel and Parrilo Reference Recht, Fazel and Parrilo2010, Cai, Candès and Shen Reference Cai, Candès and Shen2010, Cai and Osher Reference Cai and Osher2013, Yang, Ma and Osher Reference Yang, Ma and Osher2013)

(3.17) $$\begin{eqnarray}\Vert U\Vert _{\ast }=\sum \unicode[STIX]{x1D70E}_{i},\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{i}$ are the singular values of $U$ .

In many applications the low-rank part alone does not suffice to model the structure of solutions, frequently a low-rank plus sparsity (L+S) model is employed instead (Candès, Li, Ma and Wright Reference Candès, Li, Ma and Wright2011, Otazo, Candès and Sodickson Reference Otazo, Candès and Sodickson2015), which is again based on the decomposition

(3.18) $$\begin{eqnarray}J(u)=\inf _{u_{1}+u_{2}=u}(\Vert u_{1}\Vert _{\ast }+\Vert Tu_{2}\Vert _{1}),\end{eqnarray}$$

with a sparsifying transform $T$ (often some derivative as in total variation). Particularly in videos the low-rank part captures background and certain slow dynamics, while the sparse part captures the key changes.

For inverse problems an infinite-dimensional function space setting would be more appropriate, which has not yet been investigated. In particular, a formulation in a space of trace class operators between Hilbert spaces $H_{1}$ and $H_{2}$ (Reed and Simon Reference Reed and Simon1978) would be natural. Let us mention that the choice of Hilbert spaces $H_{i}$ opens novel opportunities for improved regularization that have not been exploited as yet, even in the finite-dimensional case.

3.4 Infimal convolutions

As we have seen above, infimal convolution is a versatile tool to combine different regularization approaches, and to define a novel functional that combines their advantages. We want to emphasize this approach in the following by providing formal definitions.

Definition 3.1. Let $J_{i}:{\mathcal{U}}\rightarrow \mathbb{R}\cup \{+\infty \}$ , $i=1,2$ be proper convex functionals. Then their infimal convolution $J_{1}\Box J_{2}:{\mathcal{U}}\rightarrow \mathbb{R}\cup \{+\infty \}$ is defined by

(3.19) $$\begin{eqnarray}(J_{1}\Box J_{2})(u)=\inf _{v\in {\mathcal{U}}}(J_{1}(u-v)+J_{2}(v)).\end{eqnarray}$$

More generally, we can define an infimal convolution for an arbitrary number of convex functionals.

Definition 3.2. Let $J_{i}:{\mathcal{U}}\rightarrow \mathbb{R}\cup \{+\infty \}$ , $i=1,\ldots ,M$ be proper convex functionals. Then their infimal convolution $J:{\mathcal{U}}\rightarrow \mathbb{R}\cup \{+\infty \}$ is defined by

(3.20) $$\begin{eqnarray}J(u)=\inf _{u_{i}\in {\mathcal{U}},\sum u_{i}=u}\mathop{\sum }_{i=1}^{M}J_{i}(u_{i})\end{eqnarray}$$

We mention that a priori it is unclear whether the infima above are actually minima. If a minimizer $v$ exists for the infimal convolution of $J_{1}$ and $J_{2}$ , it can be used to deduce optimality conditions, since

$$\begin{eqnarray}p\in \unicode[STIX]{x2202}J(u)\quad \text{if }p\in \unicode[STIX]{x2202}J_{1}(u-v)\cap \unicode[STIX]{x2202}J_{2}(v).\end{eqnarray}$$

As the above examples for sparsity and in particular higher-order total variation show, there is great freedom in designing infimal convolution models for regularization. Consequently, a lot of options for future research remain open, and interesting results are still to be expected.

3.5 Regularization by denoising

Romano, Elad and Milanfar (Reference Romano, Elad and Milanfar2017) recently proposed an approach named regularization by denoising, where $J$ is of the form

$$\begin{eqnarray}J(u)=\langle u,u-D(u)\rangle .\end{eqnarray}$$

Here $D$ is a potentially nonlinear operator that maps an image $u$ to the output of a black-box image denoising routine. The idea is to make state-of-the-art denoising routines available for the regularization of linear inverse problems.

3.6 Bregman distances

From a single regularization functional several variants can be constructed by using a non-trivial prior $u_{0}$ and the so-called Bregman distance (originally introduced in Bregman (Reference Bregman1967) for proximal-point-type methods). Instead of shifting the functional directly from $J(u)$ to $J(u-u_{0})$ , the approach in the Bregman distance performs a shift in the convex conjugate. In the original formulation this amounts to the following.

Definition 3.3. Let $J:{\mathcal{U}}\rightarrow \mathbb{R}\cup \{+\infty \}$ be a convex functional and let $p_{0}\in \unicode[STIX]{x2202}J(u_{0})$ . Then the Bregman distance between $u\in {\mathcal{U}}$ and $u_{0}\in {\mathcal{U}}$ with subgradient $p_{0}$ is given by

(3.21) $$\begin{eqnarray}D_{J}^{p_{0}}(u,u_{0}):=J(u)-J(u_{0})-\langle p_{0},u-u_{0}\rangle\end{eqnarray}$$

Note that the Bregman distance is not a strict distance, that is, it can vanish for $u\neq u_{0}$ if $J$ is not strictly convex. It is also not symmetric, but can be made symmetric by taking a sum of one-sided distances (see Burger Reference Burger and Hiriart-Urruty2016 for a more detailed discussion). For absolutely one-homogeneous regularization functionals as above, the identity $J(u_{0})=\langle p_{0},u_{0}\rangle$ holds, so the Bregman distance becomes

(3.22) $$\begin{eqnarray}D_{J}^{p_{0}}(u,u_{0}):=J(u)-\langle p_{0},u\rangle ,\end{eqnarray}$$

and thus it is effectively independent of $u_{0}$ : only the subgradient $p_{0}$ matters. This is particularly relevant if the subdifferential of $J$ is not a singleton, or vice versa a subgradient $p_{0}\in \unicode[STIX]{x2202}J(u_{0})$ can be an element of the subdifferential also at other values of $u$ .

Note that in the typical case of $u_{0}=0$ being a minimizer of $J$ , that is, $0\in \unicode[STIX]{x2202}J(0)$ , the regularization with $J$ can be reinterpreted as penalizing the Bregman distance to $u_{0}=0$ . Bleyer and Leitao (Reference Bleyer and Leitao2009) carried out a basic analysis of such a variational regularization. The topic has received recent interest, particularly in the context of TV-type regularization in imaging, since it allows us to introduce structural information. The key insight in total variation is that the subgradient encodes information about the discontinuity set, more precisely $p=\unicode[STIX]{x1D6FB}\cdot g$ , with $g$ being equal to the unit normal vector to the discontinuity set where it is regular. This is again related to (3.22); the total variation does not depend directly on $u_{0}$ and in particular the contrast in the image. Instead, it vanishes for all $u$ of the form

$$\begin{eqnarray}u(x)=f(u_{0}(x))\end{eqnarray}$$

with a monotonically increasing function $f$ , that is, a simple contrast change (Resmerita and Scherzer Reference Resmerita and Scherzer2006). Assuming that $g$ is a vector field realizing the supremum in the dual definition of the total variation, the Bregman distance becomes

$$\begin{eqnarray}D_{TV}^{p_{0}}(u,u_{0})=|u|_{BV}-\int _{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D6FB}\cdot g_{0})u\,\text{d}x=\int _{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D6FB}\cdot (g-g_{0}))u\,\text{d}x,\end{eqnarray}$$

and if $u$ is piecewise constant with regular discontinuity set $S_{u}$ ,

$$\begin{eqnarray}D_{TV}^{p_{0}}(u,u_{0})=\int _{S_{u}}[u](g-g_{0})\cdot \unicode[STIX]{x1D708}\,\text{d}\unicode[STIX]{x1D70E}=\int _{S_{u}}[u](1-g_{0}\cdot \unicode[STIX]{x1D708})\,\text{d}\unicode[STIX]{x1D70E},\end{eqnarray}$$

where $[u]$ denotes the jump along $S_{u}$ and $\unicode[STIX]{x1D708}$ is the unit normal (oriented such that $[u]$ is positive). One thus observes that the Bregman distance measures differences in the discontinuity set and its orientation, which is perfect for imaging applications with a structural prior (Kaipio, Kolehmainen, Vauhkonen and Somersalo Reference Kaipio, Kolehmainen, Vauhkonen and Somersalo1999) that mainly yields information about edges, i.e. discontinuity sets. One suitable example is that of anatomical priors in medical imaging, where a high-resolution method such as CT or MR is used to obtain information about organ boundaries and other anatomical features, which are the natural candidates for edge sets in functional methods such as PET, SPECT or MR imaging with special contrast. In some cases a joint reconstruction is also of interest, the most obvious case being colour or hyperspectral images, where naturally intensity changes at the same locations, usually even in the same direction (Moeller, Wittman, Bertozzi and Burger Reference Moeller, Wittman, Bertozzi and Burger2012, Moeller, Brinkmann, Burger and Seybold Reference Moeller, Brinkmann, Burger and Seybold2014).

In some applications one may find contrast inversion, that is, the jump of the two images along the discontinuity set has a different sign. In such cases the normals are parallel, which means they point in opposite directions and hence lead to large values in the Bregman distance. A potential solution for avoiding such issues is the infimal convolution of Bregman distances, in this case with the two normal fields and thus subgradients of opposite sign (Moeller et al. Reference Moeller, Brinkmann, Burger and Seybold2014, Rasch et al. Reference Rasch, Kolehmainen, Nivajärvi, Kettunen, Gröhn, Burger and Brinkmann2017):

$$\begin{eqnarray}J=D_{TV}^{p_{0}}(\cdot ,u_{0})\Box D_{TV}^{-p_{0}}(\cdot ,-u_{0}).\end{eqnarray}$$

We also mention some other related approaches to modifying total variation functionals, such as parallel level set models (Ehrhardt and Arridge Reference Ehrhardt and Arridge2014, Ehrhardt et al. Reference Ehrhardt, Thielemans, Pizarro, Atkinson, Ourselin, Hutton and Arridge2014, Ehrhardt et al. Reference Ehrhardt, Markiewicz, Liljeroth, Barnes, Kolehmainen, Duncan, Pizarro, Atkinson, Hutton and Ourselin2016), which can be related to the Bregman distance for total variation (Rasch, Brinkmann and Burger Reference Rasch, Brinkmann and Burger2018), or directional / structural total variation (Bungert et al. Reference Bungert, Coomes, Ehrhardt, Rasch, Reisenhofer and Schönlieb2018, Ehrhardt and Betcke Reference Ehrhardt and Betcke2016, Hintermüller, Holler and Papafitsoros Reference Hintermüller, Holler and Papafitsoros2017, Grasmair and Lenzen Reference Grasmair and Lenzen2010), formally

$$\begin{eqnarray}TV_{g_{0}}(u)=\int _{\unicode[STIX]{x1D6FA}}|(I-g_{0}\otimes g_{0})\unicode[STIX]{x1D6FB}u|\,\text{d}x.\end{eqnarray}$$

4.1 Abstract linear regularization methods

We start our exposition with a discussion of possible limits of regularization schemes. In essentially all linear methods such as Tikhonov regularization, truncated SVD or iterative regularization in Hilbert spaces, it is clear which solutions are approximated as the regularization parameter tends to zero, namely those obtained from a generalized inverse. The following definitions are made to characterize these limiting solutions.

Definition 4.1. Let $K:{\mathcal{U}}\rightarrow {\mathcal{V}}$ be a bounded linear operator between Hilbert spaces and $f\in {\mathcal{V}}$ . We call $\hat{u} \in {\mathcal{U}}$ a best approximate solution of (1.1) if

(4.1) $$\begin{eqnarray}\Vert K\hat{u} -f\Vert _{{\mathcal{V}}}\leq \Vert Ku-f\Vert _{{\mathcal{V}}}\quad \text{for all}~u\in {\mathcal{U}}.\end{eqnarray}$$

Moreover, we call $\hat{u}$ a minimal norm solution if it is a best approximate solution and

(4.2) $$\begin{eqnarray}\Vert \hat{u} \Vert _{{\mathcal{U}}}\leq \Vert u\Vert _{{\mathcal{U}}}\quad \text{for all}~u\in {\mathcal{U}},~\Vert K\hat{u} -f\Vert _{{\mathcal{V}}}=\Vert Ku-f\Vert _{{\mathcal{V}}}.\end{eqnarray}$$

Note that due to the strict convexity of the square of a Hilbert space norm, the minimum solution – being its minimizer on a linear manifold – is a unique object. An abstract regularization method is now a collection of continuous operators approximating the (discontinuous) generalized inverse of $K$ .

Definition 4.2. The bounded linear operators $R_{\unicode[STIX]{x1D6FC}}:{\mathcal{V}}\rightarrow {\mathcal{U}}$ defined for $\unicode[STIX]{x1D6FC}$ in $(0,\unicode[STIX]{x1D6FC}_{0})$ are called linear regularization operators. Together with a parameter choice strategy $\unicode[STIX]{x1D6FC}$ depending on the noise level $\unicode[STIX]{x1D6FF}$ and the data $f^{\unicode[STIX]{x1D6FF}}$ , that is, a function

(4.3) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}:(0,\unicode[STIX]{x1D6FF}_{0})\times {\mathcal{V}}\rightarrow (0,\unicode[STIX]{x1D6FC}_{0}),\end{eqnarray}$$

it is called the linear regularization method.

A linear regularization method is called convergent if, for all $f\in {\mathcal{R}}(K)$ , the condition

(4.4) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D6FF}\rightarrow 0}\sup \{\Vert R_{\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},f^{\unicode[STIX]{x1D6FF}})}(f^{\unicode[STIX]{x1D6FF}})-u^{\dagger }\Vert _{{\mathcal{U}}}\mid f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}},\Vert f-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{V}}}\leq \unicode[STIX]{x1D6FF}\}=0\end{eqnarray}$$

holds with $u^{\dagger }$ being the minimum norm solution of (1.1).

For ill-posed problems it is well known that convergence can be arbitrarily slow (Schock Reference Schock, Hämmerlin and Hoffmann1985). Thus convergence rates can be obtained only on a restricted subset $M_{\unicode[STIX]{x1D708}}$ with the parameter $\unicode[STIX]{x1D708}>0$ measuring the smoothness or order of convergence. The standard definition is as follows.

Definition 4.3. A regularization method is convergent at order $\unicode[STIX]{x1D708}$ on a set ${\mathcal{M}}_{\unicode[STIX]{x1D708}}$ if, given any $f=Ku^{\dagger }$ , for $u^{\dagger }\in {\mathcal{M}}_{\unicode[STIX]{x1D708}}$ , there exists a constant $C_{\unicode[STIX]{x1D708}}$ such that

(4.5) $$\begin{eqnarray}\Vert R_{\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF},f^{\unicode[STIX]{x1D6FF}})}(f^{\unicode[STIX]{x1D6FF}})-u^{\dagger }\Vert \leq C_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D6FF}^{\unicode[STIX]{x1D708}}\end{eqnarray}$$

for any data $f^{\unicode[STIX]{x1D6FF}}$ satisfying $\Vert f^{\unicode[STIX]{x1D6FF}}-f\Vert \leq \unicode[STIX]{x1D6FF}$ .

It is well known that the set $M_{\unicode[STIX]{x1D708}}$ can be related to the source condition

$$\begin{eqnarray}u^{\dagger }=(K^{\ast }K)^{\unicode[STIX]{x1D707}}w\end{eqnarray}$$

for some $w\in {\mathcal{U}}$ and appropriate $\unicode[STIX]{x1D707}>0$ related to $\unicode[STIX]{x1D708}$ (Engl et al. Reference Engl, Hanke and Neubauer1996). The constant $C_{\unicode[STIX]{x1D708}}$ is then related to the norm of $w$ . The simplest cases of source conditions are $\unicode[STIX]{x1D707}=1/2$ , which can be reformulated as

$$\begin{eqnarray}u^{\dagger }=K^{\ast }\tilde{w},\end{eqnarray}$$

for some $\tilde{w}\in {\mathcal{V}}$ , and $\unicode[STIX]{x1D707}=1$ . Source conditions induce conditional well-posedness of the problem, e.g. for $\unicode[STIX]{x1D707}=1/2$ , one has for $u_{i}=K^{\ast }\tilde{w}_{i}$

$$\begin{eqnarray}\Vert u_{1}-u_{2}\Vert ^{2}=\langle u_{1}-u_{2},K^{\ast }(\tilde{w}_{1}-\tilde{w}_{2})\rangle =\langle K(u_{1}-u_{2}),\tilde{w}_{1}-\tilde{w}_{2}\rangle .\end{eqnarray}$$

The Cauchy–Schwarz and triangle inequality then imply the Hölder stability

$$\begin{eqnarray}\Vert u_{1}-u_{2}\Vert \leq C\sqrt{\Vert Ku_{1}-Ku_{2}\Vert },\end{eqnarray}$$

with $C=\sqrt{\Vert \tilde{w}_{1}\Vert +\Vert \tilde{w}_{2}\Vert }.$

4.2 Extension to nonlinear methods

The examples of variational regularization models in the previous section call for a more general theory of nonlinear regularization methods. While the concept of a best-approximate solution is fairly straightforward to generalize, other aspects of convergence and limiting solutions are less obvious. In a general variational regularization, as in the examples discussed above, it would be natural to replace the minimum norm solution with a solution minimizing the regularization functional. The latter is not necessarily unique, however, hence some possible multivaluedness needs to be introduced in the characterization. Similar issues apply to the regularized problem and hence the definition of a regularization operator. In the following we will try to provide a fundamental setting for nonlinear regularization methods. As in the case of linear regularizations we first generalize the possible types of solutions we would like to approximate. The generalization of the first notion is fairly straightforward, so we only allow for more general distance measures, for example functionals related to negative log-likelihoods for non-Gaussian distributions.

Definition 4.4. Given an error measure $F:{\mathcal{V}}\times {\mathcal{V}}\rightarrow \mathbb{R}_{+}\cup \{+\infty \}$ , we call $\hat{u} \in {\mathcal{U}}$ a best approximate solution of (1.1) with respect to $F$ if

(4.6) $$\begin{eqnarray}F(K\hat{u} ,f)\leq F(Ku,f)\quad \text{for all}~u\in {\mathcal{U}}.\end{eqnarray}$$

A suitable generalization of the definition of a minimum norm solution is more involved. In particular, we would like to give a unified concept including the selection by minimizing a regularization functional or maximizing some prior probability. We encode the selection of specific solutions due to prior knowledge in a (multivalued) selection operator.

The general set-valued definition of a selection operator, which we use in order to take care of all the possible cases in regularization methods, also needs to use set-valued ways of convergence. Hence we recall the definition of Kuratowski convergence in a metric space.

Definition 4.6. Given a metric space $X$ with metric $d$ , we abuse notation by defining

(4.7) $$\begin{eqnarray}d(u,S):=\inf _{v\in S}d(u,v),\end{eqnarray}$$

for $x\in X$ and $S\subset X$ . Then the Kuratowski limit inferior and the Kuratowski limit superior of a sequence of sets $S_{n}\subset X$ are defined as follows:

(4.8) $$\begin{eqnarray}\displaystyle K-\liminf _{n}(S_{n}) & = & \displaystyle \{x\in X\mid \limsup _{n}d(x,S_{n})=0\},\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle K-\limsup _{n}(S_{n}) & = & \displaystyle \{x\in X\mid \liminf _{n}d(x,S_{n})=0\}.\end{eqnarray}$$

For our purposes the lim sup will be of particular interest, because we will use a minimal definition of stability often adopted in the literature on nonlinear methods following Seidman and Vogel (Reference Seidman and Vogel1989) and Engl, Kunisch and Neubauer (Reference Engl, Kunisch and Neubauer1989). Stability is expressed by subsequences of selected solutions having a limit, and each limit of a subsequence being a solution of the limiting problem. The lim inf is less interesting, since there is no reason to ask that any solution of a problem can be the limit of approximate problems. We call an inverse problem stable if, for $f_{n}\rightarrow f$ (usually in terms of norm convergence in ${\mathcal{V}}$ ), we have

(4.10) $$\begin{eqnarray}K-\limsup _{n}{\mathcal{S}}(f_{n})\subset {\mathcal{S}}(f),\quad K-\limsup _{n}{\mathcal{S}}(f_{n})\neq \emptyset .\end{eqnarray}$$

The metric used for the Kuratowski lim sup will usually be a metrization of some weak or even weak-star convergence in a Banach space; one might also use an extension of the definition to other distance measures.

Having defined the solutions we would like to approximate, the obvious next step is to define a (convergent) regularization method. We start in a deterministic setting, generalizing to a vectorial regularization parameter $\unicode[STIX]{x1D736}\in \mathbb{R}_{+}^{M}$ , however, which is useful in many examples, for example the TGV and infimal convolution models with multiple parameters mentioned above. Given an error measure $F$ and $f=Ku^{\dagger }$ for some exact solution $u^{\dagger }\in {\mathcal{U}}$ , we call $\unicode[STIX]{x1D6FF}>0$ the noise level if it is the best available bound for available data $f^{\unicode[STIX]{x1D6FF}}$ , that is,

(4.11) $$\begin{eqnarray}F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}.\end{eqnarray}$$

We will be interested in the convergence of regularized solutions to prior selected solutions as the noise level tends to zero. For ease of presentation and since this is available in almost any known example, we restrict ourselves to convergence with respect to a metric topology $\unicode[STIX]{x1D70F}$ , which is usually a weak or weak-star topology (on some bounded set in the Banach space).

Definition 4.7. The multivalued operators $R(\cdot ,\unicode[STIX]{x1D736}):{\mathcal{V}}\rightrightarrows {\mathcal{U}}$ defined for $\unicode[STIX]{x1D736}$ in a subset $A$ of $\mathbb{R}^{M}$ are called regularization operators if, for each $\unicode[STIX]{x1D736}\in A$ , the operator $R$ satisfies the stability property

(4.12) $$\begin{eqnarray}\emptyset \neq K-\liminf _{n}R(f^{\unicode[STIX]{x1D6FF}_{n}},\unicode[STIX]{x1D736})\subset R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736})\end{eqnarray}$$

for all $f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}}$ and sequences $f^{\unicode[STIX]{x1D6FF}_{n}}\in {\mathcal{V}}$ converging to $f^{\unicode[STIX]{x1D6FF}}$ . Together with a parameter choice strategy $\unicode[STIX]{x1D736}$ depending on the noise level $\unicode[STIX]{x1D6FF}$ and the data $f^{\unicode[STIX]{x1D6FF}}$ , that is, a function

(4.13) $$\begin{eqnarray}\unicode[STIX]{x1D736}:(0,\unicode[STIX]{x1D6FF}_{0})\times {\mathcal{V}}\rightarrow A,\end{eqnarray}$$

it is called a regularization method.

A regularization method is called convergent if, for any sequence $\unicode[STIX]{x1D6FF}_{n}\rightarrow 0$ and data $f^{\unicode[STIX]{x1D6FF}_{n}}$ satisfying

(4.14) $$\begin{eqnarray}F(f,f^{\unicode[STIX]{x1D6FF}_{n}})\leq \unicode[STIX]{x1D6FF}_{n},\end{eqnarray}$$

we have

(4.15) $$\begin{eqnarray}\emptyset \neq K-\liminf _{n}R(f^{\unicode[STIX]{x1D6FF}_{n}},\unicode[STIX]{x1D736}(\unicode[STIX]{x1D6FF}_{n},f^{\unicode[STIX]{x1D6FF}_{n}}))\subset {\mathcal{S}}(f).\end{eqnarray}$$

We mention that – besides the very general set-up – our definition of a regularization method deviates from the usual theory, since we do not assume any kind of convergence of the regularization parameter $\unicode[STIX]{x1D736}$ . In the classical theory and most examples, $\unicode[STIX]{x1D736}$ is a scalar positive value and assumed to converge to zero (or to infinity) as the noise level tends to zero. However, apart from the convenience there seems to be no reason to put such convergence into the definition. Note that in order to approximate a really ill-posed problem, each clustering point of $\unicode[STIX]{x1D736}(\unicode[STIX]{x1D6FF}_{n},f^{\unicode[STIX]{x1D6FF}_{n}})$ will automatically lie outside $A$ . The canonical examples are $A=(0,\unicode[STIX]{x1D6FC}_{0})$ for variational regularization or $A=\mathbb{N}$ for iterative regularization, where the limiting parameter will converge to zero or infinity. However, we may also consider multi-parameter regularization, where it depends on the formulation whether each component of $\unicode[STIX]{x1D736}$ has a limit outside the admissible set. Take for example an infimal convolution of two functionals $R_{1}$ and $R_{2}$ . If $\unicode[STIX]{x1D736}=(\unicode[STIX]{x1D6FC}_{1},\unicode[STIX]{x1D6FC}_{2})$ are the coefficients of $R_{1}$ and $R_{2}$ , then obviously both should tend to zero in the limit. If, however, $\unicode[STIX]{x1D6FC}_{2}$ is a relative parameter, that is, $\unicode[STIX]{x1D6FC}_{1}$ is the coefficient of $R_{1}$ and $\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D6FC}_{2}$ the coefficient of $R_{2}$ , then it is natural to have a positive limit of $\unicode[STIX]{x1D6FC}_{2}$ . Further motivation for our general definition comes from recent approaches to learning regularization methods for inverse problems, where the $\unicode[STIX]{x1D736}$ can represent the parameters of the learning scheme. To get a consistent infinite-dimensional theory one could even generalize to non-parametric learning, which would amount to choosing $\unicode[STIX]{x1D736}$ in some Banach space. Note that in the remainder of this article we will often write $\unicode[STIX]{x1D6FC}$ instead of $\unicode[STIX]{x1D736}$ if $\unicode[STIX]{x1D736}$ is only a scalar.

In order to define convergence rates we will also need an error measure $D:{\mathcal{U}}\times {\mathcal{U}}\rightarrow \mathbb{R}_{+}\cup \{+\infty \}$ , since there is no natural norm measure as in the Hilbert space. Moreover, we need a restriction to appropriate classes of smoothness, which we denote by $M_{\unicode[STIX]{x1D708}}$ with a parameter $\unicode[STIX]{x1D708}>0$ measuring the smoothness.

Definition 4.8. A regularization method is called $D$ -convergent if

(4.16) $$\begin{eqnarray}\lim _{\unicode[STIX]{x1D6FF}\rightarrow 0}\sup \{D(u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D736}},u^{\dagger })\mid u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D736}}\in R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736}),f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}},F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}\}=0.\end{eqnarray}$$

A regularization method is called convergent at order $\unicode[STIX]{x1D708}$ on a set if, for all $f=Ku^{\dagger }$ , $u^{\dagger }\in {\mathcal{M}}_{\unicode[STIX]{x1D708}}$ , there exists a constant $C_{\unicode[STIX]{x1D708}}$ such that, for all data $g$ satisfying (4.11), we obtain the estimate

(4.17) $$\begin{eqnarray}D(R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736}(\unicode[STIX]{x1D6FF},f^{\unicode[STIX]{x1D6FF}})),u^{\dagger })\leq C_{\unicode[STIX]{x1D708}}\unicode[STIX]{x1D6FF}^{\unicode[STIX]{x1D708}}.\end{eqnarray}$$

Of course the above definition only makes sense for suitable choices of the distance functional and smoothness classes. Remember that in the classical linear Hilbert space theory, these were just norms and spaces obtained by source conditions. We will discuss generalizations of these in the nonlinear setting, in particular related to variational and iterative regularization methods in Banach spaces related to convex regularization functionals. Note also that more general rates than just polynomial ones have been considered in the literature (e.g. Hohage Reference Hohage1997, Kaltenbacher Reference Kaltenbacher2008).

From an abstract point of view the key insight for generalizing source conditions is the range of the regularization operator. For many linear regularization methods in Hilbert spaces, it is easy to see that the source condition $u^{\dagger }=K^{\ast }\tilde{w}$ means that there exist some data $f^{\dagger }$ with $u^{\dagger }=R(f^{\dagger },\unicode[STIX]{x1D736})$ . As examples, take Tikhonov regularization:

$$\begin{eqnarray}R(\cdot ,\unicode[STIX]{x1D736})=(K^{\ast }K+\unicode[STIX]{x1D6FC}I)^{-1}K^{\ast }=K^{\ast }(KK^{\ast }+\unicode[STIX]{x1D6FC}I)^{-1}.\end{eqnarray}$$

Due to the invertibility of $(KK^{\ast }+\unicode[STIX]{x1D6FC}I)^{-1}$ , the range of the regularization operator coincides with the range of $K^{\ast }$ . Instead of defining source conditions at an abstract level, we thus make the following definition.

Definition 4.9 (range condition).

An element $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ for $f\in {\mathcal{R}}(K)$ satisfies the range condition if $u^{\dagger }\in {\mathcal{R}}(R(\cdot ,\unicode[STIX]{x1D736}))$ , that is, there exists $f_{\unicode[STIX]{x1D6FC}}^{\dagger }$ such that

$$\begin{eqnarray}u^{\dagger }\in R(f_{\unicode[STIX]{x1D736}}^{\dagger },\unicode[STIX]{x1D736}).\end{eqnarray}$$

We mention that in the case of nonlinear variational methods (with quadratic fidelity), the equivalence of a nonlinear source condition and the range condition was shown in Burger and Osher (Reference Burger and Osher2004), again confirming the appropriateness of this definition.

Roughly speaking, error estimates can now be obtained by some continuity property of the regularization operator, which implies

$$\begin{eqnarray}d_{{\mathcal{U}}}(u_{\unicode[STIX]{x1D736}}^{\unicode[STIX]{x1D6FF}},u^{\dagger })\leq C(\unicode[STIX]{x1D6FC})d_{{\mathcal{V}}}(f^{\unicode[STIX]{x1D6FF}},f_{\unicode[STIX]{x1D736}}^{\dagger }),\end{eqnarray}$$

with appropriate distances $d_{{\mathcal{U}}}$ and $d_{{\mathcal{V}}}$ . With some kind of triangle inequality the right-hand side can be estimated by a distance between $f$ and $f^{\unicode[STIX]{x1D6FF}}$ , which is related to the noise level as well as a distance between $f$ and $f_{\unicode[STIX]{x1D736}}^{\dagger }$ , which is related to the bias of the regularization. This will be discussed in detail for the case of variational regularization methods in Section 5. A weaker concept is that of approximate source conditions (e.g. Schuster, Kaltenbacher, Hofmann and Kazimierski Reference Schuster, Kaltenbacher, Hofmann and Kazimierski2012, Burger, Helin and Kekkonen Reference Burger, Helin and Kekkonen2016b ) that effectively measure how well the range condition can be approximated. On the other hand stronger conditions can be obtained if $f_{\unicode[STIX]{x1D736}}^{\dagger }$ above is not arbitrary but in the range of the forward operator $K$ .

4.3 Stochastic approaches

In addition to the deterministic viewpoint, a statistical approach has recently become popular in infinite-dimensional problems as well (Bissantz, Hohage and Munk Reference Bissantz, Hohage and Munk2004, Bissantz, Hohage, Munk and Ruymgaart Reference Bissantz, Hohage, Munk and Ruymgaart2007, Cavalier Reference Cavalier2008, Kekkonen, Lassas and Siltanen Reference Kekkonen, Lassas and Siltanen2014, Giné and Nickl Reference Giné and Nickl2015, Hohage and Werner Reference Hohage and Werner2016). In such a set-up the data $f^{\unicode[STIX]{x1D6FF}}$ are considered to be random variables drawn from a measure $\unicode[STIX]{x1D707}_{f}$ centred around the exact data $f$ (often representing the expected value and $\unicode[STIX]{x1D6FF}$ being some kind of variance). A regularization operator can then still be applied to each realization and defined in the same way, but we need a different definition of the noise level and the convergence of the regularization method. As a generalization of variance we use the statistical noise level in the mean

(4.18) $$\begin{eqnarray}\mathbb{E}(F(f,f^{\unicode[STIX]{x1D6FF}}))=\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

Definition 4.10. A regularization operator $R$ with a parameter choice strategy $\unicode[STIX]{x1D736}$ depending on the statistical noise level $\unicode[STIX]{x1D6FF}$ and the data $f^{\unicode[STIX]{x1D6FF}}$ , that is, a function

(4.19) $$\begin{eqnarray}\unicode[STIX]{x1D736}:(0,\unicode[STIX]{x1D6FF}_{0})\times {\mathcal{V}}\rightarrow A,\end{eqnarray}$$

is called a statistical regularization method.

A statistical regularization method is called convergent if, for all sequences $\unicode[STIX]{x1D6FF}_{n}\rightarrow 0$ , random variables $f^{\unicode[STIX]{x1D6FF}_{n}}$ satisfying

(4.20) $$\begin{eqnarray}\mathbb{E}(F(f,f^{\unicode[STIX]{x1D6FF}_{n}}))\leq \unicode[STIX]{x1D6FF}_{n},\end{eqnarray}$$

and each choice of random variables $u_{n}\in R_{\unicode[STIX]{x1D736}_{n}}(f^{\unicode[STIX]{x1D6FF}_{n}})$ , there exists a convergent subsequence $u_{n_{k}}$ in probability in the topology $\unicode[STIX]{x1D70F}$ , and the limiting random variable $u^{\dagger }$ satisfies $u^{\dagger }\in {\mathcal{S}}(f)$ with probability one.

An extension of this viewpoint is the Bayesian approach to inverse problems, which deals not only with point estimates but with an analogous question for the full posterior distributions. This topic is beyond the scope of this survey: we refer to Kaipio et al. (Reference Kaipio, Kolehmainen, Vauhkonen and Somersalo1999), Neubauer and Pikkarainen (Reference Neubauer and Pikkarainen2008), Stuart (Reference Stuart2010), Kolehmainen et al. (Reference Kolehmainen, Lassas, Niinimäki and Siltanen2012), Castillo et al. (Reference Castillo and Nickl2014), Kekkonen, Lassas and Siltanen (Reference Kekkonen, Lassas and Siltanen2016), Burger, Helin and Kekkonen (Reference Burger, Helin and Kekkonen2016b ) and Nickl et al. (Reference Nickl and Söhl2017) for further details.

5.1 Analysis of variational regularization

In the following we will discuss the basic analysis of variational regularization methods; again we try to give a fairly general perspective that covers most of the results in the literature (but due to its generality does not simply reproduce them). Since we focus on the nonlinear regularization we will make the assumption that ${\mathcal{V}}$ is a separable Hilbert space. A first key issue is the existence of minimizers, which of course depends strongly on the choice of the regularization functional $J$ and possibly also the operator $K$ and the fidelity $F$ . As usual the key issues are lower semicontinuity and compactness in some topology. The latter is always obtained by coercivity in a Banach space norm, which is concluded from the boundedness of the fidelity and in particular the regularization functional. Consequently, the type of compactness is always weak or weak-star, since it is derived from the Banach–Alaoglu theorem (Rudin Reference Rudin2006).

A natural assumption to make for an existence proof is the following.

Note that the above assumptions on $K$ and $F$ are reminiscent of the set-up used by Bredies and Pikkarainen (Reference Bredies and Pikkarainen2013) and later by Brinkmann, Burger, Rasch and Sutour (Reference Brinkmann, Burger, Rasch and Sutour2017). An alternative set-up is to use a compactness assumption on $K$ or some condition on the range of $K$ . Moreover, the assumption on $J$ to be the polar of a proper functional implies convexity, which is predominant in most approaches in regularization theory. With these assumptions we can first verify well-posedness of the selection operator.

Proof. If $f\in {\mathcal{R}}_{F}(K)$ then there exists a minimizer $u^{\ast }$ of $F(Ku,f)$ with $J(u^{\ast },\unicode[STIX]{x1D736})<\infty$ . Since the minimization in the definition of ${\mathcal{S}}$ can be restricted to the set of $u$ such that $F(Ku,f)=F(Ku^{\ast },f)=:a$ , we obtain an upper bound on the fidelity. On this non-empty set we look for $u$ with $J(u)\leq J(u^{\ast })=:b$ . Thus, for the set of such $u$ , the norm in ${\mathcal{U}}$ is bounded due to Assumption 5.4 and for each minimizing sequence there exists a weak-star convergent subsequence $u_{n}$ (we can use the metric version of the Banach–Alaoglu theorem due to the assumption of metrizability on bounded sets). Moreover, from our assumptions above it is straightforward to see that $J(\cdot ,\unicode[STIX]{x1D736})$ is sequentially weak-star lower semicontinuous and $F(\cdot ,f)$ is weakly lower semicontinuous. From our assumption on $K$ being the adjoint of $L$ , we see that it is continuous from the weak-star topology of ${\mathcal{U}}$ to the weak topology of ${\mathcal{V}}$ , since for $g\in {\mathcal{V}}$ we have

$$\begin{eqnarray}\langle Ku_{n},g\rangle =\langle u_{n},Lg\rangle\end{eqnarray}$$

and $Lg\in Z$ . As a consequence, the full functional $F(\cdot ,f)+J(\cdot ,\unicode[STIX]{x1D736})$ is weak-star lower semicontinuous. Hence, the weak-star limit of $u_{n}$ is a minimizer, that is, ${\mathcal{S}}$ is not empty. ☐

The next step is to verify well-definedness of the regularization operator.

Proof. In order to obtain an a priori bound, we use the assumption that there exists $\tilde{u}$ with

$$\begin{eqnarray}a:=F(K\tilde{u} ,f)+J(\tilde{u} ,\unicode[STIX]{x1D736})<\infty .\end{eqnarray}$$

Hence, we can restrict the minimization to those $u$ with functional value less than or equal to $a$ . Setting $b=a$ and using the non-negativity of both terms, we obtain the boundedness of the norm on this subset due to Assumption 5.4. The remaining weak-star compactness and lower semicontinuous arguments to verify the existence of a minimizer are analogous to the proof of Lemma 5.5. The convexity of $R(f,\unicode[STIX]{x1D6FC})$ follows from the convexity of the set of minimizers of a convex functional. ☐

In order to verify the generalized stability as well as the convergence of the variational regularization, a further condition on $F$ with respect to the second variable is needed. There are several options, the easiest one being satisfied by standard examples such as squared norms is continuity.

Proof. By definition of the regularization operator we find for $u_{n}\in R(f_{n},\unicode[STIX]{x1D736})$ that, for any $u\in {\mathcal{U}}$ ,

$$\begin{eqnarray}F(Ku_{n},f_{n})+J(u_{n},\unicode[STIX]{x1D736})\leq F(Ku,f_{n})+J(u,\unicode[STIX]{x1D736}).\end{eqnarray}$$

Due to the convergence of $f_{n}$ and the continuity of $F$ in the second argument, the right-hand side in the last estimate is uniformly bounded by some constant $a$ , which again provides uniform bounds for both terms on the left-hand side. Consequently

$$\begin{eqnarray}\Vert u_{n}\Vert _{{\mathcal{U}}}\leq c(a,a,\Vert f_{n}\Vert _{{\mathcal{V}}}).\end{eqnarray}$$

The boundedness of $\Vert f_{n}\Vert _{{\mathcal{V}}}$ and monotone dependence of $c$ yields a uniform bound on $\Vert u_{n}\Vert _{{\mathcal{U}}}$ , thus a weakly converging subsequence. Using lower semicontinuity arguments as in the results above and the continuity of $F$ with respect to the second variable, we see that for the limit $u^{\ast }$ the inequality

$$\begin{eqnarray}\displaystyle F(Ku^{\ast },f)+J(u,\unicode[STIX]{x1D736}) & \leq & \displaystyle \liminf F(Ku_{n_{k}},f_{n_{k}})+J(u_{n_{k}},\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & \leq & \displaystyle \lim F(Ku,f_{n_{k}})+J(u,\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & = & \displaystyle F(Ku,f)+J(u,\unicode[STIX]{x1D736}).\nonumber\end{eqnarray}$$

Hence $u^{\ast }\in R(f,\unicode[STIX]{x1D736}).$  ☐

As mentioned earlier, the type of convergence in Theorem 5.7 corresponds exactly to the type of stability in the Kuratowski limit superior. We finally provide a comment on the convergence of the regularization method. The proof is very similar to the stability result, and an a priori bound is obtained by the estimate

$$\begin{eqnarray}F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})+J(u^{\unicode[STIX]{x1D736}},\unicode[STIX]{x1D736})\leq F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})+J(u^{\dagger },\unicode[STIX]{x1D736})\leq \unicode[STIX]{x1D6FF}+J(u^{\dagger },\unicode[STIX]{x1D736})\end{eqnarray}$$

for $u^{\unicode[STIX]{x1D736}}\in R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736})$ and any element $u^{\dagger }\in {\mathcal{S}}(f,\boldsymbol{\unicode[STIX]{x1D6FC}})$ . Depending on the specific dependence on $\unicode[STIX]{x1D736}$ , some condition on the interplay of the noise level and the limit of $\unicode[STIX]{x1D736}$ is needed in order to pass to the limit in

$$\begin{eqnarray}J(u^{\unicode[STIX]{x1D736}},\unicode[STIX]{x1D736})\leq \unicode[STIX]{x1D6FF}+J(u^{\dagger },\unicode[STIX]{x1D736}).\end{eqnarray}$$

An abstract condition as $\unicode[STIX]{x1D736}$ converges to $\unicode[STIX]{x1D736}^{\ast }$ outside $A$ is

$$\begin{eqnarray}\lim _{\unicode[STIX]{x1D736}\rightarrow \unicode[STIX]{x1D736}^{\ast }}\frac{\unicode[STIX]{x1D6FF}}{J(u^{\dagger },\unicode[STIX]{x1D736})}=0,\end{eqnarray}$$

in which case

$$\begin{eqnarray}\limsup _{\unicode[STIX]{x1D736}\rightarrow \unicode[STIX]{x1D736}^{\ast }}\displaystyle \frac{J(u^{\unicode[STIX]{x1D736}},\unicode[STIX]{x1D736})}{J(u^{\dagger },\unicode[STIX]{x1D736})}\leq 1.\end{eqnarray}$$

In the standard case $J(u,\unicode[STIX]{x1D736})=\unicode[STIX]{x1D6FC}J(u)$ the condition is simply $\unicode[STIX]{x1D6FF}/\unicode[STIX]{x1D6FC}\rightarrow 0$ . Hence, for such parameter choices, variational regularization methods indeed define convergent regularization operators.

5.2 Error estimates

When it comes to the solution of ill-posed, inverse problems, an important question to address is the question of how errors in the measurement data are being propagated in the regularization process; in particular, convergence with respect to the noise level $\unicode[STIX]{x1D6FF}$ and the rate of convergence are of major interest. Following up on Definition 4.8, we look into $D$ -convergence in the case of $D$ being a Bregman distance.

In order to derive error estimates, we restrict ourselves to the following smoothness class ${\mathcal{M}}_{\unicode[STIX]{x1D708}}$ . Given some unknown ground truth solution $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ , we ensure $u^{\dagger }\in {\mathcal{R}}(R(\cdot ,\unicode[STIX]{x1D736}))$ , that is, we have to ensure that there exists data $f_{\unicode[STIX]{x1D736}}^{\dagger }$ such that $u^{\dagger }\in R(f_{\unicode[STIX]{x1D736}}^{\dagger },\unicode[STIX]{x1D736})$ is a solution of the corresponding variational regularization problem.

Definition 5.8 ((variational) range condition).

We say an element $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ for $f\in {\mathcal{R}}_{F}(K)$ satisfies the range condition if $u^{\dagger }\in {\mathcal{R}}(R(\cdot ,\unicode[STIX]{x1D736}))$ . If $K\in {\mathcal{L}}({\mathcal{U}},{\mathcal{V}})$ , $F$ is convex and Fréchet-differentiable with respect to its first argument, and $J(\cdot ,\unicode[STIX]{x1D736})$ is proper, convex and lower semicontinuous, then this is equivalent to the existence of $p^{\dagger }\in \unicode[STIX]{x2202}J(u^{\dagger },\unicode[STIX]{x1D736})$ and $f_{\unicode[STIX]{x1D736}}^{\dagger }\in {\mathcal{V}}$ such that

(RC) $$\begin{eqnarray}p^{\dagger }=-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger }).\end{eqnarray}$$

From now on we assume $K\in {\mathcal{L}}({\mathcal{U}},{\mathcal{V}})$ , convexity and Fréchet-differentiability of $F$ in its first argument, and properness, convexity and lower semicontinuity of $R(\cdot ,\unicode[STIX]{x1D736})$ for the remainder of this section, which will allow us to use an appropriate optimality condition.

Let us sketch the basic idea in the case of a quadratic fidelity $F(f,g)={\textstyle \frac{1}{2}}\Vert f-g\Vert ^{2}$ with some norm in a Hilbert space and $J(u,\unicode[STIX]{x1D736})=\unicode[STIX]{x1D6FC}J(u)$ . The optimality condition (3.1) is given by

$$\begin{eqnarray}K^{\ast }(Ku^{\unicode[STIX]{x1D6FC}}-f^{\unicode[STIX]{x1D6FF}})+\unicode[STIX]{x1D6FC}p^{\unicode[STIX]{x1D6FC}}=0,\quad p^{\unicode[STIX]{x1D6FC}}\in \unicode[STIX]{x2202}J(u^{\unicode[STIX]{x1D6FC}}).\end{eqnarray}$$

In order to satisfy the range condition for $u^{\dagger }$ we need to assume the existence $f_{\unicode[STIX]{x1D6FC}}^{\dagger }$ such that $p^{\dagger }\in \unicode[STIX]{x2202}J(u^{\dagger })$ and

$$\begin{eqnarray}K^{\ast }(Ku^{\dagger }-f_{\unicode[STIX]{x1D6FC}}^{\dagger })+\unicode[STIX]{x1D6FC}p^{\dagger }=0.\end{eqnarray}$$

We see that this equation implies the condition $p^{\dagger }=K^{\ast }v$ for some $v$ (noting $K^{\ast }(Ku^{\dagger }-f^{\dagger })=0$ ). On the other hand, if this condition is satisfied we can construct $f_{\unicode[STIX]{x1D6FC}}^{\dagger }=f-\unicode[STIX]{x1D6FC}v$ , that is, $p^{\dagger }=K^{\ast }v$ is equivalent to the range condition (RC). An error estimate can then be obtained by subtracting both optimality conditions

$$\begin{eqnarray}K^{\ast }K(u^{\unicode[STIX]{x1D6FC}}-u^{\dagger })+\unicode[STIX]{x1D6FC}(p^{\unicode[STIX]{x1D6FC}}-p^{\dagger })=K^{\ast }(f^{\unicode[STIX]{x1D6FF}}-f_{\unicode[STIX]{x1D6FC}}^{\dagger }).\end{eqnarray}$$

Taking a duality product with $u^{\unicode[STIX]{x1D6FC}}-u^{\dagger }$ yields

$$\begin{eqnarray}\Vert K(u^{\unicode[STIX]{x1D6FC}}-u^{\dagger })\Vert ^{2}+\unicode[STIX]{x1D6FC}D_{J}^{p^{\unicode[STIX]{x1D6FC}}}(u^{\dagger },u^{\unicode[STIX]{x1D6FC}})+\unicode[STIX]{x1D6FC}D_{J}^{p^{\dagger }}(u^{\unicode[STIX]{x1D6FC}},u^{\dagger })=\langle K(u^{\unicode[STIX]{x1D6FC}}-u^{\dagger }),f^{\unicode[STIX]{x1D6FF}}-f_{\unicode[STIX]{x1D6FC}}^{\dagger }\rangle .\end{eqnarray}$$

Applying Young’s inequality on the right-hand side and inserting the special form of $f_{\unicode[STIX]{x1D6FC}}^{\dagger }$ then immediately yields an error estimate (Burger Reference Burger and Hiriart-Urruty2016). Note that we obtain an upper bound on the residual as well as the symmetric Bregman distance

(5.5) $$\begin{eqnarray}D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})=D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{\unicode[STIX]{x1D6FC}}}(u^{\dagger },u^{\unicode[STIX]{x1D6FC}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{\dagger }}(u^{\unicode[STIX]{x1D6FC}},u^{\dagger }).\end{eqnarray}$$

For further interpretations of the error estimates see Burger and Osher (Reference Burger and Osher2004), Burger, Resmerita and He (Reference Burger, Resmerita and He2007b ), Resmerita and Scherzer (Reference Resmerita and Scherzer2006) and Burger (Reference Burger and Hiriart-Urruty2016).

We now want to show that (RC) coincides with the well-known source condition (Chavent and Kunisch Reference Chavent and Kunisch1997, Burger and Osher Reference Burger and Osher2004) for a certain class of fidelity functionals. Before we proceed, we have to define this source condition first.

Definition 5.9 (source condition).

An element $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ for $f\in {\mathcal{R}}_{K}(F)$ satisfies the source condition if

$$\begin{eqnarray}{\mathcal{R}}(K^{\ast })\cap \unicode[STIX]{x2202}J(u^{\dagger },\unicode[STIX]{x1D736})\neq \emptyset .\end{eqnarray}$$

This is equivalent to the existence of $p^{\dagger }\in \unicode[STIX]{x2202}J(u^{\dagger },\unicode[STIX]{x1D736})$ and $v\in {\mathcal{V}}^{\ast }\setminus \{0\}$ such that

(SC) $$\begin{eqnarray}p^{\dagger }=K^{\ast }v.\end{eqnarray}$$

It is obvious that (RC) implies (SC). However, we want to go one step further and show that (RC) and (SC) are actually equivalent conditions for fidelity functionals $F(Ku,f^{\unicode[STIX]{x1D6FF}}):=G(Ku-f^{\unicode[STIX]{x1D6FF}})$ , where $G$ is a Legendre functional. Legendre functionals are defined as follows.

Now we show that (RC) and (SC) are equivalent when $G$ is a Legendre functional.

Proof. ‘ $\boldsymbol{\Rightarrow }$ ’ Condition (RC) trivially implies (SC) if we define

$$\begin{eqnarray}v:=-\unicode[STIX]{x2202}_{x}F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger })=-G^{\prime }(Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }).\end{eqnarray}$$

‘ $\boldsymbol{\Leftarrow }$ ’ The source condition (SC) can be written as

$$\begin{eqnarray}\displaystyle 0 & = & \displaystyle p^{\dagger }-K^{\ast }v,\nonumber\\ \displaystyle \Leftrightarrow 0 & = & \displaystyle p^{\dagger }+K^{\ast }G^{\prime }((G^{\ast })^{\prime }(-v)),\nonumber\end{eqnarray}$$

where $G^{\ast }:{\mathcal{V}}^{\ast }\rightarrow \mathbb{R}\cup \{\infty \}$ denotes the convex conjugate of $G$ . Note that $G^{\ast }$ is also a Legendre functional since ${\mathcal{V}}$ is reflexive (Bauschke, Borwein and Combettes Reference Bauschke, Borwein and Combettes2001, Corollary 5.5), and that the last equality is valid for all $v\in \text{dom}(G)$ due to Bauschke et al. (Reference Bauschke, Borwein and Combettes2001, Theorem 5.9). Hence, if we define

$$\begin{eqnarray}f_{\unicode[STIX]{x1D736}}^{\dagger }:=Ku^{\dagger }-(G^{\ast })^{\prime }(-v),\end{eqnarray}$$

we ensure that the range condition (RC) is satisfied. ☐

The range condition (RC) allows us to derive error estimates in a Bregman distance setting for these very generic variational regularization methods. The following lemma builds the basis by estimating Bregman distances between $u^{\unicode[STIX]{x1D736}}$ and $u^{\dagger }$ in terms of differences in data fidelity.

Lemma 5.13. Let (RC) be satisfied. Then we observe

(5.6) $$\begin{eqnarray}\displaystyle & & \displaystyle D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+D_{F(K\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(u^{\unicode[STIX]{x1D736}},u^{\dagger })+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})\nonumber\\ \displaystyle & & \displaystyle \quad =F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})-F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger })+F(Ku^{\unicode[STIX]{x1D736}},f_{\unicode[STIX]{x1D736}}^{\dagger })-F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})\end{eqnarray}$$

for every $u^{\unicode[STIX]{x1D736}}\in R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736})$ .

Proof. Computing the optimality condition (5.4) of (5.3) and subtracting $p^{\dagger }\in \unicode[STIX]{x2202}J(u^{\dagger },\unicode[STIX]{x1D736})$ from both sides of the equality yields

$$\begin{eqnarray}p_{\unicode[STIX]{x1D736}}-p^{\dagger }=-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})-p^{\dagger },\end{eqnarray}$$

for any $p_{\unicode[STIX]{x1D736}}\in \unicode[STIX]{x2202}J(u^{\unicode[STIX]{x1D736}},\unicode[STIX]{x1D736})$ . Taking a duality product with $u^{\unicode[STIX]{x1D736}}-u^{\dagger }$ then yields

$$\begin{eqnarray}D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\unicode[STIX]{x1D736}},u^{\dagger })=\underbrace{\langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}}),u^{\dagger }-u^{\unicode[STIX]{x1D736}}\rangle }_{=F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})-F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})-D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})}-\langle p^{\dagger },u^{\unicode[STIX]{x1D736}}-u^{\dagger }\rangle .\end{eqnarray}$$

Hence, we conclude

(5.7) $$\begin{eqnarray}\displaystyle & & \displaystyle D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\unicode[STIX]{x1D736}},u^{\dagger })\nonumber\\ \displaystyle & & \displaystyle \quad =F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})-F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})-\langle p^{\dagger },u^{\unicode[STIX]{x1D736}}-u^{\dagger }\rangle .\end{eqnarray}$$

If we now choose $p^{\dagger }=-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger })$ – which is possible since (RC) holds true – we obtain the equality

(5.8) $$\begin{eqnarray}\displaystyle -\langle p^{\dagger },u^{\unicode[STIX]{x1D736}}-u^{\dagger }\rangle & = & \displaystyle \langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger }),u^{\unicode[STIX]{x1D736}}-u^{\dagger }\rangle \nonumber\\ \displaystyle & = & \displaystyle F(Ku^{\unicode[STIX]{x1D736}},f_{\unicode[STIX]{x1D736}}^{\dagger })-F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger })-D_{F(K\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(u^{\unicode[STIX]{x1D736}},u^{\dagger }).\end{eqnarray}$$

Inserting (5.8) into (5.7) then yields (5.6). ☐

Before we proceed, we make the following observation for data fidelities $F$ that are also Bregman distances.

Corollary 5.14. Let $F:{\mathcal{V}}\times {\mathcal{V}}\rightarrow \mathbb{R}$ be a Bregman distance, that is,

$$\begin{eqnarray}F(f,f_{\unicode[STIX]{x1D736}}^{\dagger })=G(f)-G(f_{\unicode[STIX]{x1D736}}^{\dagger })-\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f_{\unicode[STIX]{x1D736}}^{\dagger }\rangle \geq 0,\end{eqnarray}$$

for all $f,f_{\unicode[STIX]{x1D736}}^{\dagger }\in {\mathcal{V}}$ , and some functional $G:{\mathcal{V}}\rightarrow \mathbb{R}$ . Then we already observe

$$\begin{eqnarray}D_{F(\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(f,f^{\unicode[STIX]{x1D6FF}})=F(f,f^{\unicode[STIX]{x1D6FF}}),\end{eqnarray}$$

for all $f,f_{\unicode[STIX]{x1D736}}^{\dagger },f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}}$ .

Proof. We simply compute

$$\begin{eqnarray}\displaystyle D_{F(\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(f,f^{\unicode[STIX]{x1D6FF}}) & = & \displaystyle F(f,f_{\unicode[STIX]{x1D736}}^{\dagger })-F(f^{\unicode[STIX]{x1D6FF}},f_{\unicode[STIX]{x1D736}}^{\dagger })-\langle \unicode[STIX]{x2202}_{x}F(f^{\unicode[STIX]{x1D6FF}},f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f^{\unicode[STIX]{x1D6FF}}\rangle \nonumber\\ \displaystyle & = & \displaystyle G(f)-G(f_{\unicode[STIX]{x1D736}}^{\dagger })-\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f_{\unicode[STIX]{x1D736}}^{\dagger }\rangle \nonumber\\ \displaystyle & & \displaystyle \quad -\,G(f^{\unicode[STIX]{x1D6FF}})+G(f_{\unicode[STIX]{x1D736}}^{\dagger })+\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f^{\unicode[STIX]{x1D6FF}}-f_{\unicode[STIX]{x1D736}}^{\dagger }\rangle \nonumber\\ \displaystyle & & \displaystyle \quad -\,\langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}})-G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f^{\unicode[STIX]{x1D6FF}}\rangle ,\nonumber\\ \displaystyle & = & \displaystyle G(f)-G(f^{\unicode[STIX]{x1D6FF}})-\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f^{\unicode[STIX]{x1D6FF}}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad -\,\langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}})-G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),f-f^{\unicode[STIX]{x1D6FF}}\rangle ,\nonumber\\ \displaystyle & = & \displaystyle D_{G}^{}(f,f^{\unicode[STIX]{x1D6FF}})=F(f,f^{\unicode[STIX]{x1D6FF}}),\nonumber\end{eqnarray}$$

and hence prove the result. ☐

As a consequence, Lemma 5.13 reads as follows for data fidelities that are also Bregman distances.

Lemma 5.15. Let the assumptions of Lemma 5.13 and Corollary 5.14 hold true. Then we have

$$\begin{eqnarray}D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})=\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger })-G^{\prime }(Ku^{\dagger })-(G^{\prime }(f^{\unicode[STIX]{x1D6FF}})-G^{\prime }(Ku^{\unicode[STIX]{x1D736}})),Ku^{\dagger }-Ku^{\unicode[STIX]{x1D736}}\rangle .\end{eqnarray}$$

Proof. From Corollary 5.14 we know that

$$\begin{eqnarray}D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})=F(Ku^{\dagger },Ku^{\unicode[STIX]{x1D736}}),\quad D_{F(K\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(u^{\unicode[STIX]{x1D736}},u^{\dagger })=F(Ku^{\unicode[STIX]{x1D736}},Ku^{\dagger }).\end{eqnarray}$$

Hence, we observe

$$\begin{eqnarray}\displaystyle & & \displaystyle D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+D_{F(K\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(u^{\unicode[STIX]{x1D736}},u^{\dagger })\nonumber\\ \displaystyle & & \displaystyle \quad =F(Ku^{\dagger },Ku^{\unicode[STIX]{x1D736}})+F(Ku^{\unicode[STIX]{x1D736}},Ku^{\dagger })\nonumber\\ \displaystyle & & \displaystyle \quad =D_{G}^{\text{symm}}(Ku^{\dagger },Ku^{\unicode[STIX]{x1D736}})\nonumber\\ \displaystyle & & \displaystyle \quad =\langle G^{\prime }(Ku^{\dagger })-G^{\prime }(Ku^{\unicode[STIX]{x1D736}}),Ku^{\dagger }-Ku^{\unicode[STIX]{x1D736}}\rangle .\nonumber\end{eqnarray}$$

We also discover

$$\begin{eqnarray}\displaystyle & & \displaystyle F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})-F(Ku^{\dagger },f_{\unicode[STIX]{x1D736}}^{\dagger })+F(Ku^{\unicode[STIX]{x1D736}},f_{\unicode[STIX]{x1D736}}^{\dagger })-F(Ku^{\unicode[STIX]{x1D736}},f^{\unicode[STIX]{x1D6FF}})\nonumber\\ \displaystyle & & \displaystyle \quad =G(Ku^{\dagger })-G(f^{\unicode[STIX]{x1D6FF}})-\langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}}),Ku^{\dagger }-f^{\unicode[STIX]{x1D6FF}}\rangle \nonumber\\ \displaystyle & & \displaystyle \qquad -\,(G(Ku^{\dagger })-G(f_{\unicode[STIX]{x1D736}}^{\dagger })-\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\rangle )\nonumber\\ \displaystyle & & \displaystyle \qquad +\,G(Ku^{\unicode[STIX]{x1D736}})-G(f_{\unicode[STIX]{x1D736}}^{\dagger })-\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),Ku^{\unicode[STIX]{x1D736}}-f_{\unicode[STIX]{x1D736}}^{\dagger }\rangle \nonumber\\ \displaystyle & & \displaystyle \qquad -\,(G(Ku^{\unicode[STIX]{x1D736}})-G(f^{\unicode[STIX]{x1D6FF}})-\langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}}),Ku^{\unicode[STIX]{x1D736}}-f^{\unicode[STIX]{x1D6FF}}\rangle )\nonumber\\ \displaystyle & = & \displaystyle \langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}}),Ku^{\unicode[STIX]{x1D736}}-Ku^{\dagger }\rangle +\langle G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),Ku^{\dagger }-Ku^{\unicode[STIX]{x1D736}}\rangle \nonumber\\ \displaystyle & = & \displaystyle \langle G^{\prime }(f^{\unicode[STIX]{x1D6FF}})-G^{\prime }(f_{\unicode[STIX]{x1D736}}^{\dagger }),Ku^{\unicode[STIX]{x1D736}}-Ku^{\dagger }\rangle .\nonumber\end{eqnarray}$$

Combining these two equalities with (5.6) yields the desired result. ☐

Example 5.16. We can use (5.6) to derive the same error estimates presented in Burger et al. (Reference Burger, Resmerita and He2007b ) for the choice $F(Ku,f^{\unicode[STIX]{x1D6FF}})={\textstyle \frac{1}{2}}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2},$ where ${\mathcal{H}}$ is a Hilbert space. In this case we observe

$$\begin{eqnarray}D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{\dagger },u^{\unicode[STIX]{x1D736}})=D_{F(K\cdot ,f_{\unicode[STIX]{x1D736}}^{\dagger })}^{}(u^{\unicode[STIX]{x1D736}},u^{\dagger })=\displaystyle \frac{1}{2}\Vert K(u^{\dagger }-u^{\unicode[STIX]{x1D736}})\Vert _{{\mathcal{H}}}^{2}.\end{eqnarray}$$

Hence, equation (5.6) reads as

(5.9) $$\begin{eqnarray}\displaystyle & & \displaystyle \Vert K(u^{\dagger }-u^{\unicode[STIX]{x1D736}})\Vert _{{\mathcal{H}}}^{2}+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+\displaystyle \frac{1}{2}\Vert Ku^{\unicode[STIX]{x1D736}}-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad =\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}+\displaystyle \frac{1}{2}\Vert Ku^{\unicode[STIX]{x1D736}}-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}-\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}.\end{eqnarray}$$

If we make use of the estimate

$$\begin{eqnarray}\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\leq \displaystyle \frac{1}{2}\Vert f-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\leq \unicode[STIX]{x1D6FF},\end{eqnarray}$$

then equality (5.9) becomes the inequality

$$\begin{eqnarray}\displaystyle & & \displaystyle \Vert K(u^{\dagger }-u^{\unicode[STIX]{x1D736}})\Vert _{{\mathcal{H}}}^{2}+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+\displaystyle \frac{1}{2}\Vert Ku^{\unicode[STIX]{x1D736}}-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leq \,\unicode[STIX]{x1D6FF}+\displaystyle \frac{1}{2}\Vert K(u^{\unicode[STIX]{x1D736}}-u^{\dagger })+(Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger })\Vert _{{\mathcal{H}}}^{2}-\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leq \,\unicode[STIX]{x1D6FF}+\Vert K(u^{\dagger }-u^{\unicode[STIX]{x1D736}})\Vert _{{\mathcal{H}}}^{2}+\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}-\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}.\nonumber\end{eqnarray}$$

Subtracting $\Vert K(u^{\dagger }-u^{\unicode[STIX]{x1D736}})\Vert _{{\mathcal{H}}}^{2}$ on both sides of the inequality then yields the error estimate

(5.10) $$\begin{eqnarray}\displaystyle D_{J(\cdot ,\unicode[STIX]{x1D736})}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+\displaystyle \frac{1}{2}\Vert Ku^{\unicode[STIX]{x1D736}}-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\leq \unicode[STIX]{x1D6FF}+\displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}. & & \displaystyle\end{eqnarray}$$

We want to emphasize that the constant ${\textstyle \frac{1}{2}}\Vert Ku^{\dagger }-f_{\unicode[STIX]{x1D736}}^{\dagger }\Vert _{{\mathcal{H}}}^{2}$ on the right-hand-side of the inequality depends on the choice of $\unicode[STIX]{x1D736}$ . From Remark 5.10 and the proof of Theorem 5.12 it follows that if we consider regularizations of the form $J(u,\unicode[STIX]{x1D736})=\unicode[STIX]{x1D6FC}J(u)$ , the source condition (SC) and the range condition (RC) are linked via the relation $f_{\unicode[STIX]{x1D736}}^{\dagger }=Ku^{\dagger }+\unicode[STIX]{x1D6FC}v$ , where $v$ is the source condition element for $\unicode[STIX]{x1D6FC}=1$ , i.e.  $K^{\ast }v\in \unicode[STIX]{x2202}J(u^{\dagger },1)=\unicode[STIX]{x2202}J(u^{\dagger })$ . In this setting, the error estimate (5.10) then reads as

$$\begin{eqnarray}D_{J}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})+\displaystyle \frac{1}{2\unicode[STIX]{x1D6FC}}\Vert Ku^{\unicode[STIX]{x1D736}}-f^{\unicode[STIX]{x1D6FF}}\Vert _{{\mathcal{H}}}^{2}\leq \displaystyle \frac{\unicode[STIX]{x1D6FF}}{\unicode[STIX]{x1D6FC}}+\displaystyle \frac{\unicode[STIX]{x1D6FC}}{2}\Vert v\Vert _{{\mathcal{H}}}^{2}.\end{eqnarray}$$

Hence, choosing $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D6FF})=\sqrt{2\unicode[STIX]{x1D6FF}}/\Vert v\Vert _{{\mathcal{H}}}$ then yields $D_{J}^{\text{symm}}(u^{\dagger },u^{\unicode[STIX]{x1D736}})=O(\sqrt{\unicode[STIX]{x1D6FF}})$ .

There are various routes and generalizations that can be taken from these types of estimates, for example to weaker source conditions with the concepts of approximate or variational source conditions (Schuster, Kaltenbacher, Hofmann and Kazimierski Reference Schuster, Kaltenbacher, Hofmann and Kazimierski2012, Flemming and Hofmann Reference Flemming and Hofmann2010, Flemming Reference Flemming2013, Flemming Reference Flemming2017b , Hohage and Weidling Reference Hohage and Weidling2017), improved estimates for stronger conditions (Resmerita Reference Resmerita2005, Grasmair Reference Grasmair2013), or large noise that is not necessarily in ${\mathcal{V}}$ (Burger et al. Reference Burger, Helin and Kekkonen2016b ). In special cases such as $\ell ^{1}$ -regularization, improved results can be obtained, because the effective finite-dimensionality implies that this case is almost well-posed (Grasmair, Scherzer and Haltmeier Reference Grasmair, Scherzer and Haltmeier2011, Grasmair Reference Grasmair2011, Burger, Flemming and Hofmann Reference Burger, Flemming and Hofmann2013a , Flemming, Hofmann and Veselić Reference Flemming, Hofmann and Veselić2015, Flemming, Hofmann and Veselić Reference Flemming, Hofmann and Veselić2016, Flemming and Gerth Reference Flemming and Gerth2017). Converse results have also been recently obtained (Flemming Reference Flemming2017a , Hohage and Weidling Reference Hohage and Weidling2017).

5.3 Variational eigenvalue problems

The standard tool for the analysis of linear regularization methods is singular value decomposition. In the case of nonlinear regularization no analogue of singular values and singular vectors was known for a long time. A generalization for nonlinear variational methods was made in Benning and Burger (Reference Benning and Burger2013), which we discuss in the following. We generalize singular vectors as eigenvectors of the variational regularization operator $R$ as defined in Definition 5.2, that is, we look for functions $u_{\unicode[STIX]{x1D706}}$ that satisfy

(5.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}u_{\unicode[STIX]{x1D706}}\in R(\unicode[STIX]{x1D70E}Ku_{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D736}), & & \displaystyle\end{eqnarray}$$

for constants $\unicode[STIX]{x1D706},\unicode[STIX]{x1D70E}\in [0,\infty )$ , typically $\unicode[STIX]{x1D70E}=1$ . For simplicity we focus on the case where $\unicode[STIX]{x1D736}=\unicode[STIX]{x1D6FC}$ is a scalar, and $F(Ku,f^{\unicode[STIX]{x1D6FF}})=G(Ku-f^{\unicode[STIX]{x1D6FF}})$ , where $G$ is a Legendre functional for the remainder of this section. If we consider the optimality condition (5.4) of (5.3), we immediately observe that any $u_{\unicode[STIX]{x1D706}}$ satisfying (5.11) also has to satisfy

(5.12) $$\begin{eqnarray}\displaystyle -K^{\ast }G^{\prime }((\unicode[STIX]{x1D706}-\unicode[STIX]{x1D70E})Ku_{\unicode[STIX]{x1D706}})\in \unicode[STIX]{x2202}J(\unicode[STIX]{x1D706}u_{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D6FC}). & & \displaystyle\end{eqnarray}$$

We now assume that both $G^{\prime }$ and $\unicode[STIX]{x2202}J$ are homogeneous in the sense that they satisfy $G^{\prime }(cu)=s_{1}(c)G^{\prime }(u)$ and $\unicode[STIX]{x2202}J(cu,\unicode[STIX]{x1D736})=s_{2}(c,\unicode[STIX]{x1D6FC})\unicode[STIX]{x2202}J(u)$ for constants $c\in \mathbb{R}$ and functions $s_{1},s_{2}:\mathbb{R}\rightarrow \mathbb{R}$ . Then (5.12) simplifies to

(5.13) $$\begin{eqnarray}\displaystyle -\displaystyle \frac{s_{1}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D70E})}{s_{2}(\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6FC})}K^{\ast }G^{\prime }(Ku_{\unicode[STIX]{x1D706}})\in \unicode[STIX]{x2202}J(\unicode[STIX]{x1D70E}u_{\unicode[STIX]{x1D706}}). & & \displaystyle\end{eqnarray}$$

Equation (5.13) paves the way for the following definition of generalized singular vectors.

Definition 5.17 (generalized singular system).

Let $\{u_{\unicode[STIX]{x1D70E}},v_{\unicode[STIX]{x1D70E}},\unicode[STIX]{x1D70E}\}$ satisfy

(5.14) $$\begin{eqnarray}\displaystyle Ku_{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}v_{\unicode[STIX]{x1D70E}}\quad \text{and}\quad K^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}})\in \unicode[STIX]{x2202}J(\unicode[STIX]{x1D70E}u_{\unicode[STIX]{x1D70E}}) & & \displaystyle\end{eqnarray}$$

for $\unicode[STIX]{x1D70E}>0$ . Then $\{u_{\unicode[STIX]{x1D70E}},v_{\unicode[STIX]{x1D70E}},\unicode[STIX]{x1D70E}\}$ is called a generalized singular system.

Remark 5.18. For

$$\begin{eqnarray}G(v)=\displaystyle \frac{1}{2}\Vert v\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2}\quad \text{and}\quad J(u,\unicode[STIX]{x1D6FC})=\displaystyle \frac{\unicode[STIX]{x1D6FC}}{2}\Vert u\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2},\end{eqnarray}$$

this definition is consistent with the classical singular vector theory for compact operators.

Example 5.19. Suppose

$$\begin{eqnarray}G(v)=\displaystyle \frac{1}{2}\Vert v\Vert _{L^{2}([0,1])^{2}},\quad K:\text{BV}([0,1])\rightarrow L^{2}([0,1])\end{eqnarray}$$

is the embedding operator and $J(u,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}\text{TV}_{\ast }(u)$ , where $\text{TV}_{\ast }$ denotes the (one-dimensional) total variation with Dirichlet-zero boundary conditions. Benning and Burger (Reference Benning and Burger2013) have shown that Haar wavelets are generalized singular vectors of $\text{TV}_{\ast }$ . Precisely, we have $v_{\unicode[STIX]{x1D70E}_{n},k}=\unicode[STIX]{x1D70E}_{n}u_{\unicode[STIX]{x1D70E}_{n},k}\in \unicode[STIX]{x2202}\text{TV}_{\ast }(\unicode[STIX]{x1D70E}_{n}u_{\unicode[STIX]{x1D70E}_{n},k})=\unicode[STIX]{x2202}\text{TV}_{\ast }(u_{\unicode[STIX]{x1D70E}_{n},k})$ for $\unicode[STIX]{x1D70E}_{n}:=2^{-(n+4)/2}$ and $u_{\unicode[STIX]{x1D70E}_{n},k}$ defined by

$$\begin{eqnarray}u_{\unicode[STIX]{x1D70E}_{n},k}(x):=2^{n/2}\unicode[STIX]{x1D6F9}(2^{n}x-j)\quad \text{with }\unicode[STIX]{x1D6F9}(x):=\left\{\begin{array}{@{}ll@{}}1\quad & x\in [0,1/2),\\ -1\quad & x\in [1/2,1),\\ 0\quad & \text{else}.\end{array}\right.\end{eqnarray}$$

The singular value $\unicode[STIX]{x1D70E}_{n}$ is determined via (5.14). The dual singular vector $v_{\unicode[STIX]{x1D70E}_{n},k}$ has to satisfy $v_{\unicode[STIX]{x1D70E}_{n},k}\in \unicode[STIX]{x2202}\text{TV}_{\ast }(\unicode[STIX]{x1D70E}_{n}u_{\unicode[STIX]{x1D70E}_{n},k})=\text{TV}_{\ast }(u_{\unicode[STIX]{x1D70E}_{n},k})$ . If we make use of $v_{\unicode[STIX]{x1D70E}_{n},k}=u_{\unicode[STIX]{x1D70E}_{n},k}/\unicode[STIX]{x1D70E}_{n}$ and take a dual product with $u_{\unicode[STIX]{x1D70E}_{n},k}$ , we immediately observe $\unicode[STIX]{x1D70E}_{n}=\Vert u_{\unicode[STIX]{x1D70E}_{n},k}\Vert _{L^{2}([0,1])}^{2}/\text{TV}_{\ast }(u_{\unicode[STIX]{x1D70E}_{n},k})$ . In Figure 5.1 we see the Haar wavelet $u_{\unicode[STIX]{x1D70E}_{1},1/2}$ and its scaled version $v_{\unicode[STIX]{x1D70E}_{1},1/2}=\unicode[STIX]{x1D70E}_{1}u_{\unicode[STIX]{x1D70E}_{1},1/2}=2^{-5/2}u_{\unicode[STIX]{x1D70E}_{1},1/2}$ .

Figure 5.1.

The Haar wavelet $u_{\unicode[STIX]{x1D70E}_{1},1/2}$ and its scaled version $v_{\unicode[STIX]{x1D70E}_{1},1/2}$ . Benning and Burger (Reference Benning and Burger2013) have shown that together with $\unicode[STIX]{x1D70E}_{1}=2^{-5/2}$ they form a generalized singular system in the sense of (5.14) with $K$ being the identity in $L^{2}$ .

The generalized singular system is defined so that (5.13) and (5.14) coincide for

(5.15) $$\begin{eqnarray}\displaystyle s_{1}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D70E})=-s_{2}(\unicode[STIX]{x1D706}/\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6FC}). & & \displaystyle\end{eqnarray}$$

Hence, if we choose $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D706}$ such that (5.15) holds true, we already know that (5.11) is satisfied for these particular choices of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D6FC}$ .

Example 5.20. For

$$\begin{eqnarray}G(v)=\displaystyle \frac{1}{2}\Vert v\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2}\quad \text{and}\quad J(u,\unicode[STIX]{x1D6FC})=\displaystyle \frac{\unicode[STIX]{x1D6FC}}{2}\Vert u\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2},\end{eqnarray}$$

with $\unicode[STIX]{x1D6F4}$ and $\unicode[STIX]{x1D6FA}$ being domains in $\mathbb{R}^{d_{1}}$ and $\mathbb{R}^{d_{2}}$ , respectively, we observe $s_{1}(x)=x$ and $s_{2}(x,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}x$ . Hence, (5.15) simplifies to $\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D706}=(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D706})/\unicode[STIX]{x1D70E}$ . Solving for $\unicode[STIX]{x1D706}$ then yields

$$\begin{eqnarray}\unicode[STIX]{x1D706}=\displaystyle \frac{\unicode[STIX]{x1D70E}}{\unicode[STIX]{x1D70E}^{2}+\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

which perfectly coincides with the singular value decomposition representation of Tikhonov regularization.

Example 5.21. For $G(v)={\textstyle \frac{1}{2}}\Vert v\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2}$ and $J(u,\unicode[STIX]{x1D736})=\unicode[STIX]{x1D6FC}\text{TV}(u)$ we have $s_{1}(x)=x$ and $s_{2}(x,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}$ . Consequently, (5.15) solved for $\unicode[STIX]{x1D706}$ reads as

$$\begin{eqnarray}\unicode[STIX]{x1D706}=\displaystyle \frac{1-\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D70E}}.\end{eqnarray}$$

This eigenvalue of this particular regularization operator is consistent with classical singular value theory in the sense that it satisfies $\lim _{\unicode[STIX]{x1D6FC}\downarrow 0}\unicode[STIX]{x1D706}=1/\unicode[STIX]{x1D70E}$ .

An interesting observation from Examples 5.20 and 5.21 is that $\unicode[STIX]{x1D6FC}>0$ automatically implies $\unicode[STIX]{x1D706}<1/\unicode[STIX]{x1D70E}$ (unless $u_{\unicode[STIX]{x1D70E}}\in \ker (J)$ ). This implies that there is always a systematic error when it comes to recovering singular vectors with variational regularization methods that have quadratic fidelity. This is also true for input data that are not given in terms of a singular vector: see Benning and Burger (Reference Benning and Burger2013, Theorem 7). In the next section we see that iterative regularization methods can overcome this systematic reconstruction bias.

6.1 Bregman iteration

The concept of Bregman iteration – also known as proximal minimization algorithm – introduces an iteration into the variational regularization framework by replacing the regularization functional $J(u,\unicode[STIX]{x1D736})$ with the corresponding generalized Bregman distance $D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p}(u,v)$ , for $v\in {\mathcal{U}}$ and $p\in \unicode[STIX]{x2202}J(v)$ . For the choice

$$\begin{eqnarray}J(u,\unicode[STIX]{x1D736})=\displaystyle \frac{\unicode[STIX]{x1D736}}{2}\Vert u\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2},\end{eqnarray}$$

it is also known as iterated Tikhonov regularization, which dates back to the works of Kryanev (Reference Kryanev1974), further analysed in Groetsch (Reference Groetsch1977) and Thomas King and Chillingworth (Reference Thomas King and Chillingworth1979), for example. The extension to more general choices of Bregman distances was first proposed by Censor and Zenios (Reference Censor and Zenios1992), shortly followed by Teboulle (Reference Teboulle1992), and has since been the subject of extensive research (Eckstein Reference Eckstein1993, Kiwiel Reference Kiwiel1997). Notably, Osher et al. (Reference Osher, Burger, Goldfarb, Xu and Yin2005) have extended it to generalized Bregman distances that allow for subdifferentiable rather than differentiable functionals. Note that in such cases there is no one-to-one relation between $u^{k-1}$ and its subgradient $p^{k-1}$ , so we set $v^{k-1}=(u^{k-1},p^{k-1})$ . With a set-valued iteration operator, the Bregman iteration can be written as

$$\begin{eqnarray}\displaystyle u^{k}\in R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736}) & = & \displaystyle \underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\{F(Ku,f^{\unicode[STIX]{x1D6FF}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u,u^{k-1})\},\nonumber\\ \displaystyle p^{k} & = & \displaystyle p^{k-1}-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}),\vphantom{\operatorname{argmin}_{u\in {\mathcal{U}}}\{F(Ku,f^{\unicode[STIX]{x1D6FF}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u,u^{k-1})\}}\nonumber\end{eqnarray}$$

for $p^{0}\in \unicode[STIX]{x2202}J(u^{0},\unicode[STIX]{x1D736})$ . The entire method is summarized in Algorithm 1.

Remark 6.2. The update for the subgradient can also be written as

(6.3) $$\begin{eqnarray}\displaystyle p^{k}=p^{0}-\mathop{\sum }_{n=1}^{k}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f^{\unicode[STIX]{x1D6FF}}). & & \displaystyle\end{eqnarray}$$

Hence, we can rewrite the primal update to

(6.4) $$\begin{eqnarray}\displaystyle & & \displaystyle R_{I}(f^{\unicode[STIX]{x1D6FF}},\{u^{n}\}_{n=1}^{k-1},p^{0},\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\biggl\{F(Ku,f^{\unicode[STIX]{x1D6FF}})+J(u,\unicode[STIX]{x1D736})-\langle p^{0}-\mathop{\sum }_{n=1}^{k-1}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f^{\unicode[STIX]{x1D6FF}}),u\rangle \biggr\}.\end{eqnarray}$$

In the following we want to recall (or derive) a few important properties of Algorithm 1. We start with a trivial result establishing monotonic decrease of the data fidelity.

Corollary 6.3 (monotonic decrease of data fidelity).

Suppose that $u^{0}$ satisfies $F(Ku^{0},f^{\unicode[STIX]{x1D6FF}})<\infty$ . Then the iterates of Algorithm 1 satisfy

$$\begin{eqnarray}F(Ku^{k+1},f^{\unicode[STIX]{x1D6FF}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{k+1},u^{k})\leq F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}),\end{eqnarray}$$

and

$$\begin{eqnarray}\lim _{k\rightarrow \infty }D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{k+1},u^{k})=0,\end{eqnarray}$$

for $u^{k}\in R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})$ and all $k\in \mathbb{N}$ .

Proof. The first statement follows trivially from the convexity of $F$ (in its first argument) and $J$ , and the fact that $u^{k+1}$ is a minimizer of $E(u):=F(Ku,f^{\unicode[STIX]{x1D6FF}})+D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u,u^{k})$ . The first statement then implies

$$\begin{eqnarray}\displaystyle \mathop{\sum }_{k=0}^{N-1}D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{k+1},u^{k}) & \leq & \displaystyle F(Ku^{0},f^{\unicode[STIX]{x1D6FF}})-F(Ku^{N},f^{\unicode[STIX]{x1D6FF}})\nonumber\\ \displaystyle & \leq & \displaystyle F(Ku^{0},f^{\unicode[STIX]{x1D6FF}})<\infty .\nonumber\end{eqnarray}$$

Taking the limit $N\rightarrow \infty$ then yields the second statement. ☐

If we want to show that the Bregman iteration is a convergent regularization method in the sense of Definition 4.8, a first step towards this result would be the following monotonicity lemma.

Lemma 6.4 (Fejér monotonicity of Algorithm 1).

Let $f\in {\mathcal{R}}_{F}(K)$ , $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ and let $f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}}$ with $F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}$ . We further assume that the iterates of Algorithm 1 satisfy Definition 6.1 for $\unicode[STIX]{x1D702}=1$ . Then the iterates also satisfy the strict Fejér monotonicity

$$\begin{eqnarray}D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{\dagger },u^{k})<D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u^{\dagger },u^{k-1}),\end{eqnarray}$$

for $u^{k}\in R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})$ and all $k<k^{\ast }$ .

Proof. Through straightforward computations we obtain

$$\begin{eqnarray}\displaystyle D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{\dagger },u^{k})-D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u^{\dagger },u^{k-1}) & = & \displaystyle \underbrace{-D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u^{k},u^{k-1})}_{{<}0}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\langle p^{k}-p^{k-1},u^{\dagger }-u^{k}\rangle \nonumber\\ \displaystyle & \leq & \displaystyle \langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}),u^{\dagger }-u^{k}\rangle \nonumber\\ \displaystyle & \leq & \displaystyle (\unicode[STIX]{x1D6FF}-F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}))\nonumber\\ \displaystyle & {<} & \displaystyle 0\nonumber\end{eqnarray}$$

for $k<k^{\ast }$ , where we have made use of the convexity of $F$ in its first argument, and $F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})\leq F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}$ . ☐

Corollary 6.5. Let $f\in {\mathcal{R}}_{F}(K)$ and $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ . Then the iterates of Algorithm 1 satisfy

(6.5) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{k=0}^{\infty }F(Ku^{k},f)<\infty & & \displaystyle\end{eqnarray}$$

for $\unicode[STIX]{x1D6FF}=0$ (and, thus, $f^{\unicode[STIX]{x1D6FF}}=f$ ) and $u^{0}$ (with $p^{0}\in \unicode[STIX]{x2202}J(u^{0},\unicode[STIX]{x1D736})$ ) chosen such that $D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{0}}(u^{\dagger },u^{0})<\infty$ .

Proof. For $\unicode[STIX]{x1D6FF}=0$ we conclude

$$\begin{eqnarray}F(Ku^{k},f)\leq D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k-1}}(u^{\dagger },u^{k-1})-D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k}}(u^{\dagger },u^{k})\end{eqnarray}$$

from Lemma 6.4. Summing up from $k=0$ to some $k=k^{\ast }$ therefore yields

$$\begin{eqnarray}\mathop{\sum }_{k=0}^{k^{\ast }}F(Ku^{k},f)\leq D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{0}}(u^{\dagger },u^{0})-D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k^{\ast }}}(u^{\dagger },u^{k^{\ast }})\leq D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{0}}(u^{\dagger },u^{0})<\infty .\end{eqnarray}$$

Taking the limit $k^{\ast }\rightarrow \infty$ yields the assertion. ☐

Remark 6.6. Given the continuity of $F$ and $K\in {\mathcal{L}}({\mathcal{U}},{\mathcal{V}})$ , equation (6.5) already implies

(6.6) $$\begin{eqnarray}\displaystyle Ku^{k}\underset{\unicode[STIX]{x1D70F}_{{\mathcal{V}}}}{\longrightarrow }f, & & \displaystyle\end{eqnarray}$$

if $\unicode[STIX]{x1D70F}_{{\mathcal{V}}}$ is an appropriate topology on ${\mathcal{V}}$ related to $F$ .

Proof. We know $D_{J(\cdot ,\unicode[STIX]{x1D736})}^{p^{k^{\ast }}}(u,u^{k^{\ast }})\geq 0$ for all $u\in {\mathcal{U}}$ and $u^{k^{\ast }}\in R_{I}(f,v^{k^{\ast }-1},\unicode[STIX]{x1D736})$ , since $J$ is convex; this in particular holds true for any $\hat{u} \in \{u\mid Ku=f\}$ . Hence, we observe

$$\begin{eqnarray}\displaystyle J(u^{k^{\ast }}) & \leq & \displaystyle J(\hat{u} )-\langle p^{k^{\ast }},\hat{u} -u^{k^{\ast }}\rangle \nonumber\\ \displaystyle & = & \displaystyle J(\hat{u} )-\langle p^{0},\hat{u} -u^{k^{\ast }}\rangle +\mathop{\sum }_{n=1}^{k^{\ast }}\langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f),\hat{u} -u^{k^{\ast }}\rangle ,\nonumber\\ \displaystyle & \leq & \displaystyle J(\hat{u} )-\langle q^{0},\underbrace{K\hat{u} -Ku^{k^{\ast }}}_{=0}\rangle +\mathop{\sum }_{n=1}^{k^{\ast }}\langle \unicode[STIX]{x2202}_{x}F(Ku^{n},f),\vphantom{K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f),\hat{u} -u^{k^{\ast }}}\underbrace{K\hat{u} -Ku^{k^{\ast }}}_{=0}\vphantom{K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f),\hat{u} -u^{k^{\ast }}}\rangle ,\nonumber\\ \displaystyle & = & \displaystyle J(\hat{u} ),\nonumber\end{eqnarray}$$

for the substitution $p^{0}:=K^{\ast }q^{0}$ , which is possible since $p^{0}\in {\mathcal{R}}(K^{\ast })$ . Here we have made use of equation (6.3). Consequently, we conclude $u^{k^{\ast }}\in {\mathcal{S}}(f,\unicode[STIX]{x1D736})$ . ☐

In the limiting case $k^{\ast }\rightarrow \infty$ the selection is less clear: one cannot prove in general that the limit minimizes $J$ . To make this more apparent consider the case of a least-squares fidelity functional $F(Ku,f)={\textstyle \frac{1}{2}}\Vert Ku-f\Vert ^{2}$ , with the initial value being a minimizer of the regularization, i.e.  $p^{0}=0$ . Then the estimate, as in the last proof (at arbitrary index $m$ ) with $\hat{u} =u^{\dagger }$ , becomes

$$\begin{eqnarray}\displaystyle J(u^{m}) & \leq & \displaystyle J(u^{\dagger })-\mathop{\sum }_{n=1}^{m}\langle K^{\ast }(f^{\unicode[STIX]{x1D6FF}}-Ku^{n}),\hat{u} -u^{m}\rangle \nonumber\\ \displaystyle & = & \displaystyle J(\hat{u} )-\mathop{\sum }_{n=1}^{m}\langle f^{\unicode[STIX]{x1D6FF}}-Ku^{n},f-Ku^{m}\rangle .\nonumber\end{eqnarray}$$

Using Young’s inequality and monotonicity of the residual ( $\Vert Ku^{m}-f^{\unicode[STIX]{x1D6FF}}\Vert \leq \Vert Ku^{n}-f^{\unicode[STIX]{x1D6FF}}\Vert$ ), we conclude that

$$\begin{eqnarray}J(u^{m})\leq J(u^{\dagger })+\displaystyle \frac{3}{2}\mathop{\sum }_{n=1}^{m}\Vert Ku^{n}-f^{\unicode[STIX]{x1D6FF}}\Vert ^{2}+m\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

Summing the estimate in the proof of the Fejér monotonicity, we further find

$$\begin{eqnarray}\mathop{\sum }_{n=1}^{m}\Vert Ku^{n}-f^{\unicode[STIX]{x1D6FF}}\Vert ^{2}\leq D_{J}^{p^{0}}(u^{\dagger },u^{0})=J(u^{\dagger }).\end{eqnarray}$$

Thus, we find

$$\begin{eqnarray}J(u^{k^{\ast }})\leq \displaystyle \frac{5}{2}J(u^{\dagger })+k^{\ast }\unicode[STIX]{x1D6FF}.\end{eqnarray}$$

Since for the discrepancy principle one can show that $k^{\ast }\unicode[STIX]{x1D6FF}^{2}\rightarrow 0$ in the limit $\unicode[STIX]{x1D6FF}\rightarrow 0$ (Osher et al. Reference Osher, Burger, Goldfarb, Xu and Yin2005), the limit of the regularization has a functional value $J$ bounded by ${\textstyle \frac{5}{2}}J(u^{\dagger })$ . Using a more careful argument based on Young’s inequality, this upper bound can be decreased to $2J(u^{\dagger })$ but not to $J(u^{\dagger })$ . On the other hand this might be advantageous, since an estimate of $J(u^{k^{\ast }})$ smaller than $J(u^{\dagger })$ might mean a bias depending on $J$ , since the value of the regularization functional is actually underestimated. In the case of total variation, for example, this means that the contrast is underestimated by variational methods, which is improved by iterative regularization (Osher et al. Reference Osher, Burger, Goldfarb, Xu and Yin2005). To conclude this section, we show numerical results of Bregman-iterative regularization in the context of deconvolution, which demonstrates the effect of total variation regularization.

Figure 6.1.

(a) Image $u^{\dagger }\in \mathbb{R}^{400\times 300}$ of Pixel, a Gambian pouched rat. (b) Degraded and noisy version $f^{\unicode[STIX]{x1D6FF}}\in \mathbb{R}^{400\times 300}$ of the original image $u^{\dagger }$ . The degradation stems from a discretized version of the convolution (see (6.7)) with periodic boundary conditions and the convolution kernel depicted in the bottom left corner of (b). (c–h) Different iterates of Algorithm 1 for $F(Ku,f^{\unicode[STIX]{x1D6FF}})={\textstyle \frac{1}{2}}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\mathbb{R}^{2})}^{2}$ , $J(u,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}\text{TV}(u)$ and $\unicode[STIX]{x1D6FC}=1/4$ . The 96th iterate visualized in (h) is the first that violates Definition 6.1, for $\unicode[STIX]{x1D702}=1$ and $\unicode[STIX]{x1D6FF}=5.95$ .

Example 6.8. We consider the inverse problem of the convolution operation, that is, $Ku=f$ with

(6.7) $$\begin{eqnarray}\displaystyle (Ku)(y):=\int _{\mathbb{R}^{2}}u(x)h(x-y)\,\text{d}x, & & \displaystyle\end{eqnarray}$$

which is therefore also known as deconvolution. Here, $h$ denotes the convolution kernel, that we assume to be known a priori. Since we cannot expect to know $f$ but just $f^{\unicode[STIX]{x1D6FF}}$ with $F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}$ , we need to approximate the inverse problem solution through regularization. In Figure 6.1 we can see selected iterates of Algorithm 1 for a single parameter $\unicode[STIX]{x1D6FC}=1/4$ , the data fidelity term

$$\begin{eqnarray}F(Ku,f^{\unicode[STIX]{x1D6FF}})=\displaystyle \frac{1}{2}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\mathbb{R}^{2})}^{2},\end{eqnarray}$$

and the regularization functional $J(u,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}\text{TV}(u)$ . The data $f^{\unicode[STIX]{x1D6FF}}=f+n$ are the sum of $f$ , created via a discretized version of the exact forward model (6.7), and noise $n\in {\mathcal{N}}(0,0.05)$ . For the particular example used here, the fidelity-noise-bound is $F(f,f^{\unicode[STIX]{x1D6FF}})=5.95$ . The inner variational regularization method is solved via the primal–dual hybrid gradient method (PDHGM): see Zhu and Chan (Reference Zhu and Chan2008), Pock, Cremers, Bischof and Chambolle (Reference Pock, Cremers, Bischof and Chambolle2009), Esser, Zhang and Chan (Reference Esser, Zhang and Chan2010) and Chambolle and Pock (Reference Chambolle and Pock2011, Reference Chambolle and Pock2016). We clearly observe the inverse scale-space nature of the Bregman iteration. The first iterate only contains features at a very coarse scale, and then ever more features at finer and finer scales are introduced throughout the course of the iteration.

Debiasing generalized eigenfunctions

We want to continue the analysis of the generalized eigenvalue problem introduced in Section 5.3. We have figured out that there is always a systematic bias of variational regularization methods for

$$\begin{eqnarray}F(Ku,f^{\unicode[STIX]{x1D6FF}})=\displaystyle \frac{1}{2}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2},\end{eqnarray}$$

i.e.  $\unicode[STIX]{x1D706}<1$ in (5.11) for $f^{\unicode[STIX]{x1D6FF}}=v_{\unicode[STIX]{x1D70E}}$ . Benning (Reference Benning2011) showed that this systematic bias can be corrected with the help of Bregman iterations in the case of scalar $\unicode[STIX]{x1D736}=\unicode[STIX]{x1D6FC}$ and $J$ with $\unicode[STIX]{x2202}J(cu,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}\unicode[STIX]{x2202}J(u)$ . Assume that $\unicode[STIX]{x1D6FC}$ is chosen such that $u^{k}=0$ for all $k<k^{\ast }-1$ , and $u^{k^{\ast }-1}=(1-\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D70E}^{-1}u_{\unicode[STIX]{x1D70E}}$ for some $k^{\ast }\in \mathbb{N}$ . Then we can easily conclude from (5.4) and (6.4) that $u^{k^{\ast }}$ has to satisfy

$$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D6FC}}\biggl(\unicode[STIX]{x1D706}-\displaystyle \frac{1-\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D70E}}\biggr)K^{\ast }Ku^{k^{\ast }}\in \unicode[STIX]{x2202}J(\unicode[STIX]{x1D706}u^{k^{\ast }}).\end{eqnarray}$$

We easily calculate that the above equation simplifies to the singular vector condition (5.14) for the choice $\unicode[STIX]{x1D706}=1/\unicode[STIX]{x1D70E}$ . Consequently, $u^{k^{\ast }}=R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC})=u_{\unicode[STIX]{x1D70E}}/\unicode[STIX]{x1D70E}$ , and we have corrected for the bias of the previous iterate.

These computations demonstrate that Bregman iterations correct for the systematic bias of variational regularization reconstructions of generalized singular vectors for one-homogeneous regularization functionals $J$ . However, the phenomenon is not limited to singular vectors. The following numerical toy example shows that the average reconstruction bias can be significantly reduced with the help of Bregman iterations. Assume the following set-up. Our forward model $K\in \mathbb{R}^{m\times n}$ , for $m=128$ and $n=512$ , is a matrix with its entries drawn randomly from ${\mathcal{N}}(0,1)$ . We define a sparse vector $u^{\dagger }\in \mathbb{R}^{n}$ with nine non-zero entries, drawn randomly from ${\mathcal{N}}(0,1)$ , and set $f=Ku^{\dagger }$ . Subsequently, we create one hundred instances of noisy data via $f_{j}^{\unicode[STIX]{x1D6FF}}:=f+n_{j}$ , for $n\in {\mathcal{N}}(0,0.5)$ and $j\in \{1,\ldots ,100\}$ . We now compute reconstructions for each of the one hundred instances with the following two regularization methods:

(6.8) $$\begin{eqnarray}\displaystyle R_{\text{Morozov}}(f_{j}^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FF}_{j})=\underset{u\in \mathbb{R}^{n}}{\operatorname{argmin}}\,\{\Vert u\Vert _{1}\quad \text{subject to }\Vert Ku-f_{j}^{\unicode[STIX]{x1D6FF}}\Vert _{2}\leq \unicode[STIX]{x1D6FF}_{j}\}, & & \displaystyle\end{eqnarray}$$

and

(6.9) $$\begin{eqnarray}\displaystyle & & \displaystyle R_{\text{Bregman}}(f_{j}^{\unicode[STIX]{x1D6FF}},\{u_{j}^{n}\}_{n=1}^{k-1},\unicode[STIX]{x1D6FC})\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in \mathbb{R}^{n}}{\operatorname{argmin}}\,\biggl\{\displaystyle \frac{1}{2}\biggl\|Ku-\biggl(kf^{\unicode[STIX]{x1D6FF}}-\mathop{\sum }_{n=1}^{k-1}Ku_{j}^{n}\biggr)\biggr\|_{2}^{2}+\unicode[STIX]{x1D6FC}\Vert u\Vert _{1}\biggr\},\end{eqnarray}$$

for

$$\begin{eqnarray}\displaystyle u_{j}^{\text{Morozov}} & \in & \displaystyle R_{\text{Morozov}}(f_{j}^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FF}_{j}),\nonumber\\ \displaystyle \unicode[STIX]{x1D6FF}_{j} & := & \displaystyle \displaystyle \frac{1}{2}\Vert Ku^{\dagger }-f_{j}^{\unicode[STIX]{x1D6FF}}\Vert _{2}^{2},\nonumber\\ \displaystyle u_{j}^{\text{Bregman}} & \in & \displaystyle R_{\text{Bregman}}(f_{j}^{\unicode[STIX]{x1D6FF}},\{u_{j}^{n}\}_{n=1}^{k^{\ast }-1},\unicode[STIX]{x1D6FC}),\nonumber\end{eqnarray}$$

and $k^{\ast }$ chosen according to Definition 6.1 for $\unicode[STIX]{x1D702}=1$ and $\unicode[STIX]{x1D6FF}=\unicode[STIX]{x1D6FF}_{j}$ . We then compute the average absolute bias and the standard deviation of the reconstructions, that is, we compute

(6.10) $$\begin{eqnarray}\displaystyle \biggl|u^{\dagger }-\displaystyle \frac{1}{100}\mathop{\sum }_{j=1}^{100}\hat{u} _{j}\biggr|\quad \text{and}\quad \sqrt{\displaystyle \frac{1}{99}\mathop{\sum }_{j=1}^{100}\biggl(\hat{u} _{j}-\displaystyle \frac{1}{100}\mathop{\sum }_{j=1}^{100}\hat{u} _{j}\biggr)^{2}} & & \displaystyle\end{eqnarray}$$

for $\hat{u} _{j}\in \{u_{j}^{\text{Morozov}},u_{j}^{\text{Bregman}}\}$ . Both average absolute bias and standard deviation are visualized for each of the $n=512$ coefficients in Figure 6.2. We clearly observe that with similar standard deviation, the average absolute bias is significantly reduced by the Bregman iteration in comparison to the Morozov regularization model.

Figure 6.2.

Equation (6.10) for the compressed sensing toy example in Section 6.1. The blue circles represent the standard deviation and average absolute bias values for all coefficients recovered with (6.8). The red circles show the same but for quantities for all coefficients recovered with (6.9). It becomes evident that for this example the average bias is significantly reduced, while the standard deviation of the reconstructed coefficients is comparable.

6.2 Linearized Bregman iteration

As the name suggests, the linearized Bregman iteration can be derived from Algorithm 1 by replacing the term $F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})$ with its linearization

$$\begin{eqnarray}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})\approx F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})+\langle \unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),Ku^{k}-Ku^{k-1}\rangle .\end{eqnarray}$$

Hence, if we replace $F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})$ in Algorithm 1 with this linearization multiplied by some constant $\unicode[STIX]{x1D70F}>0$ , we obtain

$$\begin{eqnarray}\displaystyle & & \displaystyle R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\{\unicode[STIX]{x1D70F}(F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})+\langle \unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),Ku^{k}-Ku^{k-1}\rangle )\nonumber\\ \displaystyle & & \displaystyle \qquad +\,D_{J(\cdot ,\unicode[STIX]{x1D6FC})}^{p^{k-1}}(u,u^{k-1})\},\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\{\unicode[STIX]{x1D70F}\langle \unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),Ku^{k}-Ku^{k-1}\rangle +D_{J(\cdot ,\unicode[STIX]{x1D6FC})}^{p^{k-1}}(u,u^{k-1})\},\nonumber\\ \displaystyle & & \displaystyle u^{k}\in R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & & \displaystyle p^{k}=p^{k-1}-\unicode[STIX]{x1D70F}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),\nonumber\end{eqnarray}$$

for $\unicode[STIX]{x1D736}=(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D6FC})$ and $v^{k-1}:=(u^{k},p^{k})$ for all $k\in \mathbb{N}$ . These equations are summarized in Algorithm 2.

The linearized Bregman iteration is a generalization of the Landweber regularization (Landweber Reference Landweber1951) for the choices

$$\begin{eqnarray}F(Ku,f^{\unicode[STIX]{x1D6FF}})=\displaystyle \frac{1}{2}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2}\quad \text{and}\quad J(u,\unicode[STIX]{x1D6FC})=\displaystyle \frac{\unicode[STIX]{x1D6FC}}{2}\Vert u\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2},\end{eqnarray}$$

for some signal domains $\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D6F4}$ . It is also a generalization of the mirror descent algorithm proposed in Nemirovskii and Yudin (Reference Nemirovskii and Yudin1983), where $J(\cdot ,\unicode[STIX]{x1D6FC})$ is a Legendre functional in the sense of Definition 5.11. This connection for convex, differentiable $F$ and strongly convex and differentiable $J(\cdot ,\unicode[STIX]{x1D6FC})$ was made in Beck and Teboulle (Reference Beck and Teboulle2003). The extension to subdifferentiable convex $J(\cdot ,\unicode[STIX]{x1D6FC})$ was first proposed in Darbon and Osher (Reference Darbon and Osher2007) and has since been studied extensively (Yin, Osher, Goldfarb and Darbon Reference Yin, Osher, Goldfarb and Darbon2008, Cai, Osher and Shen Reference Cai, Osher and Shen2009b , Cai, Osher and Shen Reference Cai, Osher and Shen2009a , Yin Reference Yin2010).

As in Remark 6.2, we can rewrite the dual update of the linearized Bregman iteration as

(6.11) $$\begin{eqnarray}\displaystyle p^{k}=p^{0}-\mathop{\sum }_{n=0}^{k-1}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f^{\unicode[STIX]{x1D6FF}}), & & \displaystyle\end{eqnarray}$$

and the primal update as

(6.12) $$\begin{eqnarray}\displaystyle & & \displaystyle R_{I}(f^{\unicode[STIX]{x1D6FF}},\{u^{n}\}_{n=0}^{k-1},p^{0},\unicode[STIX]{x1D6FC})\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\biggl\{J(u,\unicode[STIX]{x1D6FC})-\langle p^{0}-\mathop{\sum }_{n=0}^{k-1}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{n},f^{\unicode[STIX]{x1D6FF}}),u\rangle \biggr\}.\end{eqnarray}$$

In order to carry out a convergence analysis similar to the analysis for the standard Bregman iteration, we define the surrogate functional

(6.13) $$\begin{eqnarray}\displaystyle J_{\unicode[STIX]{x1D70F}}(u,\unicode[STIX]{x1D6FC}):=J(u,\unicode[STIX]{x1D6FC})-\unicode[STIX]{x1D70F}F(Ku,f^{\unicode[STIX]{x1D6FF}}). & & \displaystyle\end{eqnarray}$$

We further assume for the remainder of this section that $J$ and $\unicode[STIX]{x1D70F}$ are chosen such that $J_{\unicode[STIX]{x1D70F}}$ is convex. In practice, this requires strong convexity properties of $J$ , which can simply be established by adding a sufficiently strongly convex functional to the original choice of $J$ .

Example 6.9. Let

$$\begin{eqnarray}K\in {\mathcal{L}}(L^{2}(\unicode[STIX]{x1D6FA}),L^{2}(\unicode[STIX]{x1D6F4})),\quad F(Ku,f^{\unicode[STIX]{x1D6FF}})=\displaystyle \frac{1}{2}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2},\end{eqnarray}$$

for domains $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{n}$ and $\unicode[STIX]{x1D6F4}\subset \mathbb{R}^{m}$ , and let $J_{1}$ be a proper, lower semicontinuous and convex functional. Then the functional

$$\begin{eqnarray}J_{\unicode[STIX]{x1D70F}}(u,\unicode[STIX]{x1D6FC}):=J(u,\unicode[STIX]{x1D6FC})-\displaystyle \frac{\unicode[STIX]{x1D70F}}{2}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\unicode[STIX]{x1D6F4})}^{2}\end{eqnarray}$$

is convex for the choices

$$\begin{eqnarray}J(u,\unicode[STIX]{x1D6FC}):=\displaystyle \frac{1}{2}\Vert u\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2}+J_{1}(u)\quad \text{and}\quad \unicode[STIX]{x1D70F}<\displaystyle \frac{1}{\Vert K\Vert _{{\mathcal{L}}(L^{2}(\unicode[STIX]{x1D6FA}),L^{2}(\unicode[STIX]{x1D6F4}))}^{2}}.\end{eqnarray}$$

As for the Bregman iteration analysis, we start with a statement about the monotonic decrease of data fidelity.

Corollary 6.10 (monotonic decrease of data fidelity).

Suppose that $u^{0}$ satisfies $F(Ku^{0},f^{\unicode[STIX]{x1D6FF}})<\infty$ . Then the iterates of Algorithm 2 satisfy

(6.14) $$\begin{eqnarray}\displaystyle F(Ku^{k+1},f^{\unicode[STIX]{x1D6FF}})+\displaystyle \frac{1}{\unicode[STIX]{x1D70F}}D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k}}(u^{k+1},u^{k})\leq F(Ku^{k+1},f^{\unicode[STIX]{x1D6FF}}) & & \displaystyle\end{eqnarray}$$

and

$$\begin{eqnarray}\lim _{k\rightarrow \infty }D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k}}(u^{k+1},u^{k})=0,\end{eqnarray}$$

for $u^{k}\in R(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})$ and $q^{k}\in \unicode[STIX]{x2202}J_{\unicode[STIX]{x1D70F}}(u^{k},\unicode[STIX]{x1D6FC})$ .

Proof. First of all we emphasize that

$$\begin{eqnarray}\langle r^{k-1},u^{k}-u^{k-1}\rangle =\langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),u^{k}-u^{k-1}\rangle\end{eqnarray}$$

can be written as

$$\begin{eqnarray}\displaystyle & & \displaystyle \langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}),u^{k}-u^{k-1}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})-F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})-D_{F(K\cdot ,f^{\unicode[STIX]{x1D6FF}})}^{}(u^{k},u^{k-1}),\nonumber\end{eqnarray}$$

for all $k\in \mathbb{N}$ . Hence, the (primal) update of the linearized Bregman iteration can be rewritten as

$$\begin{eqnarray}\displaystyle & & \displaystyle R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})\nonumber\\ \displaystyle & & \displaystyle \quad =\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\{\unicode[STIX]{x1D70F}(F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})-F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}))+D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{k},u^{k-1})\},\nonumber\end{eqnarray}$$

for

$$\begin{eqnarray}q^{k-1}=p^{k-1}-\unicode[STIX]{x1D70F}K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})\in \unicode[STIX]{x2202}J_{\unicode[STIX]{x1D70F}}(u^{k-1},\unicode[STIX]{x1D6FC}),\quad p^{k-1}\in \unicode[STIX]{x2202}J(u^{k-1},\unicode[STIX]{x1D6FC}).\end{eqnarray}$$

Hence, we conclude

$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70F}(F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})-F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}))+D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{k},u^{k-1})\nonumber\\ \displaystyle & & \displaystyle \quad \leq \underbrace{\unicode[STIX]{x1D70F}(F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})-F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}}))}_{=0}+\underbrace{D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{k-1},u^{k-1})}_{=0},\nonumber\end{eqnarray}$$

and thus equation (6.14). In the same fashion as in the proof of Corollary 6.3, we further conclude $\lim _{k\rightarrow \infty }D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k}}(u^{k+1},u^{k})=0$ . ☐

As in the case of the standard Bregman iteration, the linearized Bregman iteration also satisfies Fejér monotonicity when the discrepancy principle is not violated.

Lemma 6.11 (Fejér monotonicity of Algorithm 2).

Let $f\in {\mathcal{R}}_{F}(K)$ , $u^{\dagger }\in {\mathcal{S}}(f,\unicode[STIX]{x1D6FC})$ and let $f^{\unicode[STIX]{x1D6FF}}\in {\mathcal{V}}$ with $F(f,f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}$ . We further assume that the iterates of Algorithm 2 satisfy Definition 6.1 for $\unicode[STIX]{x1D702}=1$ . Then the iterates also satisfy the strict Fejér monotonicity

$$\begin{eqnarray}D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k}}(u^{\dagger },u^{k})<D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{\dagger },u^{k-1}),\end{eqnarray}$$

for all $u^{k}\in R(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})$ and $q^{k}\in \unicode[STIX]{x2202}J_{\unicode[STIX]{x1D70F}}(u^{k},\unicode[STIX]{x1D6FC})$ , for all $k\leq k^{\ast }$ .

Proof. Through straightforward computations we obtain

$$\begin{eqnarray}\displaystyle & & \displaystyle D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k}}(u^{\dagger },u^{k})-D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{\dagger },u^{k-1})\nonumber\\ \displaystyle & & \displaystyle \quad =\underbrace{-D_{J_{\unicode[STIX]{x1D70F}}(\cdot ,\unicode[STIX]{x1D6FC})}^{q^{k-1}}(u^{k},u^{k-1})}_{{<}0}-\langle q^{k}-q^{k-1},u^{\dagger }-u^{k}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad \leq \langle p^{k-1}-p^{k}+\unicode[STIX]{x1D70F}K^{\ast }(\unicode[STIX]{x2202}_{x}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})-\unicode[STIX]{x2202}_{x}F(Ku^{k-1},f^{\unicode[STIX]{x1D6FF}})),u^{\dagger }-u^{k}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70F}\langle K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}),u^{\dagger }-u^{k}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad \leq \unicode[STIX]{x1D70F}(\unicode[STIX]{x1D6FF}-F(Ku^{k},f^{\unicode[STIX]{x1D6FF}}))\nonumber\\ \displaystyle & & \displaystyle \quad <0\nonumber\end{eqnarray}$$

for $k\leq k^{\ast }$ , where we have made use of the convexity of $F$ in its first argument, and $F(Ku^{\dagger },f^{\unicode[STIX]{x1D6FF}})\leq \unicode[STIX]{x1D6FF}$ . ☐

In analogy to Corollary 6.5, we can show the same result for the linearized Bregman iteration.

As in the case of Bregman iteration, Remark 6.6 follows from this result.

The following result guarantees convergence to a solution in ${\mathcal{S}}(f,\unicode[STIX]{x1D6FC})$ when $Ku^{k^{\ast }}=f$ is satisfied after a finite number $k^{\ast }$ of iterations of Algorithm 2.

Figure 6.3.

(a) Image $u^{\dagger }\in \mathbb{R}^{400\times 300}$ of Pixel the Gambian pouched rat, introduced in Figure 6.1(a). (b) The same degraded and noisy version $f^{\unicode[STIX]{x1D6FF}}\in \mathbb{R}^{400\times 300}$ together with the convolution kernel $h$ as shown in Figure 6.1(b). (c–h) Different iterates of Algorithm 2 for $F(Ku,f^{\unicode[STIX]{x1D6FF}})={\textstyle \frac{1}{2}}\Vert Ku-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\mathbb{R}^{2})}^{2}$ , $J(u,\unicode[STIX]{x1D6FC})={\textstyle \frac{1}{2}}\Vert u\Vert _{L^{2}(\mathbb{R}^{2})}^{2}+\unicode[STIX]{x1D6FC}\text{TV}(u)$ and $\unicode[STIX]{x1D6FC}=1/4$ . The 128th iterate visualized in (h) is the first that violates Definition 6.1, for $\unicode[STIX]{x1D6FF}=5.95$ .

We conclude this section with numerical results for the same deconvolution example introduced in Example 6.8. We observe that with the same choice of regularization parameter and the same initialization, Algorithm 2 requires more iterations in order to converge to a solution that violates the discrepancy principle with the same error bound. See Figure 6.3. On the other hand, the variational subproblems are computationally cheaper to solve compared to the standard Bregman iteration case, at least with the (accelerated) PDHGM used for this example.

6.3 Coupled and modified Bregman iterations

The Bregman iteration (as well as its linearized variant) leave some freedom for modifications. One obvious modification is the choice of the subgradient $p^{k-1}$ . The update from the optimality condition is of course the obvious one and particularly suitable for a convergence proof. However, one may use different ways to determine a subgradient $p^{k}$ from $u^{k}$ . As an example one may solve some variational problem

$$\begin{eqnarray}p^{k}\in \underset{p}{\operatorname{argmin}}\,\{H(p,p^{k-1})\mid p\in \unicode[STIX]{x2202}J(u^{k},\boldsymbol{\unicode[STIX]{x1D6FC}})\},\end{eqnarray}$$

with a convex functional $H$ . In the case of $\ell ^{1}$ -minimization one might choose $H(p,p^{k-1})=\Vert p\Vert _{2}$ , which yields the minimal subgradient, that is, again choosing sign $_{0}(p)$ (the single-valued version with sign $_{0}(0)=0$ ) in the case of a multivalued sign.

Another option for choosing subgradients has been investigated by Moeller et al. (Reference Moeller, Brinkmann, Burger and Seybold2014), when one solves joint reconstruction problems for multiple unknowns $u_{1},\ldots ,u_{M}$ . Moeller et al. (Reference Moeller, Brinkmann, Burger and Seybold2014) proposed and analysed a coupled Bregman iteration, which is based on choosing a new subgradient for the Bregman iteration in the $i$ th image $u_{i}$ from a linear combination of the subgradients in the other channels. In this way a joint subgradient for all the channels is approximated, which means a structural joint sparsity in the case of the $\ell ^{1}$ -norm or joint edge information in the total variation case. Rasch et al. (Reference Rasch, Brinkmann and Burger2018) have investigated an infimal convolution version of the coupled Bregman iteration for an application to PET-MR imaging.

7.1 Inverse scale space

For regularization functionals of the form $J(u,\unicode[STIX]{x1D736})=\unicode[STIX]{x1D6FC}J_{1}(u)$ we can write the dual Bregman iteration update as

$$\begin{eqnarray}\displaystyle \frac{p^{k}-p^{k-1}}{\unicode[STIX]{x1D6E5}t}=-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku^{k},f^{\unicode[STIX]{x1D6FF}})\end{eqnarray}$$

for $\unicode[STIX]{x1D6E5}t:=1/\unicode[STIX]{x1D6FC}$ and $p^{k}\in \unicode[STIX]{x2202}J_{1}(u^{k})$ , for all $k\in \mathbb{N}$ . Thus, taking the limit $\unicode[STIX]{x1D6FC}\rightarrow \infty$ , so that $\unicode[STIX]{x1D6E5}t\rightarrow 0$ , yields the following time-continuous formulation of the Bregman iteration, also known as the inverse scale space flow (Burger, Osher, Xu and Gilboa Reference Burger, Osher, Xu, Gilboa and Paragios2005, Burger et al. Reference Burger, Gilboa, Osher and Xu2006, Burger, Frick, Osher and Scherzer Reference Burger, Frick, Osher and Scherzer2007a ),

(7.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}p(t)=-K^{\ast }\unicode[STIX]{x2202}_{x}F(Ku(t),f^{\unicode[STIX]{x1D6FF}}), & & \displaystyle\end{eqnarray}$$

for $p(t)\in \unicode[STIX]{x2202}J_{1}(u(t))$ .

For many typical choices of regularization functionals $J_{1}$ , it is difficult to numerically compute solutions of (7.1), with the $\ell ^{1}$ -norm and in general any polyhedral regularization functional being the exception (Burger, Moeller, Benning and Osher Reference Burger, Moeller, Benning and Osher2013c , Moeller Reference Moeller2012, Moeller and Burger Reference Moeller and Burger2013). Nevertheless, (7.1) is very useful for studying theoretical properties of iterative regularizations in the limiting case.

Unsurprisingly, it is straightforward to carry out an eigenanalysis similar to the one discussed in Sections 5.3 and 6.1 for the regularization operator $R(f^{\unicode[STIX]{x1D6FF}},t)=u(t)$ with $u(t)$ satisfying (7.1) for $F(Ku,f^{\unicode[STIX]{x1D6FF}})=G(Ku-f^{\unicode[STIX]{x1D6FF}})$ . The following result is a generalization of Benning and Burger (Reference Benning and Burger2013, Theorem 9).

Theorem 7.1. Let $(u_{\unicode[STIX]{x1D70E}},v_{\unicode[STIX]{x1D70E}})$ be a pair of generalized singular vectors with singular value $\unicode[STIX]{x1D70E}$ , $f=v_{\unicode[STIX]{x1D70E}}$ and suppose $J_{1}$ is (absolutely) one-homogeneous, that is, $J_{1}(cu)=|c|J_{1}(u)$ for all $c\in \mathbb{R}$ . Then $0\in R(v_{\unicode[STIX]{x1D70E}},t)$ for $0\leq t<t_{\ast }$ and

$$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D70E}}u_{\unicode[STIX]{x1D70E}}\in R(v_{\unicode[STIX]{x1D70E}},t)\end{eqnarray}$$

for $t\geq t_{\ast }=1$ .

Proof. First we verify $0\in R(v_{\unicode[STIX]{x1D70E}},t)$ for $0\in [0,t_{\ast })$ . From (5.14) and the absolute one-homogeneity of $J_{1}$ , we observe $J_{1}(u_{\unicode[STIX]{x1D70E}})=\langle G^{\prime }(v_{\unicode[STIX]{x1D70E}}),Ku_{\unicode[STIX]{x1D70E}}\rangle$ . We further see from the definition of the subdifferential that $t\leq 1=\langle G^{\prime }(v_{\unicode[STIX]{x1D70E}}),Ku_{\unicode[STIX]{x1D70E}}\rangle /J_{1}(u_{\unicode[STIX]{x1D70E}})$ implies $p(t):=tK^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}})\in \unicode[STIX]{x2202}J_{1}(0)$ . Since $\unicode[STIX]{x2202}_{t}p(t)=K^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}})$ and $p(0)=0$ , we have shown that $u(t)=0$ is a solution of (7.1).

For $t\geq t_{\ast }$ a continuous extension of $p(t)$ is

$$\begin{eqnarray}p(t)=p(t_{\ast })+(t_{\ast }-t)K^{\ast }G^{\prime }(Ku(t)-v_{\unicode[STIX]{x1D70E}}).\end{eqnarray}$$

We immediately see that $u(t)=u_{\unicode[STIX]{x1D70E}}/\unicode[STIX]{x1D70E}$ is a solution for $t\geq t_{\ast }$ , since $p(t_{\ast })=t_{\ast }K^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}})\in \unicode[STIX]{x2202}J_{1}(u_{\unicode[STIX]{x1D70E}}/\unicode[STIX]{x1D70E})$ and $\unicode[STIX]{x2202}_{t}p(t)=0$ . ☐

Hence, the inverse scale space reconstruction also has no bias (for input data $v_{\unicode[STIX]{x1D70E}}$ satisfying (5.14)), compared to the variational regularization method.

A similar result can be derived even in the case of noisy data $f^{\unicode[STIX]{x1D6FF}}=v_{\unicode[STIX]{x1D70E}}+n$ , where $n$ is an error term that satisfies the specific source condition

$$\begin{eqnarray}\unicode[STIX]{x1D707}K^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}})+\unicode[STIX]{x1D702}K^{\ast }n\in \unicode[STIX]{x2202}J(\unicode[STIX]{x1D70E}u_{\unicode[STIX]{x1D70E}}),\end{eqnarray}$$

for constants $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D702}$ . For more details we refer to Benning and Burger (Reference Benning and Burger2013, Theorem 10).

In the following we briefly want to discuss reconstruction guarantees for linear combinations of multiple singular vectors. More precisely, we ask when we can guarantee

$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FE}_{j}}{\unicode[STIX]{x1D70E}_{j}}u_{\unicode[STIX]{x1D70E}_{j}}\in R\biggl(\mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D6FE}_{j}v_{\unicode[STIX]{x1D70E}_{j}},t\biggr),\end{eqnarray}$$

for coefficients $\{\unicode[STIX]{x1D6FE}_{j}\}_{j\in \mathbb{N}}$ . Due to the nonlinearity of $J_{1}$ , in general there is no such decomposition. If we restrict ourselves to the following two conditions, however, such a result can be guaranteed (Schmidt, Benning and Schönlieb Reference Schmidt, Benning and Schönlieb2018, Theorem 3.14). The first condition is $K$ -orthogonality of the singular vectors, that is,

(OC) $$\begin{eqnarray}\langle Ku_{\unicode[STIX]{x1D70E}_{i}},Ku_{\unicode[STIX]{x1D70E}_{j}}\rangle =\left\{\begin{array}{@{}ll@{}}1\quad & i=j,\\ 0\quad & i\neq j,\end{array}\right.\end{eqnarray}$$

for $i,j\in \{1,\ldots ,n\}$ . The second condition is the so-called (SUB0)-condition, which reads as follows.

Definition 7.2 (Schmidt et al. Reference Schmidt, Benning and Schönlieb2018, Definition 3.1).

Suppose that $(u_{\unicode[STIX]{x1D70E}_{1}},u_{\unicode[STIX]{x1D70E}_{2}},\ldots ,u_{\unicode[STIX]{x1D70E}_{n}})$ are an ordered set of primal singular vectors of $J_{1}$ with corresponding dual singular vectors $(v_{\unicode[STIX]{x1D70E}_{1}},v_{\unicode[STIX]{x1D70E}_{2}},\ldots ,v_{\unicode[STIX]{x1D70E}_{n}})$ and singular values $(\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{n})$ . Then the singular vectors satisfy the (SUB0) condition if

(SUB0) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{j=1}^{k}K^{\ast }G^{\prime }(v_{\unicode[STIX]{x1D70E}_{j}})\in \unicode[STIX]{x2202}J_{1}(0), & & \displaystyle\end{eqnarray}$$

for all $k\in \{1,\ldots ,n\}$ .

Given (OC) and (SUB0), we can guarantee the following decomposition result, which is a direct generalization of Schmidt et al. (Reference Schmidt, Benning and Schönlieb2018, Theorem 3.14).

Theorem 7.3. Let $(u_{\unicode[STIX]{x1D70E}_{1}},u_{\unicode[STIX]{x1D70E}_{2}},\ldots ,u_{\unicode[STIX]{x1D70E}_{n}})$ , $(v_{\unicode[STIX]{x1D70E}_{1}},v_{\unicode[STIX]{x1D70E}_{2}},\ldots ,v_{\unicode[STIX]{x1D70E}_{n}})$ be a system of ordered singular vectors, for singular values $(\unicode[STIX]{x1D70E}_{1},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{n})$ , for which the $v_{j}$ are normalized and (OC) and (SUB0) are satisfied. Then, for data $f=\sum _{j=1}^{n}\unicode[STIX]{x1D6FE}_{j}v_{\unicode[STIX]{x1D70E}_{j}}$ with positive coefficients $(\unicode[STIX]{x1D6FE}_{1},\ldots ,\unicode[STIX]{x1D6FE}_{n})$ we have $u(t)\in R(f,t)$ , with

$$\begin{eqnarray}u(t)=\left\{\begin{array}{@{}ll@{}}0\quad & 0\leq t\leq t_{1},\\ \mathop{\sum }_{j=1}^{k}\unicode[STIX]{x1D6FE}_{j}{\unicode[STIX]{x1D70E}_{j}}^{-1}\,u_{\unicode[STIX]{x1D70E}_{j}}\quad & t_{k}\leq t<t_{k+1},\,\text{for all }k=1,\ldots ,n-1,\\ \mathop{\sum }_{j=1}^{n}\unicode[STIX]{x1D6FE}_{j}{\unicode[STIX]{x1D70E}_{j}}^{-1}\,u_{\unicode[STIX]{x1D70E}_{j}}\quad & t_{n}\leq t,\end{array}\right.\end{eqnarray}$$

where $t_{k}=\unicode[STIX]{x1D6FE}_{k}$ and $t_{k}<t_{k+1}$ for all $k\in \{1,\ldots ,n\}$ .

We refer to Benning and Burger (Reference Benning and Burger2013) for more information on individual generalized singular vectors and the inverse scale space flow. For more theoretical results and analytical as well as numerical examples of ordered sets of singular vectors that satisfy (OC) and (SUB0), we refer to Schmidt et al. (Reference Schmidt, Benning and Schönlieb2018).

7.2 Two-step debiasing

While Bregman iterations and inverse scale space methods perform debiasing in an iterative fashion (and effectively change the variational model), one may also consider two-step procedures that first solve the original variational model and then perform a second step to reduce the bias (Deledalle, Papadakis and Salmon Reference Deledalle, Papadakis, Salmon and Aujol2015, Deledalle, Papadakis, Salmon and Vaiter Reference Deledalle, Papadakis, Salmon and Vaiter2017). The first and simplest case where this idea was brought up is regularization with the $\ell ^{1}$ -norm, where a so-called refitting strategy (Lederer Reference Lederer2013) is quite natural. After the variational problem

(7.2) $$\begin{eqnarray}u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}\in \operatorname{argmin}F(Ku,f^{\unicode[STIX]{x1D6FF}})+\unicode[STIX]{x1D6FC}\Vert u\Vert _{\ell ^{1}}\end{eqnarray}$$

is solved, the second step simply consists in minimizing $F(Ku,f^{\unicode[STIX]{x1D6FF}})$ over the set of all $u$ sharing the support of $u_{\unicode[STIX]{x1D6FC}}^{\unicode[STIX]{x1D6FF}}$ . Since this procedure throws away information about the sign of the entries of $u$ , one can further improve to define the regularization operator via

(7.3) $$\begin{eqnarray}R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC})=\operatorname{argmin}\{F(Ku,f^{\unicode[STIX]{x1D6FF}})\mid \text{sign}_{0}(u_{i})=\text{sign}_{0}((u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}})_{i}),~\text{for all}~i\}.\end{eqnarray}$$

where sign $_{0}(u_{i})$ is the single-valued sign (i.e. zero for $u_{i}=0$ ). Since the sign corresponds to a subgradient of the $\ell ^{1}$ -norm, we can reinterpret the debiased regularization operator in a variational way: we minimize the fidelity subject to the constraint of $u$ sharing a subgradient with $u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}$ . This is a key observation towards a generalization for arbitrary convex regularizations, as noted by Brinkmann et al. (Reference Brinkmann, Burger, Rasch and Sutour2017). The general debiasing problem can be rephrased as a two-step procedure,

(7.4) $$\begin{eqnarray}u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}\in \underset{u}{\operatorname{argmin}}\,F(Ku,f^{\unicode[STIX]{x1D6FF}})+\unicode[STIX]{x1D6FC}J(u),\end{eqnarray}$$

followed by

(7.5) $$\begin{eqnarray}R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC})=\underset{u}{\operatorname{argmin}}\,\{F(Ku,f^{\unicode[STIX]{x1D6FF}})\mid p\in \unicode[STIX]{x2202}J(u)\cap \unicode[STIX]{x2202}J(u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}})\}.\end{eqnarray}$$

For computational purposes the arbitrary choice of the subgradient $p\in \unicode[STIX]{x2202}J(u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}})$ is not suitable, but we can indeed use the subgradient from the first step. Noting that for differentiable fidelities the optimality condition reads

(7.6) $$\begin{eqnarray}p_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}=-\displaystyle \frac{1}{\unicode[STIX]{x1D6FC}}K^{\ast }\unicode[STIX]{x2202}F(Ku,f^{\unicode[STIX]{x1D6FF}})\in \unicode[STIX]{x2202}J(u_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}),\end{eqnarray}$$

we can use the debiasing procedure

(7.7) $$\begin{eqnarray}R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC})=\underset{u}{\operatorname{argmin}}\,\{F(Ku,f^{\unicode[STIX]{x1D6FF}})\mid p_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}\in \unicode[STIX]{x2202}J(u)\}.\end{eqnarray}$$

The condition $p_{\unicode[STIX]{x1D6FF}}^{\unicode[STIX]{x1D6FC}}\in \unicode[STIX]{x2202}J(u)$ can be reformulated as a vanishing Bregman distance between $u$ and $u_{\unicode[STIX]{x1D6FC}}^{\unicode[STIX]{x1D6FF}}$ ; thus we observe some relations to the Bregman iteration. The second step can be interpreted as a Bregman iteration step in the limit of the regularization parameter to infinity. We refer to Brinkmann et al. (Reference Brinkmann, Burger, Rasch and Sutour2017) for a detailed analysis of this debiasing approach.

The effect of the debiasing is illustrated for the simple case of total variation denoising, that is, the solution of

(7.8) $$\begin{eqnarray}R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC})=\underset{u\in BV(\unicode[STIX]{x1D6FA})}{\operatorname{argmin}}\,\biggl(\displaystyle \frac{1}{2}\Vert u-f^{\unicode[STIX]{x1D6FF}}\Vert _{L^{2}(\unicode[STIX]{x1D6FA})}^{2}+\unicode[STIX]{x1D6FC}|u|_{BV}\biggr).\end{eqnarray}$$

Figure 7.1 compares the solution of the variational problem in (c) with the one obtained in the two-step debiasing procedure (d) and the Bregman iteration (e). Both methods reduce the contrast loss of the TV regularization (which is difficult to see in the image, but becomes more apparent in the small background buildings). Overall, however, the Bregman iteration seems to restore more of the small details such as the grass structure. Figure 7.2 demonstrates the debiasing effect for increasing regularization parameters, where the variational model destroys ever more detail. In particular, for larger $\unicode[STIX]{x1D6FC}$ one observes the effect of restoring smaller structures apparently contained in the subgradient but not the primal variable of the variational model.

Figure 7.1.

Cameraman (256 $\times$ 256): comparison of TV denoising for $\unicode[STIX]{x1D6FC}=0.1$ , with the two-step debiasing, and Bregman iterations ( $\unicode[STIX]{x1D6FC}=0.5$ and 7 Bregman iterations).

Figure 7.2.

Cameraman: TV denoising and debiasing for different values of the regularization parameter.

7.3 Nonlinear spectral transform

The iterative regularization methods presented in Section 6 can easily be extended to nonlinear spectral decomposition methods via the following trivial observation. Every iterate $u^{k}\in R_{I}(f^{\unicode[STIX]{x1D6FF}},v^{k-1},\unicode[STIX]{x1D736})$ can be represented as the sum of the differences of two subsequent iterates, that is,

$$\begin{eqnarray}u^{k}=u^{0}+\mathop{\sum }_{n=1}^{k}u^{n}-u^{n-1}.\end{eqnarray}$$

If we define $\unicode[STIX]{x1D711}^{0}:=u^{0}$ and $\unicode[STIX]{x1D711}^{n}:=u^{n}-u^{n-1}$ for $n>1$ , and equip the sum with coefficients $\{c^{n}\}_{n=0}^{k}$ , we can write $u^{k}$ as

$$\begin{eqnarray}u^{k}=\mathop{\sum }_{n=0}^{k}c^{n}\unicode[STIX]{x1D711}^{n}.\end{eqnarray}$$

In the following we are going to motivate why such a decomposition is useful for localizing individual scales if the underlying regularization functional is (absolutely) one-homogeneous and where we have a scalar parameter $\unicode[STIX]{x1D6FC}$ . Following up on the bias correction example for generalized singular vectors in Section 6.1, we know that for the Bregman iteration $R_{I}(f,v^{k-1},\unicode[STIX]{x1D6FC})$ with $f=v_{\unicode[STIX]{x1D70E}}$ we observe

$$\begin{eqnarray}u^{k}=\left\{\begin{array}{@{}ll@{}}0\quad & k<k^{\ast },\\ (k^{\ast }-\unicode[STIX]{x1D6FC})\unicode[STIX]{x1D70E}^{-1}\,u_{\unicode[STIX]{x1D70E}}\quad & k=k^{\ast },\\ \unicode[STIX]{x1D70E}^{-1}\,u_{\unicode[STIX]{x1D70E}}\quad & k\geq k^{\ast }+1.\end{array}\right.\end{eqnarray}$$

Replacing $f=v_{\unicode[STIX]{x1D70E}}$ with $f=\unicode[STIX]{x1D70E}v_{\unicode[STIX]{x1D70E}}=Ku_{\unicode[STIX]{x1D70E}}$ therefore yields

$$\begin{eqnarray}u^{k}=\left\{\begin{array}{@{}ll@{}}0\quad & k<k^{\ast },\\ (k^{\ast }-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70E}^{-1})\,u_{\unicode[STIX]{x1D70E}}\quad & k=k^{\ast },\\ u_{\unicode[STIX]{x1D70E}}\quad & k\geq k^{\ast }+1,\end{array}\right.\end{eqnarray}$$

and consequently we observe

$$\begin{eqnarray}\unicode[STIX]{x1D711}^{n}=\left\{\begin{array}{@{}ll@{}}0\quad & n\not \in \{k^{\ast },k^{\ast }+1\},\\ (k^{\ast }-\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70E}^{-1})\,u_{\unicode[STIX]{x1D70E}}\quad & n=k^{\ast },\\ (1+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70E}^{-1}-k^{\ast })\,u_{\unicode[STIX]{x1D70E}}\quad & n=k^{\ast }+1.\end{array}\right.\end{eqnarray}$$

The last equation implies that if the input datum is given in terms of the forward model applied to a (primal) singular vector, this primal singular vector is localized in only two components $\unicode[STIX]{x1D711}^{k^{\ast }}$ and $\unicode[STIX]{x1D711}^{k^{\ast }+1}$ . The index $k^{\ast }$ depends on the choice of $\unicode[STIX]{x1D6FC}$ and on the singular value $\unicode[STIX]{x1D70E}$ . Hence, singular vectors with different scales, or different values of $\unicode[STIX]{x1D70E}$ , will be localized in $\unicode[STIX]{x1D711}^{\hat{k}}$ and $\unicode[STIX]{x1D711}^{\hat{k}+1}$ for $\hat{k}\neq k^{\ast }$ . This is visualized in Figures 7.3 and 7.4. It is therefore fair to call $\{\unicode[STIX]{x1D711}^{n}\}_{n=1}^{k}$ a spectrum and the individual $\unicode[STIX]{x1D711}^{n}$ , for $n\in \{1,\ldots ,k\}$ , the spectral components.

Figure 7.3.

Two singular vectors of $J=\text{TV}_{\ast }$ with different $\unicode[STIX]{x1D70E}$ -values, and excerpts of their corresponding (analytically computed) spectra, for $\unicode[STIX]{x1D6FC}=1.24$ . We clearly observe that both vectors are located at different positions of the spectrum. Hence, both singular vectors could be isolated from a sum of the two by applying a band-pass filter to the spectrum.

Figure 7.4.

Spectral decomposition of the image of a bee. From Benning et al. (Reference Benning, Möller, Nossek, Burger, Cremers, Gilboa, Schönlieb and Lauze2017d ).

Consequently, the operator ${\mathcal{S}}:{\mathcal{V}}\times {\mathcal{U}}^{k}\times \mathbb{R}^{k}\times A\rightarrow {\mathcal{U}}$ with

$$\begin{eqnarray}{\mathcal{S}}(f,(u_{n})_{n=0}^{k},(c_{n})_{n=0}^{k},\unicode[STIX]{x1D6FC}):=\mathop{\sum }_{n=0}^{k}c_{n}\unicode[STIX]{x1D711}^{n}\quad \text{with }\unicode[STIX]{x1D711}^{n}:=\left\{\begin{array}{@{}ll@{}}u^{n}-u^{n-1}\quad & n>1,\\ u^{0}\quad & n=1,\end{array}\right.\end{eqnarray}$$

for $u^{n}\in R_{I}(f,v^{n-1},\unicode[STIX]{x1D6FC})$ can be seen as a spectral transform of the input signal $f^{\unicode[STIX]{x1D6FF}}$ . For $K:\text{BV}(\unicode[STIX]{x1D6FA})\rightarrow L^{2}(\unicode[STIX]{x1D6FA})$ this type of spectral transform is a discretization of the inverse-scale-space based spectral transform defined in Burger, Eckardt, Gilboa and Moeller (Reference Burger, Eckardt, Gilboa, Moeller and Aujol2015a ) and Burger et al. (Reference Burger, Gilboa, Moeller, Eckardt and Cremers2016a ). For $K:\text{BV}(\unicode[STIX]{x1D6FA})\rightarrow L^{2}(\unicode[STIX]{x1D6FA})$ and $J(u,\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D6FC}\text{TV}(u)$ , the idea of generalized spectral transforms goes back to Gilboa (Reference Gilboa2014a , Reference Gilboa2014b ). For a detailed overview of this form of nonlinear spectral transform we refer to Gilboa, Moeller and Burger (Reference Gilboa, Moeller and Burger2016). Another interesting recent extension is the spectral transform in the context of image segmentation (Zeune et al. Reference Zeune, van Dalum, Terstappen, van Gils and Brune2017).

8.1 Velocity-encoded magnetic resonance imaging

Magnetic resonance imaging (MRI) is an imaging technique that allows us to visualize the chemical composition of humans/animals or materials. MRI scanners utilize strong magnetic fields and radio waves to excite subatomic particles, such as protons, that subsequently emit radio frequency signals which can be measured with the radio frequency coils that initially excited those radio waves: see for example Callaghan (Reference Callaghan1993). MRI is often used to measure contrast in tissue. However, due to shear, endless possibilities of radio-frequency pulse sequence design, and programming of the gradient coils, MRI is a versatile imaging tool with capabilities beyond imaging contrast in tissue. A potential, more sophisticated application is phase-encoded magnetic resonance velocity imaging, which in medical imaging is used to study the distribution and variation in blood flow (Gatehouse et al. Reference Gatehouse, Keegan, Crowe, Masood, Mohiaddin, Kreitner and Firmin2005). In the physical sciences it is used to study the rheology of complex fluids (Callaghan Reference Callaghan1999), liquids and gases flowing through packed beds (Sederman, Johns, Alexander and Gladden Reference Sederman, Johns, Alexander and Gladden1998, Holland et al. Reference Holland, Malioutov, Blake, Sederman and Gladden2010), granular flows (Holland et al. Reference Holland, Müller, Dennis, Gladden and Sederman2008) and multiphase turbulence flows (Tayler, Holland, Sederman and Gladden Reference Tayler, Holland, Sederman and Gladden2012). The main advantage of MRI over other methods when it comes to studying flow is that it is possible to image flows non-invasively. However, the main drawback of the technique is the acquisition time of the measurement.

Figure 8.1.

A simulated spiral on a Cartesian grid. From Benning et al. (Reference Benning, Gladden, Holland, Schönlieb and Valkonen2014).

Figure 8.2.

Magnitude images of the velocity dataset used in Benning et al. (Reference Benning, Gladden, Holland, Schönlieb and Valkonen2014), courtesy of Andrew J. Sederman. (a) Magnitude image derived from applying the inverse of the Fourier transform to the fully sampled Fourier data of the velocity-encoded MRI measurement and subsequently taking the modulus. (b) Magnitude image obtained if we set to zero all Fourier samples that are not part of the spiral visualized in Figure 8.1, and subsequently proceed as for the fully sampled data. (c) Magnitude reconstructions from the $\text{TGV}_{\unicode[STIX]{x1D6FD}}^{2}$ -based variational regularization reconstruction (8.1).

Figure 8.3.

The different velocity reconstructions corresponding to the magnitude reconstructions in Figure 8.2.

Lustig, Donoho and Pauly (Reference Lustig, Donoho and Pauly2007) exploited the idea of sub-sampling in the spatial data domain to overcome this limitation and to speed up the MRI acquisition process. Due to fewer measurements being taken, compared to the number of unknowns to be recovered, some form of regularization needs to be integrated into the reconstruction process. Sparsity-promoting variational regularization methods are suitable candidates and, most prominently, total variation regularization has been successfully deployed to increase the temporal resolution of MRI acquisitions. Since measurement noise in MRI data can be modelled as being normally distributed, a standard variational regularization approach is

(8.1) $$\begin{eqnarray}\displaystyle R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D736})=\underset{u\in {\mathcal{U}}}{\operatorname{argmin}}\,\biggl\{\displaystyle \frac{1}{2}\Vert {\mathcal{F}}u-f^{\unicode[STIX]{x1D6FF}}\Vert _{2}^{2}+J(u,\unicode[STIX]{x1D736})\biggr\}, & & \displaystyle\end{eqnarray}$$

where ${\mathcal{F}}$ is the operator

$$\begin{eqnarray}({\mathcal{F}}u)(t^{k}):=(2\unicode[STIX]{x1D70B})^{-n/2}\int _{\mathbb{R}^{n}}u(x)\exp \biggl(-\text{i}\int _{t^{k-1}}^{t^{k}}x(t)\cdot g(t)\,\text{d}t\biggr)\,\text{d}x,\end{eqnarray}$$

and $n\in \{2,3\}$ denotes the dimension of the signal and $g:[0,T]\rightarrow \mathbb{R}^{n}$ represents the function that controls the gradient coils of the MRI machine. We observe that ${\mathcal{F}}$ is almost identical to the Fourier transform sampled at discrete locations, if we can approximate

$$\begin{eqnarray}\int _{t^{k-1}}^{t^{k}}x(t)\cdot g(t)\,\text{d}t\approx x\cdot \int _{t^{k-1}}^{t^{k}}g(t)\,\text{d}t.\end{eqnarray}$$

This can be achieved by adequate programming of the gradient coils. However, $\int _{t^{k-1}}^{t^{k}}x(t)\cdot g(t)\,\text{d}t$ can be approximated more generally via the Taylor series

$$\begin{eqnarray}\int _{t^{k-1}}^{t^{k}}x(t)\cdot g(t)\,\text{d}t\approx \mathop{\sum }_{r=0}^{\infty }\displaystyle \frac{x^{(r)}(t^{k-1})}{r!}\cdot \int _{t^{k-1}}^{t^{k}}g(t)t^{r}\,\text{d}t,\end{eqnarray}$$

and with clever programming of $g$ , other moments such as velocity or acceleration can be so encoded. In the following, we assume that the radio-frequency pulse sequence and the gradient coils are programmed such that we first encode the velocity in the $z$ -direction, that is, for $x=(x_{1},x_{2},x_{3})$ and $g(t)=(g_{1}(t),g_{2}(t),g_{3}(t))$ we have

$$\begin{eqnarray}\int _{0}^{t_{0}}x(t)\cdot g(t)\,\text{d}t\approx \underbrace{x_{3}^{\prime }(0)}_{=:v_{z}}\int _{0}^{t_{0}}g_{3}(t)\,t\,\text{d}t,\end{eqnarray}$$

in the interval $[0,t_{0}]$ , and then perform the spatial encoding such that

$$\begin{eqnarray}\int _{t_{k-1}}^{t^{k}}x(t)\cdot g(t)\,\text{d}t\approx \left(\begin{array}{@{}c@{}}x_{1}(t_{k-1})\\ x_{2}(t_{k-1})\end{array}\right)\cdot \int _{t^{k-1}}^{t^{k}}\left(\begin{array}{@{}c@{}}g_{1}(t)\\ g_{2}(t)\end{array}\right)\,\text{d}t\end{eqnarray}$$

holds true for $t_{0}<t_{1}<\cdots <t_{m}=T$ . Then, with $x=(x_{1}(t_{k-1}),x_{2}(t_{k-1}))$ and $g=(g_{1},g_{2})$ as an abuse of notation, ${\mathcal{F}}$ reads as

(8.2) $$\begin{eqnarray}\displaystyle ({\mathcal{F}}(u,v_{z}))(t^{k})=\displaystyle \frac{1}{2\unicode[STIX]{x1D70B}}\int _{\mathbb{R}^{2}}u(x)\exp (-\text{i}\unicode[STIX]{x1D70E}v_{z}(x))\exp \biggl(-\text{i}x\cdot \int _{t^{k-1}}^{t^{k}}g(t)\,\text{d}t\biggr)\,\text{d}x, & & \displaystyle\end{eqnarray}$$

for some constant $\unicode[STIX]{x1D70E}$ . In order to avoid non-linearity of the forward model, we couple $u$ and $v_{z}$ by simply defining $w:=u\exp (-\text{i}\unicode[STIX]{x1D70E}v_{z})$ . Then the forward model ${\mathcal{F}}$ simply reduces to the (sub-sampled) Fourier transform.

Benning et al. (Reference Benning, Gladden, Holland, Schönlieb and Valkonen2014) have investigated three choices for regularization functionals: assuming $\unicode[STIX]{x1D736}=(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ , we have

(8.3) $$\begin{eqnarray}\displaystyle J(w,\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=\unicode[STIX]{x1D6FC}\left\{\begin{array}{@{}l@{}}\text{TV}(w),\quad \\ \text{TGV}_{\unicode[STIX]{x1D6FD}}^{2}(w),\quad \\ \mathop{\sum }_{j=1}^{\infty }|\langle w,\unicode[STIX]{x1D711}_{j}\rangle |.\quad \end{array}\right. & & \displaystyle\end{eqnarray}$$

Here $\{\unicode[STIX]{x1D711}_{j}\}_{j\in \mathbb{Z}}$ denotes a wavelet basis, and $\text{TGV}_{\unicode[STIX]{x1D6FD}}^{2}$ is the second-order total generalized variation in the sense of Bredies and Valkonen (Reference Bredies and Valkonen2011, Theorem 3.1). In Figure 8.2 we see computational solutions of (8.1) for the choice $J(w,\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=\text{TGV}_{\unicode[STIX]{x1D6FD}}^{2}(w)$ , a spiral sub-sampling strategy on a Cartesian grid (see Figure 8.1), and the parameter choices $\unicode[STIX]{x1D6FC}=0.1$ and $\unicode[STIX]{x1D6FD}=3$ . These results have again been computed with the PDHGM. Subsequently, $v_{z}$ was extracted as the principal value of the reconstruction $w\in R(f^{\unicode[STIX]{x1D6FF}},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})$ . The reconstructed $z$ -velocity $v_{z}$ is subsequently unwrapped by solving the linear system

$$\begin{eqnarray}\unicode[STIX]{x1D6E5}\hat{v}_{z}=\cos (v_{z})\unicode[STIX]{x1D6E5}\sin (v_{z})-\sin (v_{z})\unicode[STIX]{x1D6E5}\cos (v_{z})\end{eqnarray}$$

for $\hat{v}_{z}$ . Here $\unicode[STIX]{x1D6E5}$ denotes the Laplace operator. The unwrapped reconstructed velocity $\hat{v}_{z}$ is visualized in Figure 8.3.