Augmented NETT Regularization of Inverse Problems

TL;DR

aNETT introduces a novel data-driven regularization framework for inverse problems, combining encoder-decoder networks with Tikhonov regularization, offering stability, convergence guarantees, and efficient reconstruction especially for complex models.

Contribution

It presents a new regularization method using neural networks that does not require repeated forward model evaluations and can adapt to different sampling rates without retraining.

Findings

Achieves state-of-the-art results in sparse-view and low-dose CT reconstruction.

Provides rigorous convergence analysis including stability and rates.

Enables efficient inverse problem solving with expensive forward models.

Abstract

We propose aNETT (augmented NETwork Tikhonov) regularization as a novel data-driven reconstruction framework for solving inverse problems. An encoder-decoder type network defines a regularizer consisting of a penalty term that enforces regularity in the encoder domain, augmented by a penalty that penalizes the distance to the data manifold. We present a rigorous convergence analysis including stability estimates and convergence rates. For that purpose, we prove the coercivity of the regularizer used without requiring explicit coercivity assumptions for the networks involved. We propose a possible realization together with a network architecture and a modular training strategy. Applications to sparse-view and low-dose CT show that aNETT achieves results comparable to state-of-the-art deep-learning-based reconstruction methods. Unlike learned iterative methods, aNETT does not require repeated application of the forward and adjoint models, which enables the use of aNETT for inverse problems with numerically expensive forward models. Furthermore, we show that aNETT trained on coarsely sampled data can leverage an increased sampling rate without the need for retraining.

Tables

Table 1. Table 3.1 : Parameter specifications for proposed aNETT functional and it numerical minimization.

	$\alpha$	$c$	$\gamma$	$N_{\mathrm{iter}}$	$N_{\phi }$	noise model	$𝒟$
Sparse view	${10}^{-4}$	${10}^2$	$5\cdot {10}^{-1}$	50	$40$	Gaussian, $\sigma =0.02$	${𝒟}_{{\parallel \cdot \parallel }^2}$
Low dose	$5\cdot {10}^{-3}$	${10}^2$	${10}^{-3}$	$20$	1138	Poisson, $p={10}^4$	${𝒟}_{\mathrm{KL}}$
Universality	${10}^{-4}$	${10}^2$	$5\cdot {10}^{-1}$	$50$	160	Gaussian, $\sigma =0.02$	${𝒟}_{{\parallel \cdot \parallel }^2}$

Table 2. Table 4.1 : Overview of metric results evaluated on the test-set. The values shown are the average of the PSNR ± plus-or-minus \pm the standard deviation calculated over the test dataset. The values in bold show the best results. The na entry means that LPD cannot be applied after changing the angular sampling pattern.

PSNR	FBP	LPD	Post	aNETT
Sparse view	$23.8\pm 1.3$	$37.9\pm 1.2$	$37.1\pm 0.9$	$37.1\pm 1.0$
Low dose	$36.9\pm 1.6$	$43.6\pm 1.3$	$44.1\pm 1.5$	$43.9\pm 1.3$
Universality	$32.4\pm 1.6$	na	$37.7\pm 0.8$	$38.3\pm 1.0$

Equations

$$y^{\delta }=\mathbf{K}x+{\eta }^{\delta },$$

$${\mathcal{N}}_{\alpha ,y^{\delta }}(x):=\mathcal{D}(\mathbf{K}x,y^{\delta })+\alpha \mathcal{Q}(\mathbf{E}(x)).$$

$${\mathcal{A}}_{\alpha ,y^{\delta }}(x):=\mathcal{D}(\mathbf{K}x,y^{\delta })+\alpha \left(\mathcal{Q}(\mathbf{E}(x))+\frac{c}{2}\Vert x-(\mathbf{D}\circ \mathbf{E})(x){\Vert }_2^2\right).$$

$$\mathcal{R}:\mathbb{X}\to [0,\infty ]:x\mapsto \mathcal{Q}(\mathbf{E}(x))+\frac{c}{2}\Vert x-(\mathbf{D}\circ \mathbf{E})(x){\Vert }_2^2,$$

$$\begin{gathered} {\mathcal{A}}_{\alpha ,y}(x^{‡})\le \underset{n\to \infty }{\mathrm{liminf}}\mathcal{D}(\mathbf{K}{x}_{\tau (n)},y_{\tau (n)})+\alpha \underset{n\to \infty }{\mathrm{liminf}}\mathcal{R}(x_{\tau (n)}) \\ \le \underset{n\to \infty }{\mathrm{limsup}}\mathcal{D}(\mathbf{K}{x}_{\tau (n)},y_{\tau (n)})+\alpha \mathcal{R}(x_{\tau (n)})\le \underset{n\to \infty }{\mathrm{limsup}}{\mathcal{A}}_{\alpha ,y_{\tau (n)}}(x)={\mathcal{A}}_{\alpha ,y}(x). \end{gathered}$$

$$\begin{gathered} \underset{n\to \infty }{\mathrm{limsup}}\alpha \mathcal{R}(x_{\tau (n)})\le \underset{n\to \infty }{\mathrm{limsup}}{\mathcal{A}}_{\alpha ,y_{\tau (n)}}(x_{\tau (n)})-\underset{n\to \infty }{\mathrm{liminf}}\mathcal{D}(\mathbf{K}{x}_{\tau (n)},y_{\tau (n)}) \\ \le {\mathcal{A}}_{\alpha ,y}(x^{‡})-\mathcal{D}(\mathbf{K}{x}^{‡},y)=\alpha \mathcal{R}(x^{‡}). \end{gathered}$$

$$\forall y_{0},y_{1},y_{2}\in\mathbb{Y}\colon\quad\mathcal{D}(y_{0},y_{1})\leq q\cdot\bigl{(}\mathcal{D}(y_{0},y_{2})+\mathcal{D}(y_{2},y_{1})\bigr{)}\,.$$

$$\mathcal{D}(y_0,y_1):=H(y_1)\Vert y_0-y_1{\Vert }^2,$$

$$x^{‡}\in \mathrm{argmin}\{\mathcal{R}(x):x\in \mathrm{dom}(\mathcal{R})\wedge \mathbf{K}x=y\}.$$

$$\begin{gathered} {\alpha }_n\mathcal{R}(x_n)\le \mathcal{D}(\mathbf{K}{x}_n,y_n)+{\alpha }_n\mathcal{R}(x_n)\le \mathcal{D}(\mathbf{K}{x}^{‡},y_n)+{\alpha }_n\mathcal{R}(x^{‡}) \\ =\mathcal{D}(y,y_n)+{\alpha }_n\mathcal{R}(x^{‡})\le {\delta }_n+{\alpha }_n\mathcal{R}(x^{‡}). \end{gathered}$$

$$\mathcal{R}(x^{‡})\le \mathcal{R}(x_{\star })\le \underset{n\to \infty }{\mathrm{liminf}}\mathcal{R}(x_{\tau (n)})\le \underset{n\to \infty }{\mathrm{liminf}}\frac{{\delta }_{\tau (n)}}{{\alpha }_{\tau (n)}}+\mathcal{R}(x^{‡})=\mathcal{R}(x^{‡})$$

$$\displaystyle\forall x\in\mathbb{X}\colon\quad\Delta_{\mathcal{F}}(x,x_{\star})\coloneqq\bigl{|}\mathcal{F}(x)-\mathcal{F}(x_{\star})-\mathcal{F}^{\prime}(x)(x-x_{\star})\bigr{|}\,.$$

$$\forall x\in \mathbb{X}:\Vert x-x^{‡}\Vert \le ϵ\Rightarrow {\Delta }_{\mathcal{R}}(x,x^{‡})\le \mathcal{R}(x)-\mathcal{R}(x^{‡})+c\sqrt{\mathcal{D}(\mathbf{K}x,\mathbf{K}{x}^{‡})}.$$

$$\alpha {\Delta }_{\mathcal{R}}(x_{\alpha }^{\delta },x^{‡})$$

$$\displaystyle=\mathcal{R}(x_{\alpha}^{\delta})-\mathcal{D}(\mathbf{K}x_{\alpha}^{\delta},y^{\delta})-\Bigl{(}\mathcal{R}(x_{\alpha}^{\delta})-\mathcal{D}(\mathbf{K}x_{\alpha}^{\delta},y^{\delta})\Bigr{)}+c\alpha\sqrt{\mathcal{D}(\mathbf{K}x_{\alpha}^{\delta},\mathbf{K}x^{\ddagger})}$$

$$\le {\delta }^2-\mathcal{D}(\mathbf{K}{x}_{\alpha }^{\delta },y^{\delta })+c\alpha \delta \sqrt{q}+c\alpha \sqrt{q\mathcal{D}(\mathbf{K}{x}_{\alpha }^{\delta },y^{\delta })}.$$

$$\forall x\in \mathbb{X}:\mathcal{R}(x^{‡})-\mathcal{R}(x)\le (L/a)\Vert \mathbf{K}x-\mathbf{K}{x}^{‡}\Vert$$

$$\forall x\in \mathbb{X}:|⟨{\mathcal{R}}^{\prime }(x^{‡}),x-x^{‡}⟩|\le (\Vert {\mathcal{R}}^{\prime }(x^{‡})\Vert /a)\Vert \mathbf{K}x-\mathbf{K}{x}^{‡}\Vert .$$

$${\Delta }_{\mathcal{R}}(x,x^{‡})$$

$$\le |\mathcal{R}(x)-\mathcal{R}(x^{‡})|+|⟨{\mathcal{R}}^{\prime }(x^{‡}),x^{‡}-x⟩|$$

$$\le \mathcal{R}(x)-\mathcal{R}(x^{‡})+(2L+\Vert {\mathcal{R}}^{\prime }(x^{‡})\Vert )a^{-1}\Vert \mathbf{K}{x}^{‡}-\mathbf{K}x\Vert ,$$

$${\theta }^{\ast }\in \underset{\theta }{\mathrm{argmin}}\frac{1}{m}\sum _{i=1}^m\Vert ({\mathbf{D}}_{\theta }\circ {\mathbf{E}}_{\theta })(x_i+a_i{ϵ}_i)-x_i{\Vert }_2^2+\nu (1-a_i)\mathcal{Q}({\mathbf{E}}_{\theta }(x_i))+\beta \Vert \theta {\Vert }_2^2,$$

$${\theta }^{\ast }\in \mathrm{argmin}\frac{1}{m}\sum _{i=1}^m\Vert {\mathbf{U}}_{\theta }(z_i)-x_i{\Vert }_2^2+\gamma \Vert \theta {\Vert }_2^2,$$

$$\{\begin{array}{rlrl}& {\mathrm{argmin}}_{x,\xi }& & \mathcal{D}(\mathbf{K}x,y^{\delta })+\alpha \mathcal{Q}(\xi )+\alpha c/2\Vert x-\mathbf{N}(x){\Vert }^2\\ & s.t.& & \mathbf{E}(x)=\xi .\end{array}$$

$$\mathrm{PSNR}(x,x_{\mathrm{rec}}):=20{\log }_{10}\left(\frac{\max x}{\Vert x-x_{\mathrm{rec}}{\Vert }_2}\right).$$

Full text

Keywords: inverse problems, regularization, stability guarantees, convergence rates, learned regularizer, computed tomography, neural networks

1 Introduction

Various applications in medical imaging, remote sensing and elsewhere require solving inverse problems of the form

[TABLE]

where $\mathbf{K}\colon\mathbb{X}\to\mathbb{Y}$ is an operator between Hilbert spaces modeling the forward problem, $\eta^{\delta}$ is the data perturbation, $y^{\delta}\in\mathbb{Y}$ is the noisy data and $x\in\mathbb{X}$ is the sought for signal. Inverse problems are well analyzed and several established approaches for its stable solution exist [13, 39]. Recently, neural networks and deep learning appeared as new paradigms for solving inverse problems [7, 17, 31, 29, 45]. Several approaches based on deep learning have been developed, including post-processing networks [26, 21, 6, 30, 37], regularizing null-space networks [41, 42], plug-and-play priors [44, 9, 38], deep image priors [12, 43], variational networks [23, 18], network cascades [40, 24], learned iterative schemes [2, 4, 10, 20, 25, 47] and learned regularizers [5, 27, 28, 33].

Classical deep learning approaches may lack data consistency for unknowns very different from the training data. To address this issue, in [27] a deep learning approach named NETT (NETwork Tikhonov) regularization has been introduced which considers minimizers of the NETT functional

[TABLE]

Here, $\mathcal{D}$ is a similarity measure, $\mathbf{E}\colon\mathbb{X}\rightarrow\Xi$ is a trained neural network, $\mathcal{Q}\colon\Xi\rightarrow[0,\infty]$ a functional and $\alpha>0$ a regularization parameter. In [27] it is shown that under suitable assumptions, NETT yields a convergent regularization method. This in particular includes provable stability guarantees and error estimates. Moreover, a training strategy has been proposed, where $\mathbf{E}$ is trained such that $\mathcal{Q}\circ\mathbf{E}$ favors artifact-free reconstructions over reconstructions with artifacts.

1.1 The augmented NETT

One of the main assumptions for the analysis of [27] is the coercivity of the regularizer $\mathcal{Q}\circ\mathbf{E}$ which requires special care in network design and training. In order to overcome this limitation, we propose an augmented form of the regularizer for which we are able to rigorously prove coercivity. More precisely, for fixed $c>0$, we consider minimizers $x_{\alpha}^{\delta}$ of the augmented NETT functional

[TABLE]

Here, $\mathcal{D}\colon\mathbb{Y}\times\mathbb{Y}\rightarrow[0,\infty]$ is a similarity measure and $\mathbf{D}\circ\mathbf{E}\colon\mathbb{X}\rightarrow\mathbb{X}$ is an encoder-decoder network trained such that for any signal $x$ on a signal manifold we have $(\mathbf{D}\circ\mathbf{E})(x)\simeq x$ and that $\mathcal{Q}(\mathbf{E}(x))$ is small. We term this approach augmented NETT (aNETT) regularization. In this work we provide a mathematical convergence analysis for aNETT, present a novel modular training strategy and investigate its practical performance.

The term $\mathcal{Q}(\mathbf{E}(x))$ implements learned prior knowledge on the encoder coefficients, while smallness of $\lVert x-(\mathbf{D}\circ\mathbf{E})(x)\rVert_{2}^{2}$ forces $x$ to be close to the signal manifold. The latter term also guarantees the coercivity of (1.3). In the original NETT version (1.2), coercivity of the regularizer requires coercivity conditions on the network involved. Indeed, in the numerical experiments, the authors of [27] observed a semi-convergence behaviour when minimizing (1.2), so early stopping of the iterative minimization scheme has been used as additional regularization. We attribute this semi-convergence behavior to the non-coercivity of the regularization term. In the present, paper we address this issue systematically by augmentation of the NETT functional which guarantees coercivity and allows a more stable minimization. Coercivity is also one main ingredient for the mathematical convergence analysis.

An interesting practical instance of aNETT takes $\mathcal{Q}$ as a weighted $\ell^{q}$-norm enforcing sparsity of the encoding coefficients [11, 16]. An important example for the similarity measure is given by the squared norm distance, which from a statistical viewpoint can be motivated by a Gaussian white noise model. General similarity measures allow us to adapt to different noise models which can be more appropriate for certain problems.

1.2 Main contributions

The contributions of this paper are threefold. As described in more detail below, we introduce the aNETT framework, mathematically analyze its convergence, and propose a practical implementation that is applied to tomographic limited data problems.

• •

The first contribution is to introduce the structure of the aNETT regularizer $\mathcal{R}(x)=\mathcal{Q}(\mathbf{E}(x))+(c/2)\,\lVert x-(\mathbf{D}\circ\mathbf{E})(x)\rVert_{2}^{2}$. The combination of the two terms in the regularizer, even in the case of a linear encoder seems to be new. The term $\mathcal{Q}(\mathbf{E}(x))$ enforces regularity of the analysis coefficients, which is an ingredient in most of existing variational regularization techniques. For example, this includes sparse regularization in frames or dictionaries, regularization with Sobolev norms or total variation regularization. On the other hand, the augmented term $\lVert x-(\mathbf{D}\circ\mathbf{E})(x)\rVert_{2}^{2}$ penalized distance to the signal manifold. It is the combination of these two terms that results in a stable reconstruction scheme without the need of strong assumptions on the involved networks.

• •

The second main contribution is the theoretical analysis of aNETT (1.3) in the context of regularization theory. We investigate the case where the image domain of the encoder is given by $\Xi=\ell^{2}(\Lambda)$ for some countable set $\Lambda$, and $\mathcal{Q}$ is a coercive functional measuring the complexity of the encoder coefficients. The presented analysis is in the spirit of the analysis of NETT given in [27]. However, opposed to NETT, the required coercivity property is derived naturally for the class of considered regularizers. This supports the use of the regularizer $\mathcal{R}$ also from a theoretical side. Moreover, the convergence rates results presented here uses assumptions significantly different from [27]. While we present our analysis for the transform domain $\Xi=\ell^{2}(\Lambda)$ we could replace the encoder space by a general Hilbert or Banach space.

• •

As a third main contribution we propose a modular strategy for training $\mathbf{D}\circ\mathbf{E}$ together with a possible network architecture. First, independent of the given inverse problem, we train a $\mathcal{Q}$-penalized autoencoder that learns representing signals from the training data with low complexity. In the second step, we train a task-specific network which can be adapted to specific inverse problem at hand. In our numerical experiments, we empirically found this modular training strategy to be superior to directly adapting the autoencoder to the inverse problem. For the $\mathcal{Q}$-penalized autoencoder, we train the modified version described in [35] of the tight frame U-Net of [19] in a way such that $\mathcal{Q}$ poses additional constraints on the autoencoder during the training process.

1.3 Outline

In Section 2 we present the mathematical convergence analysis of aNETT. In particular, as an auxiliary result, we establish the coercivity of the regularization term. Moreover, we prove stability and derive convergence rates. Section 3 presents practical aspects for aNETT. We propose a possible architecture and training strategy for the networks, and a possible ADMM based scheme to obtain minimizers of the aNETT functional. In Section 4, we present reconstruction results and compare aNETT with other deep learning based reconstruction methods. The paper concludes with a short summary and discussion. Parts of this paper were presented at the ISBI 2020 conference and the corresponding proceedings [34]. Opposed to the proceedings, this article treats a general similarity measure $\mathcal{D}$ and considers a general complexity measure $\mathcal{Q}$. Further, all proofs and all numerical results presented in this paper are new.

2 Mathematical analysis

In this section we prove the stability and convergence of aNETT as regularization method. Moreover, we derive convergence rates in the form of quantitative error estimates between exact solutions for noise-free data and aNETT regularized solutions for noisy data.

2.1 Assumptions and coercivity results

For our convergence analysis we make use of the following assumptions on the underlying spaces and operators involved.

Condition 2.1** (General assumptions for aNETT).**

1. (A1)

$\mathbb{X}$ and $\mathbb{Y}$ are Hilbert spaces. 2. (A2)

$\Xi=\ell^{2}(\Lambda)$ for a countable $\Lambda$. 3. (A3)

$\mathbf{K}\colon\mathbb{X}\rightarrow\mathbb{Y}$ is weakly sequentially continuous. 4. (A4)

$\mathbf{E}\colon\mathbb{X}\rightarrow\Xi$ is weakly sequentially continuous. 5. (A5)

$\mathbf{D}\colon\Xi\rightarrow\mathbb{X}$ is weakly sequentially continuous. 6. (A6)

$\mathcal{Q}\colon\Xi\rightarrow[0,\infty]$ is coercive and weakly sequentially lower semi-continuous.

We set $\mathbf{N}\coloneqq\mathbf{D}\circ\mathbf{E}$ and, for given $c>0$, define

[TABLE]

which we refer to as the aNETT (or augmented NETT) regularizer.

According to (A4)-(A6), the aNETT regularizer is weakly sequentially lower semi-continuous. As a main ingredient for our analysis we next prove its coercivity.

Theorem 2.2** (Coercivity of the aNETT regularizer $\mathcal{R}$).**

If Condition 2.1 holds, then $\mathcal{R}\colon\mathbb{X}\rightarrow[0,\infty]$ is coercive.

Proof.

Let $(x_{n})_{n\in\mathds{N}}$ be some sequence in $\mathbb{X}$ such that $(\mathcal{R}(x_{n}))_{n\in\mathds{N}}$ is bounded. Then by definition of $\mathcal{R}$ it follows that $(\mathcal{Q}(\mathbf{E}(x_{n})))_{n\in\mathds{N}}$ is bounded and by coercivity of $\mathcal{Q}$ we have that $(\mathbf{E}(x_{n}))_{n\in\mathds{N}}$ is also bounded. By assumption, $\mathbf{D}$ is weakly sequentially continuous and thus $(\lVert\mathbf{N}(x_{n})\rVert)_{n\in\mathds{N}}$ must be bounded. Using that $c>0$, we obtain the inequality $\lVert x_{n}\rVert^{2}\leq 2\lVert x_{n}-\mathbf{N}(x_{n})\rVert^{2}+2\lVert\mathbf{N}(x_{n})\rVert^{2}\leq(4/c)\,\mathcal{R}(x_{n})+2\lVert\mathbf{N}(x_{n})\rVert^{2}$. This shows that $(x_{n})_{n\in\mathds{N}}$ is bounded and therefore that $\mathcal{R}$ is coercive. ∎

Example 2.3** (Sparse aNETT regularizer).**

To obtain a sparsity promoting regularizer we can choose $\mathcal{Q}(\xi)=\lVert\xi\rVert_{1;w}\coloneqq\sum_{\lambda\in\Lambda}w_{\lambda}\lvert\xi_{\lambda}\rvert^{q}$ where $q\in[1,2]$ and $\inf_{\lambda}w_{\lambda}>0$. Since $q\in[1,2]$ we have $\lVert\,\cdot\,\rVert^{2}\leq(\inf_{\lambda}w_{\lambda})^{-1/q}\,\lVert\,\cdot\,\rVert_{1;w}$ and hence $\lVert\,\cdot\,\rVert_{1;w}$ is coercive. As a sum of weakly sequentially lower semi-continuous functionals it is also weakly sequentially lower semi-continuous [16]. Therefore Condition (A6) is satisfied for the weighted $\ell^{q}$-norm. Together with Theorem 2.2, we conclude that the resulting weighted sparse aNETT regularizer $x\mapsto\lVert\mathbf{E}(x)\rVert_{1;w}+(c/2)\,\lVert x-(\mathbf{D}\circ\mathbf{E})(x)\rVert_{2}^{2}$ is a coercive and weakly sequentially lower semi-continuous functional.

For the further analysis we will make the following assumptions including the similarity measure $\mathcal{D}\colon\mathbb{Y}\times\mathbb{Y}\rightarrow[0,\infty]$.

Condition 2.4** (Similarity measure).**

1. (B1)

$\forall y_{0},y_{1}\colon\mathcal{D}(y_{0},y_{1})=0\Leftrightarrow y_{0}=y_{1}$. 2. (B2)

$\mathcal{D}$ is sequentially lower semi-continuous with respect to the weak topology in the first and the norm topology in the second argument. 3. (B3)

$\forall(y_{n})_{n\in\mathds{N}}\in\mathbb{Y}^{\mathds{N}}\colon(\mathcal{D}(y,y_{n})\rightarrow 0\Rightarrow y_{n}\rightarrow y$ as $n\to\infty$). 4. (B4)

$\mathcal{D}(y,y_{n})\rightarrow 0$ as $n\to\infty$ $\Rightarrow$ ($\forall z\in\mathbb{Y}\colon\mathcal{D}(z,y)<\infty\Rightarrow\mathcal{D}(z,y_{n})\rightarrow\mathcal{D}(z,y)$). 5. (B5)

$\forall y\in\mathbb{Y}\;\forall\alpha>0\colon$ $\exists x\in\mathbb{X}$ with $\mathcal{D}(\mathbf{K}x,y)+\alpha\mathcal{R}(x)<\infty$.

While (B1)-(B4) restrict the choice of the similarity measure, (B5) is a technical assumption involving the forward operator, the regularizer and the similarity measure, that is required for the existence of minimizers. For a more detailed discussion of these assumptions we refer to [36].

Example 2.5** (Similarity measures using the norm).**

The classical example of a similarity measure satisfying (B1)-(B4) is given by $\mathcal{D}(y_{0},y_{1})=\lVert y_{0}-y_{1}\rVert^{p}$ for some $p\geq 1$ and more generally by $\mathcal{D}(y_{0},y_{1})=\psi(\lVert y_{0}-y_{1}\rVert)$, where $\psi\colon[0,\infty)\rightarrow[0,\infty)$ is a continuous and monotonically increasing function that satisfies $\forall t\geq 0\colon\psi(t)=0\Leftrightarrow t=0$.

Taking into account Theorem 2.2, Conditions 2.1 and 2.4 imply that the aNETT functional $\mathcal{A}_{\alpha,y}$ defined by (1.3), (2.1) is proper, coercive and weakly sequentially lower semi-continuous. This in particular implies the existence of minimizers of $\mathcal{A}_{\alpha,y}$ for all data $y\in\mathbb{Y}$ and regularization parameters $\alpha>0$.

2.2 Stability

Next we prove the stability of minimizing the aNETT functional $\mathcal{A}_{\alpha,y}$ regarding perturbations of the data $y$.

Theorem 2.6** (Stability).**

Let Conditions 2.1 and 2.4 hold, $y\in\mathbb{Y}$ and $\alpha>0$. Moreover, let $(y_{n})_{n\in\mathds{N}}\in\mathbb{Y}^{\mathds{N}}$ be a sequence of perturbed data with $\mathcal{D}(y,y_{n})\rightarrow 0$ and consider minimizers $x_{n}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha,y_{n}}$. Then the sequence $(x_{n})_{n\in\mathds{N}}\in\mathbb{X}^{\mathds{N}}$ has at least one weak accumulation point and weak accumulation points are minimizers of $\mathcal{A}_{\alpha,y}$. Moreover, for any weak accumulation point $x^{\ddagger}$ of $(x_{n})_{n\in\mathds{N}}$ and any subsequence $(x_{\tau(n)})_{n\in\mathds{N}}$ with $x_{\tau(n)}\rightharpoonup x^{\ddagger}$ we have $\mathcal{R}(x_{\tau(n)})\rightarrow\mathcal{R}(x^{\ddagger})$.

Proof.

Let $x_{\star}\in\mathbb{X}$ be such that $\mathcal{A}_{\alpha,y}(x_{\star})<\infty$. By definition of $x_{n}$ we have $\mathcal{A}_{\alpha,y_{n}}(x_{n})\leq\mathcal{A}_{\alpha,y_{n}}(x_{\star})$. Since by assumption $\mathcal{D}(y,y_{n})\rightarrow 0$ and $\mathcal{D}(\mathbf{K}x_{\star},y)<\infty$, we have $\mathcal{D}(\mathbf{K}x_{\star},y_{n})\rightarrow\mathcal{D}(\mathbf{K}x_{\star},y)$. This implies that $\mathcal{A}_{\alpha,y_{n}}(x_{\star})$ is bounded by some positive constant $m_{\star}$ for sufficiently large $n$. By definition of $\mathcal{A}_{\alpha,y_{n}}$ we have $\alpha\mathcal{R}(x_{n})\leq\mathcal{A}_{\alpha,y_{n}}(x_{n})\leq m_{\star}+\alpha\mathcal{R}(x_{\star})$. Since $\mathcal{R}$ is coercive it follows that $(x_{n})_{n\in\mathds{N}}$ is a bounded sequence and hence it has a weakly convergent subsequence.

Let $(x_{\tau(n)})_{n\in\mathds{N}}$ be a weakly convergent subsequence of $(x_{n})_{n\in\mathds{N}}$ and denote its limit by $x^{\ddagger}$. By the lower semi-continuity we get $\mathcal{D}(\mathbf{K}x^{\ddagger},y)\leq\liminf_{n\to\infty}\mathcal{D}(\mathbf{K}x_{\tau(n)},y_{\tau(n)})$ and $\mathcal{R}(x^{\ddagger})\leq\liminf_{n\to\infty}\mathcal{R}(x_{\tau(n)})$. Thus for all $x\in\mathbb{X}$ with $\mathcal{A}_{\alpha,y}(x)<\infty$ we have

[TABLE]

This shows that $x^{\ddagger}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha,y}$ and, by considering $x=x^{\ddagger}$ in the above displayed equation, that $\mathcal{A}_{\alpha,y_{\tau(n)}}(x_{\tau(n)})\to\mathcal{A}_{\alpha,y}(x^{\ddagger})$. Moreover, we have

[TABLE]

This shows $\mathcal{R}(x_{\tau(n)})\to\mathcal{R}(x^{\ddagger})$ as $n\to\infty$ and concludes the proof. ∎

In the following we say that the similarity measure $\mathcal{D}$ satisfies the quasi triangle-inequality if there is some $q\geq 1$ such that

[TABLE]

While this property is essential for deriving convergence rate results, we will show below that it is not enough to guarantee stability of minimizing the augmented NETT functional in the sense of Theorem 2.6. Note that [27] assumes the quasi triangle-inequality (2.2) instead of Condition (B4). The following remarks shows that (2.2) is not sufficient for the stability result of Theorem 2.6 to hold and therefore Condition (B4) has to be added to the list of assumptions in [27] required for the stability.

Example 2.7** (Instability in the absence of Condition (B4)).**

Consider the similarity measure $\mathcal{D}\colon\mathbb{Y}\times\mathbb{Y}\to[0,\infty]$ defined by

[TABLE]

where $H\colon\mathbb{Y}\to[0,1]$ is defined by $H(y_{1})=1$ if $\lVert y_{1}\rVert\leq 1$ and $H(y_{1})=2$ otherwise. Moreover, choose $\mathbb{X}=\mathbb{Y}$, let $\mathbf{K}=\mathrm{Id}$ be the identity operator and suppose the regularizer takes the form $\mathcal{R}=\lVert\,\cdot\,\rVert^{2}$.

• •

The similarity measure defined in (2.3) satisfies (B1)-(B3): Convergence with respect to $\mathcal{D}$ is equivalent to convergence in norm which implies that (B3) is satisfied. Moreover, we have $\forall y_{0},y_{1}\colon\mathcal{D}(y_{0},y_{1})=0\Leftrightarrow y_{0}=y_{1}$, which is (B1). Consider sequences $z_{n}\to z$ and $y_{n}\rightharpoonup y$. The sequential lower semi-continuity stated in (B2) can be derived by separately looking at the cases $\lVert z\rVert\leq 1$ and $\lVert z\rVert>1$. In the first case, by the continuity of the norm we have $\forall n\in\mathds{N}\colon H(z)\leq H(z_{n})$. In the second case, we have $\lVert y_{n}\rVert>1$ for $n$ sufficiently large. In both cases, the lower semi-continuity property follows from the weak lower semi-continuity property of the norm.

• •

The similarity measure defined by (2.3) does not satisfy (B4): To see this, we define $y_{n}\coloneqq(1+\delta_{n})y$ where $\lVert y\rVert=1$ and $\delta_{n}>0$ is taken as a non-increasing sequence converging to zero. We have $y_{n}\to y$ and hence also $\mathcal{D}(y,y_{n})\to 0$ as $n\to\infty$. For any $z\in\mathbb{Y}$ we have $\mathcal{D}(z,y)<\infty$ and $\mathcal{D}(z,y_{n})=2\lVert z-y_{n}\rVert^{2}\rightarrow 2\lVert z-y\rVert^{2}=2\mathcal{D}(z,y)$ as $n\to\infty$. In particular, $\mathcal{D}(z,y_{n})$ does not converge to $\mathcal{D}(z,y)$ if $z\neq y$ and therefore (B5) does not hold. In summary, all requirement for Theorem 2.6 are satisfied, except of the continuity assumption (B4).

• •

We have $\lVert y_{0}-y_{1}\rVert^{2}\leq\mathcal{D}(y_{0},y_{1})\leq 2\lVert y_{0}-y_{1}\rVert^{2}$ which implies that the similarity measure satisfies the quasi triangle-inequality (2.2). However as shown next, this is not sufficient for stable reconstruction in the sense of Theorem 2.6. To that end, let $y_{n}\coloneqq(1+\delta_{n})y$ with $\lVert y\rVert=1$ and $\delta_{n}\downarrow 0$ be as above and let $\alpha>0$. In particular, $\mathcal{D}(x,y_{n})=2\lVert x-y_{n}\rVert^{2}$ and $\mathcal{D}(x,y)=\lVert x-y_{n}\rVert^{2}$. Therefore the minimizer of $\operatorname*{arg\,min}\mathcal{D}(x,y_{n})+\alpha\lVert x\rVert^{2}$ with perturbed data $y_{n}$ is given by $x_{n}=y_{n}/(1+\alpha/2)$ and the minimizer of $\operatorname*{arg\,min}\mathcal{D}(x,y_{n})+\alpha\lVert x\rVert^{2}$ for data $y$ is given by $x^{\ddagger}=y/(1+\alpha)$. We see that $x_{n}\to y/(1+\alpha/2)$, which is clearly different from $x^{\ddagger}$. In particular, minimizing $\mathcal{D}(\,\cdot\,,y)+\alpha\lVert\,\cdot\,\rVert^{2}$ does not stably depend on data $y$. Theorem 2.6 states that stability holds if (B4) is satisfied.

While the above example may seem somehow constructed, it shows that one has to be careful when choosing the similarity measure in order to obtain a stable reconstruction scheme.

2.3 Convergence

In this subsection we prove that minimizers of the aNETT functional for noisy data convergence to $\mathcal{R}$-minimizing solution of the equation $\mathbf{K}x=y$ as the noise level goes to zero and the regularization parameter is chosen properly. Here and below, we use the following notation.

Definiton 2.8** ($\mathcal{R}$-minimizing solutions).**

For $y\in\mathbb{Y}$, we call an element $x^{\ddagger}\in\operatorname{dom}(\mathcal{R})$ an $\mathcal{R}$-minimizing solution of the equation $\mathbf{K}x=y$ if

[TABLE]

An $\mathcal{R}$-minimizing solution always exists provided that data satisfies $y\in\mathbf{K}(\operatorname{dom}(\mathcal{R}))$, which means that the equation $\mathbf{K}x=y$ has at least on solution with finite value of $\mathcal{R}$. To see this, consider a sequence of solutions $(x_{n})_{n\in\mathds{N}}$ with $\mathcal{R}(x_{n})\rightarrow\inf\{\mathcal{R}(x)\colon\mathbf{K}x=y\wedge\mathcal{R}(x)<\infty\}$. Since $\mathcal{R}$ is coercive there exists a weakly convergent subsequence $(x_{\tau(n)})_{n\in\mathds{N}}$ with weak limit $x^{\ddagger}$. Using the weak sequential lower semi-continuity of $\mathcal{R}$ one concludes that $x^{\ddagger}$ is an $\mathcal{R}$-minimizing solution.

We first show weak convergence.

Theorem 2.9** (Weak convergence of aNETT).**

Suppose Conditions 2.1 and 2.4 are satisfied. Let $y\in\mathbf{K}(\operatorname{dom}(\mathcal{R}))$, $(\delta_{n})_{n\in\mathds{N}}\in(0,\infty)^{\mathds{N}}$ with $\delta_{n}\rightarrow 0$ and let $(y_{n})_{n\in\mathds{N}}\in\mathbb{Y}^{\mathds{N}}$ satisfy $\mathcal{D}(y,y_{n})\leq\delta_{n}$. Choose $\alpha_{n}>0$ such that $\lim_{n\to\infty}\delta_{n}/\alpha_{n}=\lim_{n\to\infty}\alpha_{n}=0$ and let $x_{n}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha_{n},y_{n}}$. Then the following hold:

1. (a)

$(x_{n})_{n\in\mathds{N}}$* has at least one weakly convergent subsequence.* 2. (b)

All accumulation points of $(x_{n})_{n\in\mathds{N}}$ are $\mathcal{R}$-minimizing solutions of $\mathbf{K}x=y$. 3. (c)

For every convergent subsequence $(x_{\tau(n)})_{n\in\mathds{N}}$ it holds $\mathcal{R}(x_{\tau(n)})\rightarrow\mathcal{R}(x^{\ddagger})$. 4. (d)

If the $\mathcal{R}$-minimizing solution $x^{\ddagger}$ is unique then $x_{n}\rightharpoonup x^{\ddagger}$.

Proof.

(a): Because $y\in\mathbf{K}(\operatorname{dom}(\mathcal{R}))$, there exists an $\mathcal{R}$-minimizing solution of the equation $\mathbf{K}x=y$ which we denote by $x^{\ddagger}$. Because $x_{n}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha_{n},y_{n}}$ we have

[TABLE]

Because $\alpha_{n},\delta_{n}\to 0$ this shows that $(\mathcal{R}(x_{n}))_{n\in\mathds{N}}$ is bounded. Due to the coercivity of the aNETT regularizer (see Theorem 2.2), this implies that $(x_{n})_{n\in\mathds{N}}$ has a weakly convergent subsequence.

(b), (c): Let $(x_{\tau(n)})_{n\in\mathds{N}}$ be a weakly convergent subsequence of $(x_{n})_{n\in\mathds{N}}$ with limit $x_{\star}$. From the weak lower semi-continuity we get $\mathcal{D}(\mathbf{K}x_{\star},y)\leq\liminf_{n\to\infty}\mathcal{D}(\mathbf{K}x_{\tau(n)},y_{\tau(n)})+\alpha_{\tau(n)}\mathcal{R}(x_{\tau(n)})\leq\liminf_{n\to\infty}\mathcal{D}(\mathbf{K}x^{\ddagger},y_{\tau(n)})+\alpha_{\tau(n)}\mathcal{R}(x^{\ddagger})=0$, which shows that $x_{\star}$ is a solution of $\mathbf{K}x=y$. Moreover,

[TABLE]

where for the second last inequality we used (2.4) and for the last equality we used that $\delta_{n}/\alpha_{n}\to 0$. Therefore, $x_{\star}$ is an $\mathcal{R}$-minimizing solution of the equation $\mathbf{K}x=y$. In a similar manner we derive $\mathcal{R}(x_{\star})\leq\liminf_{n\to\infty}\mathcal{R}(x_{\tau(n)})\leq\limsup_{n\to\infty}\delta_{\tau(n)}/\alpha_{\tau(n)}+\mathcal{R}(x^{\ddagger})=\mathcal{R}(x_{\star})$ which shows $\mathcal{R}(x_{\tau(n)})\rightarrow\mathcal{R}(x_{\star})$.

(d): If $\mathbf{K}x=y$ has a unique $\mathcal{R}$-minimizing solution $x^{\ddagger}$, then every subsequence of $(x_{n})_{n\in\mathds{N}}$ has itself a subsequence weakly converging to $x^{\ddagger}$, which implies that $(x_{n})_{n\in\mathds{N}}$ weakly converges to the $\mathcal{R}$-minimizing solution. ∎

Next we derive strong convergence of the regularized solutions. To this end we recall the absolute Bregman distance, the modulus of total nonlinearity and the total nonlinearity, defined in [27].

Definiton 2.10** (Absolute Bregman distance).**

Let $\mathcal{F}\colon\mathbb{X}\rightarrow[0,\infty]$ be Gâteaux differentiable at $x_{\star}\in\mathbb{X}$. The absolute Bregman distance $\Delta_{\mathcal{F}}(\,\cdot\,,x_{\star})\colon\mathbb{X}\rightarrow[0,\infty]$ at $x_{\star}$ with respect to $\mathcal{F}$ is defined by

[TABLE]

Here and below $\mathcal{F}^{\prime}(x_{\star})$ denotes the Gâteaux derivative of $\mathcal{F}$ at $x_{\star}$.

Definiton 2.11** (Modulus of total nonlinearity and total nonlinearity).**

Let $\mathcal{F}\colon\mathbb{X}\rightarrow[0,\infty]$ be Gâteaux differentiable at $x_{\star}\in\mathbb{X}$. We define the modulus of total nonlinearity of $\mathcal{F}$ at $x_{\star}$ as $\nu_{\mathcal{F}}(x_{\star},\,\cdot\,)\colon[0,\infty)\rightarrow[0,\infty)\colon t\mapsto\inf\{\Delta_{\mathcal{F}}(x,x_{\star})\colon\lVert x-x_{\star}\rVert=t\}$. We call $\mathcal{F}$ totally nonlinear at $x_{\star}$ if $\nu_{\mathcal{F}}(x_{\star},t)>0$ for all $t\in(0,\infty)$.

Using these definitions we get the following convergence result in the norm topology.

Theorem 2.12** (Strong convergence of aNETT).**

Let Conditions 2.1 and 2.4 hold, $y\in\mathbf{K}(\operatorname{dom}(\mathcal{R}))$ and let $\mathcal{R}$ be totally nonlinear at all $\mathcal{R}$-minimizing solutions of $\mathbf{K}x=y$. Let $(y_{n})_{n\in\mathds{N}},(x_{n})_{n\in\mathds{N}},(\alpha_{n})_{n\in\mathds{N}}$ be as in Theorem 2.9. Then there is a subsequence $(x_{\tau(n)})_{n\in\mathds{N}}$ which converges in norm to an $\mathcal{R}$-minimizing solution $x^{\ddagger}$ of $\mathbf{K}x=y$. If the $\mathcal{R}$-minimizing solution is unique, then $x_{n}\to x^{\ddagger}$ as $n\to\infty$.

Proof.

In [27, Proposition 2.9] it is shown that the total nonlinearity of $\mathcal{R}$ implies that for every bounded sequence $(z_{n})_{n\in\mathds{N}}$ with $\Delta_{\mathcal{R}}(z_{n},z)\rightarrow 0$ it holds that $z_{n}\to z$. Theorem 2.9 gives us a weakly converging subsequence $(x_{\tau(n)})_{n\in\mathds{N}}$ of $(x_{n})_{n\in\mathds{N}}$ with weak limit $x^{\ddagger}$ and $\mathcal{R}(x_{\tau(n)})\to\mathcal{R}(x^{\ddagger})$. By the definition of the absolute Bregman distance it follows that $\Delta_{\mathcal{R}}(x_{\tau(n)},x^{\ddagger})\rightarrow 0$ and hence, together with [27, Proposition 2.9], that $x_{\tau(n)}\rightarrow x^{\ddagger}$. If the $\mathcal{R}$-minimizing solution of $\mathbf{K}x=y$ is unique, then every subsequence has a subsequence converging to $x^{\ddagger}$ and hence the claim follows. ∎

2.4 Convergence rates

We will now prove convergence rates by deriving quantitative estimates for the absolute Bregman distance between $\mathcal{R}$-minimizing solutions for exact data and regularized solutions for noisy data. The convergence rates will be derived under the additional assumption that $\mathcal{D}$ satisfies the quasi triangle-inequality (2.2).

Proposition 2.13** (Convergence rates for aNETT).**

Let the assumptions of Theorem 2.12 be satisfied and suppose that $\mathcal{D}$ satisfies the quasi triangle-inequality (2.2) for some $q\geq 1$. Let $x^{\ddagger}\in\mathbb{X}$ be an $\mathcal{R}$-minimizing solution of $\mathbf{K}x=y$ such that $\mathcal{R}$ is Gâteaux differentiable at $x^{\ddagger}$ and assume there exist $\epsilon,c>0$ with

[TABLE]

For any $\delta>0$, let $y^{\delta}\in\mathbb{Y}$ be noisy data satisfying $\mathcal{D}(y^{\delta},\mathbf{K}x^{\ddagger})^{1/2}\leq\delta$, $\mathcal{D}(\mathbf{K}x^{\ddagger},y^{\delta})^{1/2}\leq\delta$, and write $x_{\alpha}^{\delta}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha,y^{\delta}}$. Then the following hold:

1. (a)

For sufficiently small $\alpha$, it holds $\Delta_{\mathcal{R}}(x_{\alpha}^{\delta},x^{\ddagger})\leq\delta^{2}/\alpha-c\delta\sqrt{q}+c^{2}\alpha q/4$. 2. (b)

If $\alpha\asymp\delta$, then \Delta_{\mathcal{R}}(x_{\alpha}^{\delta},x^{\ddagger})=\mathcal{O}\bigl{(}\delta\bigr{)} as $\delta\to 0$.

Proof.

By definition of $x_{\alpha}^{\delta}$ we have $\mathcal{D}(\mathbf{K}x_{\alpha}^{\delta},y^{\delta})+\alpha\mathcal{R}(x_{\alpha}^{\delta})-\alpha\mathcal{R}(x^{\ddagger})\leq\mathcal{D}(\mathbf{K}x^{\ddagger},y^{\delta})\leq\delta^{2}$. By Theorem 2.12 for sufficiently small $\alpha$ we can assume that $\lVert x_{\alpha}^{\delta}-x^{\ddagger}\rVert\leq\epsilon$ and hence

[TABLE]

Together with the inequality of arithmetic and geometric means $(a+b)/2\geq\sqrt{ab}$ for $a=\mathcal{D}(\mathbf{K}x_{\alpha}^{\delta},y^{\delta})$ and $b=c^{2}\alpha^{2}q/4$ this implies $\alpha\Delta_{\mathcal{R}}(x_{\alpha}^{\delta},x^{\ddagger})\leq\delta^{2}-c\alpha\delta\sqrt{q}+c^{2}\alpha^{2}q/4$ which shows (a). Item (b) is an immediate consequence of (a). ∎

The following results is our main convergence rates result. It is similar to Proposition [27, Theorem 3.1], but uses different assumptions.

Theorem 2.14** (Convergence rates for finite rank operators).**

Let the assumptions of Theorem 2.12 be satisfied, take $\mathcal{D}(y_{1},y_{2})=\lVert y_{1}-y_{2}\rVert^{2}$, assume $\mathbf{K}$ has finite dimensional range and that $\mathcal{R}$ is Lipschitz continuous and Gâteaux differentiable. For any $\delta>0$, let $y^{\delta}\in\mathbb{Y}$ be noisy data satisfying $\lVert y^{\delta}-\mathbf{K}x^{\ddagger}\rVert\leq\delta$ and write $x_{\alpha}^{\delta}\in\operatorname*{arg\,min}\mathcal{A}_{\alpha,y^{\delta}}$. Then for the parameter choice $\alpha\asymp\delta$ we have the convergence rates result \Delta_{\mathcal{R}}(x_{\alpha}^{\delta},x^{\ddagger})=\mathcal{O}\bigl{(}\delta\bigr{)} as $\delta\to 0$.

Proof.

According to Proposition 2.13, it is sufficient to show that (2.5) holds with $\lVert\mathbf{K}x-\mathbf{K}x^{\ddagger}\rVert$ in place of $\mathcal{D}(\mathbf{K}x,\mathbf{K}x^{\ddagger})^{1/2}$. For that purpose, let $\mathbf{P}$ denote the orthogonal projection onto the null-space $\ker(\mathbf{K})$ and let $L$ be a Lipschitz constant of $\mathcal{R}$. Since $\mathbf{K}$ restricted to $\ker(\mathbf{K})^{\bot}$ is injective with finite dimensional range, we can choose a constant $a>0$ such that $\forall z\in\ker(\mathbf{K})^{\bot}\colon\lVert\mathbf{K}z\rVert\geq a\lVert z\rVert$.

We first show the estimates

[TABLE]

To that end, let $x\in\mathbb{X}$ and write $x_{0}\coloneqq(x^{\ddagger}-\mathbf{P}x^{\ddagger})+\mathbf{P}x$. Then $\mathbf{K}x_{0}=\mathbf{K}x^{\ddagger}$. Since $x^{\ddagger}$ is an $\mathcal{R}$-minimizing solution, we have $\mathcal{R}(x^{\ddagger})-\mathcal{R}(x)\leq\mathcal{R}(x_{0})-\mathcal{R}(x)\leq L\lVert x_{0}-x\rVert$. Since $x_{0}-x\in\ker(\mathbf{K})^{\bot}$, we have $\lVert\mathbf{K}x^{\ddagger}-\mathbf{K}x\rVert=\lVert\mathbf{K}(x_{0}-x)\rVert\geq a\lVert x_{0}-x\rVert$. The last two estimates prove (2.6). Because $x^{\ddagger}$ is an $\mathcal{R}$-minimizing solution, we have $\langle\mathcal{R}^{\prime}(x^{\ddagger}),x^{\ddagger}-x\rangle=0$ whenever $x^{\ddagger}-x\in\ker(\mathbf{K})$. On the other hand, using that $\mathcal{R}$ is Gâteaux differentiable and that $\mathbf{K}$ has finite rank, shows $|\langle\mathcal{R}^{\prime}(x^{\ddagger}),x^{\ddagger}-x\rangle\rvert\leq\lVert\mathcal{R}^{\prime}(x^{\ddagger})\rVert\,a^{-1}\lVert\mathbf{K}(x^{\ddagger}-x)\rVert$ for $x^{\ddagger}-x\in\ker(\mathbf{K})^{\bot}$. This proves (2.7).

Inequality (2.6) implies $\lvert\mathcal{R}(x)-\mathcal{R}(x^{\ddagger})\rvert\leq\mathcal{R}(x)-\mathcal{R}(x^{\ddagger})+2\,(L/a)\,\lVert\mathbf{K}x-\mathbf{K}x^{\ddagger}\rVert$. Together with (2.7) this yields

[TABLE]

which proves (2.5) with $c=(2L+\lVert\mathcal{R}^{\prime}(x^{\ddagger})\rVert)/a$. ∎

Note that the theoretical results stated remain valid, if we replace $\mathcal{R}$ by a general coercive and weakly lower semi-continuous regularizer $\mathcal{R}\colon\mathbb{X}\rightarrow[0,\infty]$.

3 Practical realization

In this section we investigate practical aspects of aNETT. We present a possible network architecture together with a possible training strategy in the discrete setting. Further we discuss minimization of aNETT using the ADMM algorithm.

For the sake of clarity we restrict our discussion to the finite dimensional case where $\mathbb{X}=\mathds{R}^{N\times N}$ and $\ell^{2}(\Lambda)=\mathds{R}^{\Lambda}$ for a finite index set $\Lambda$.

3.1 Proposed modular aNETT training

To find a suitable network $\mathbf{D}\circ\mathbf{E}$ defining the aNETT regularizer $\mathcal{R}(x)=\mathcal{Q}(\mathbf{E}(x))+(c/2)\,\lVert x-\mathbf{D}\circ\mathbf{E}(x)\rVert^{2}$, we propose a modular data driven approach that comes in two separate steps. In a first step, we train a $\mathcal{Q}$-regularized denoising autoencoder $\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E}$ independent of the forward problem $\mathbf{K}$, whose purpose is to well represent elements of a training data set by low complexity encoder coefficients. In a second step, we train a task-specific network that increases the ability of the aNETT regularizer to distinguish between clean images and images containing problem specific artifacts.

Let $x_{1},\dots,x_{m}\in\mathcal{M}$ denote the given set of artifact-free training phantoms.

• •

$\mathcal{Q}$-regularized autoencoder:

First, an autoencoder $\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E}$ is trained such that $x_{i}$ is close to $\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E}(x_{i})$ and that $\mathcal{Q}(\mathbf{E}(x_{i}))$ is small for the given training signals. For that purpose, let $(\mathbf{D}_{\theta}\circ\mathbf{E}_{\theta})_{\theta\in\Theta}$ be a family of autoencoder networks, where $\mathbf{E}_{\theta}\colon\mathds{R}^{N\times N}\to\mathds{R}^{\Lambda}$ are encoder and $\mathbf{D}_{\theta}\colon\mathds{R}^{\Lambda}\to\mathds{R}^{N\times N}$ decoder networks, respectively.

To achieve that unperturbed images are sparsely represented by $\mathbf{E}$, whereas disrupted images are not, we apply the following training strategy. We randomly generate images $x_{i}+a_{i}\epsilon_{i}$ where $\epsilon_{i}$ is additive Gaussian white noise with a standard deviation proportional to the mean value of $x_{i}$, and $a_{i}\in\{0,1\}$ is a binary random variable that takes each value with probability $0.5$. For the numerical results below we use a standard deviation of 0.05 times the mean value of $x_{i}$. To select the particular autoencoder based on the training data, we consider the following training strategy

[TABLE]

and set $[\mathbf{D}^{\mathcal{Q}},\mathbf{E}]\coloneqq[\mathbf{D}_{\theta^{*}},\mathbf{E}_{\theta^{*}}]$. Here $\nu,\beta>0$ are regularization parameters.

Including perturbed signals $x_{i}+\epsilon_{i}$ in (3.1) increases robustness of the $\mathcal{Q}$-regularized autoencoder. To enforce regularity for the encoder coefficients only on the noise-free images, the penalty $\mathcal{Q}(\mathbf{E}(\,\cdot\,))$ is only used for the noise-free inputs, reflected by the pre-factor $1-a_{i}$. Using auto-encoders, regularity for a signal class could also be achieved by means of dimensionality reduction techniques, where $\mathds{R}^{\Lambda}$ is used as a bottleneck in the network architecture. However, in order to get a regularizer that is able to distinguish between perturbed and undisturbed perturbed signals we use $\mathds{R}^{\Lambda}$ to be of sufficiently high dimensionality.

• •

Task-specific network:

Numerical simulations showed that the $\mathcal{Q}$-regularized autoencoder alone was not able to sufficiently well distinguish between artifact-free training phantoms and images containing problem specific artifacts. In order to address this issue, we compose the operator independent network with another network $\mathbf{U}$, that is trained to distinguish between image with and without problem specific artifacts.

For that purpose, we consider randomly generated images $z_{1},\dots,z_{m}$ where either $z_{i}=(\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E})(x_{i})$ or $z_{i}=(\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E})(\mathbf{K}^{\ddagger}(\mathbf{K}x_{i}+\eta_{i}))$ with equal probability. Here $\mathbf{K}^{\ddagger}$ is an approximate right inverse and $\eta_{i}$ are error terms. We choose a network architecture $(\mathbf{U}_{\theta})_{\theta\in\Theta}$ and select $\mathbf{U}=\mathbf{U}_{\theta^{*}}$, where

[TABLE]

for some regularization parameter $\gamma>0$. In particular, the image residuals $(\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E})(\mathbf{K}^{\ddagger}(\mathbf{K}x_{i}+\eta_{i}))-(\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E})(x_{i})$ now depend on the specific inverse problem and we can consider them to consist of operator and training signal specific artifacts.

The above training procedure ensures that the network $\mathbf{U}$ adapts to the inverse problem at hand as well as to the $\mathcal{Q}$-regularized autoencoder. Training the network $\mathbf{U}$ independently of $\mathbf{D}^{\mathcal{Q}}\circ\mathbf{E}$, or directly training the auto-encoder to distinguish between images with and without problem specific artifacts, we empirically found to perform considerably worse.

The final autoencoder is then given as $\mathbf{N}=\mathbf{D}\circ\mathbf{E}$ with modular decoder $\mathbf{D}\coloneqq\mathbf{U}\circ\mathbf{D}^{\mathcal{Q}}$. For the numerical results we take $(\mathbf{U}_{\theta})_{\theta\in\Theta}$ as the tight frame U-Net of [19]. Moreover, we choose $(\mathbf{D}_{\theta}\circ\mathbf{E}_{\theta})_{\theta\in\Theta}$ as modified tight frame U-Net proposed in [35] for deep synthesis regularization. In particular, opposed to the original tight frame U-net, the modified tight frame U-Net does not involve skip connections.

3.2 Possible aNETT minimization

For minimizing the aNETT functional (1.3) we use the alternating direction method of multiplies (ADMM) with scaled dual variable [8, 14, 15]. For that purpose, the aNETT minimization problem is rewritten as the following constraint minimization problem

[TABLE]

The resulting ADMM update scheme with scaling parameter $\rho>0$ initialized by $\xi_{0}=\mathbf{E}(\mathbf{K}^{\ddagger}y^{\delta})$ and $\eta_{0}=0$ then reads as follows:

1. (S1)

$x_{k+1}=\operatorname*{arg\,min}_{x}\mathcal{D}(\mathbf{K}x,y^{\delta})+(\alpha c/2)\,\lVert x-\mathbf{N}(x)\rVert_{2}^{2}+(\rho/2)\,\lVert\mathbf{E}(x)-\xi_{k}+\eta_{k}\rVert_{2}^{2}$. 2. (S2)

$\xi_{k+1}=\operatorname*{arg\,min}_{\xi}\alpha\mathcal{Q}(\xi)+(\rho/2)\,\lVert\mathbf{E}(x_{k+1})-\xi+\eta_{k}\rVert_{2}^{2}$. 3. (S3)

\eta_{k+1}=\eta_{k+1}+\bigl{(}\mathbf{E}(x_{k+1})-\xi_{k+1}\bigr{)}.

One interesting feature of the above approach is that the signal update (S1) is independent of the possibly non-smooth penalty $\mathcal{Q}$. Moreover, the encoder update (S2) uses the proximal mapping of $\mathcal{Q}$ which in important special cases can be evaluated explicitly and therefore fast and exact. Moreover, it guarantees regular encoder coefficients during each iteration. For example, if we choose the penalty as the $\ell^{1}$-norm, then (S2) is a soft-thresholding step which results in sparse encoder coefficients. Step (S1) in typical cases has to be computed iteratively via an inner iteration. To find an approximate solution for (S1) for the results presented below we use gradient descent with at most $10$ iterations. We stop the gradient descent updates early if the difference of the functional evaluated at two consecutive iterations is below our predefined tolerance of $10^{-5}$.

The concrete implementation of the aNETT minimization requires specification of the similarity measure, the total number of outer iterations $N_{\mathrm{iter}}$, the step-size $\gamma$ for the iteration in (S1) and the parameters defining the aNETT functional. These specifications are selected dependent of the inverse problem at hand. Table 3.1 lists the particular choices for the reconstruction scenarios considered in the following section.

The ADMM scheme for aNETT minimization shares similarities with existing iterative neural network based reconstruction methods. In particular, ADMM inspired plug-and-play priors [44, 9, 38] may be most closely related. However, opposed the plug and play approach we can deduce convergence from existing results for ADMM for non-convex problems [46]. While convergence of (S1)-(S3) and relations with plug and play priors are interesting and relevant, they are beyond the scope of this work. This also applies to the comparison with other iterative minimization schemes for minimizing aNETT.

4 Application to sparse view and low dose CT

In this section we apply aNETT regularization to sparse view and low-dose computed tomography (CT). For the experiments we always choose $\mathcal{Q}$ to be the $\ell^{1}$-norm. The parameter specifications for the proposed aNETT functional and its numerical minimization are given in Table 3.1. For quantitative evaluation, we use the peak-signal-to-noise-ratio (PSNR) defined by

[TABLE]

Here $x\in\mathds{R}^{N\times N}$ is the ground truth image and $x_{\rm rec}\in\mathds{R}^{N\times N}$ its numerical reconstruction. Higher value of PSNR indicates better reconstruction.

4.1 Discretization and dataset

For sparse view CT as well as for low dose CT we work with a discretization of the Radon transform $(\mathbf{R}f)(\varphi,s)\coloneqq\int_{\mathds{R}}f(s\cos(\varphi)-t\sin(\varphi),s\sin(\varphi)+t\cos(\varphi))\mathrm{d}t$. The values $(\mathbf{R}f)(\varphi,s)$ are integrals of the function $f\colon\mathds{R}^{2}\to\mathds{R}$ over lines orthogonal to $(\cos(\varphi),\sin(\varphi))^{\intercal}$ for angle $\varphi\in[0,\pi)$ and signed distance $s\in\mathds{R}$. We discretize the Radon transform using the ODL library [1] where we assume that the function has compact support in $[-1,1]^{2}$ and sampled on an equidistant grid. We use $N_{\varphi}$ equidistant samples of $\varphi\in[0,\pi)$ and $N_{s}$ equidistant samples of $s\in[-1.5,1.5]$. In both cases, we end up with an inverse problem of the form (1.1), where $\mathbf{K}\colon\mathds{R}^{N\times N}\to\mathds{R}^{N_{\varphi}\times N_{s}}$ is the discretized linear forward operator. Elements $x\in\mathds{R}^{N\times N}$ will be referred to as CT images and the elements $y\in\mathds{R}^{N_{\varphi}\times N_{s}}$ as sinograms.

or all results presented below we work with image size $512\times 512$ and use $N_{s}=768$. The number of angular samples $N_{\varphi}$ is taken $40$ for low the dose CT and $N_{\varphi}=1138$ for the low dose example. In both cases we use the CT images from the Low Dose CT Grand Challenge dataset [32] provided by the Mayo Clinic. The dataset consists of $512\times 512$ grayscale images of $10$ different patients, where for each patient there are multiple CT scanning series available. We use the split $7/2/1$ for training, validation and testing which corresponds to $4267/1143/526$ CT images in the respective sets. We use the validation set to select networks which achieve the minimal loss on the validation set. The test set is used to evaluate the final performance. Note that by splitting the dataset according to patient we avoid validation and testing on images patients that have already be seen during training time. An example image and the corresponding simulated sparse view and low-dose data are shown in Figure 3.1.

4.2 Numerical results

We compare results of aNETT to the learned primal-dual algorithm (LPD) [3], the tight frame U-Net [19] applied as post-processing network (CNN), and the filtered back-projection (FBP) as a base-line. Minimization of the loss-function for all methods was done using Adam [22] for $100$ epochs, cosine decay learning rate $\eta_{t}=(\eta_{0}/2)\cdot(1+\cos(\pi t/100))$ with $\eta_{0}=10^{-3}$ in the $t$-th epoch, and a batch-size of $4$. For the LPD we take hyper-parameters $N_{\mathrm{primal}}=N_{\mathrm{dual}}=5$ and $N=7$ network iterations and train according to [3]. For training of the tight frame U-Net we do not follow the patch approach of [19] but instead use full images obtained with FBP as CNN inputs. Training of all the networks was done on a GTX 1080 Ti with an Intel Xeon Bronze 3104 CPU.

• •

Sparse View CT: To simulate sparse view data we evaluate the Radon transform for $N_{\varphi}=40$ directions. We generate noisy data $y^{\delta}=\mathbf{K}x+\eta^{\delta}$ by adding Gaussian white noise with standard deviation taken as $0.02$ times the mean value of $\mathbf{K}x$. We use the $\ell^{2}$-norm distance as the similarity measure. Quantitative results evaluated on the test set are shown in Table 4.1. All learning-based methods yield comparable performance in terms of PSNR and clearly outperform FBP. The reconstructions shown in Figure 4.1 indicate that aNETT reconstructions are less smooth than CNN reconstructions and less blocky than LPD reconstructions.

• •

Low Dose CT: For the low dose problem, we use a fully sampled sinogram with $N_{\varphi}=1138$ and add Poisson noise corresponding to $10^{4}$ incident photons per pixel bin. The Kullback-Leibler divergence $\mathcal{D}_{\rm KL}$ is a more appropriate than the squared $\ell^{2}$-norm distance in case of Poisson noise and the reported values and reconstructions use the Kullback-Leibler divergence as the similarity measure. Quantitative results are shown in Table 4.1. Again, all learning-based methods give similar results and significantly outperform FBP. Visual comparison of the reconstructions in Figure 4.2 shows that CNN yields cartoon like images and the LPD reconstruction again looks blocky. The aNETT reconstruction shows more texture than the CNN reconstruction and at the same time is less blocky than the LPD reconstruction.

• •

Universality: In practical applications, we may not have a fixed sampling pattern. If we have many different sampling patterns, then training network for each sampling pattern is infeasible and hence reconstruction methods should be applicable to different sampling scenarios. Additionally, it is desirable that an increased number of samples indeed increases performance. In order to test this issue, we consider the sparse view CT problem but with an increased number of angular samples without retraining the networks. LPD is not applicable in this case, as the changing the forward operator changes the network architecture. Quantitative evaluation for this scenario is given in Table 4.1. We see that aNETT are better than CNN in terms of PSNR. The advantage of aNETT over CNN, however, is best observed in Figure 4.3. One observes that CNN yields a similar reconstructions for both angular sampling patterns. On the other hand, aNETT is able to synergistically combine the increased sampling rate of the sinogram with the network trained on coarsely sampled data. Despite using the network, aNETT with 160 angular samples reconstructs small details which are not present in the reconstruction from 40 angular samples.

4.3 Discussion

The results show that the proposed aNETT regularization is competitive with prominent deep-learning methods such as LPD and post-processing CNNs. We found that the aNETT does not suffer as much from over-smoothing which is often observed in other deep-learning reconstruction methods. This can for example be seen in Figure 4.2 where the CNN yields an over-smoothed reconstruction and the aNETT reconstruction shows more texture. Besides this, aNETT reconstructions are less blocky than LPD reconstructions. Moreover, aNETT is able to leverage higher sampling rates to reconstruct small details while other deep-learning methods fail to do so. We conjecture that this is due the synergistic interplay of the aNETT regularizer with the data-consistency term in (1.3). In some scenarios, it may not be possible to retrain networks. Especially for learned iterative schemes network training is a time-consuming task. Training aNETT on the other hand is straightforward and, as demonstrated, yields a method which is robust to changes of the forward problem during testing time.

Finally, we note that aNETT relies on minimizing (1.3) iteratively. With the use of the ADMM minimization scheme presented in this article, aNETT is slower than the methods used for comparison in this article. Designing faster optimization schemes for (1.3) is beyond the scope of this work, but is an important and interesting aspect.

5 Conclusion

We have proposed the aNETT (augmented NETwork Tikhonov) for which derived coercivity of the regularizer under quite mild assumptions on the networks involved. Using this coercivity we presented a convergence analysis of aNETT with a general similarity measure $\mathcal{D}$. We proposed a modular training strategy in which we first train an $\mathcal{Q}$-regularized autoencoder independent of the problem at hand and then a network which is adapted to the problem and first autoencoder. Experimentally we found this training strategy to be superior to directly training the autoencoder on the full task. Lastly, we conducted numerical simulations demonstrating the feasibility of aNETT.

The experiments show that aNETT is able to keep up with the classical post-processing CNNs and the learned primal-dual approach for to sparse view and low dose CT. Typical deep learning methods work well for a fixed sampling pattern on which they have been trained on. However, reconstruction methods are expected to perform better if we use an increased sampling rate. We have experimentally shown that aNETT is able to leverage higher sampling rates to reconstruct small details in the images which are not visible in the other reconstructions. This universality can be advantageous in applications where one is not fixed to one sampling pattern or is not able to train a network for every sampling pattern.

Acknowledgments

D.O. and M.H. acknowledge support of the Austrian Science Fund (FWF), project P 30747-N32. The research of L.N. has been supported by the National Science Foundation (NSF) Grants DMS 1212125 and DMS 1616904.