# Invertible residual networks in the context of regularization theory for
linear inverse problems

Clemens Arndt  Center for Industrial Mathematics, University of Bremen, 28359
Bremen, Germany. Alphabetical author order. Emails: {carndt, adenker,
heilenkoetter, iskem, pmaass<EMAIL_ADDRESS>{sdittmer,
<EMAIL_ADDRESS>Alexander Denker11footnotemark: 1 Sören
Dittmer11footnotemark: 1 Nick Heilenkötter11footnotemark: 1 Meira
Iske11footnotemark: 1 Tobias Kluth11footnotemark: 1 Peter Maass11footnotemark:
1 Judith Nickel11footnotemark: 1

###### Abstract

Learned inverse problem solvers exhibit remarkable performance in applications
like image reconstruction tasks. These data-driven reconstruction methods
often follow a two-step scheme. First, one trains the often neural network-
based reconstruction scheme via a dataset. Second, one applies the scheme to
new measurements to obtain reconstructions. We follow these steps but
parameterize the reconstruction scheme with invertible residual networks
(iResNets). We demonstrate that the invertibility enables investigating the
influence of the training and architecture choices on the resulting
reconstruction scheme. For example, assuming local approximation properties of
the network, we show that these schemes become convergent regularizations. In
addition, the investigations reveal a formal link to the linear regularization
theory of linear inverse problems and provide a nonlinear spectral
regularization for particular architecture classes. On the numerical side, we
investigate the local approximation property of selected trained architectures
and present a series of experiments on the MNIST dataset that underpin and
extend our theoretical findings.

##  1 Introduction

In inverse problems, one aims to recover underlying causes from measurements
by reversing the forward measurement process. Naturally, they arise in several
fields, e.g., medical imaging, non-destructive testing, and PDE-based models.
One defines inverse problems via their ill-posed nature, i.e., reversing the
measurement process is discontinuous, and obtaining reasonable reconstructions
requires a stable reconstruction scheme. Usually, one defines the forward
problem by possibly nonlinear forward operator
$\tilde{A}\mathrel{\mathop{\mathchar 58\relax}}X\to Y$ mapping between Banach
or Hilbert spaces. Linear spectral regularizations for linear problems in
Hilbert space settings were already well studied over two decades ago [17].
Also, more general nonlinear regularizations, e.g., sophisticated variational
and iterative approaches, were the subject of extensive investigations for
linear and nonlinear problems. We refer to the review [9] for a more general
overview of this development. More recently, the linear spectral
regularization approach was generalized to diagonal frames [15]. Frames
provide larger flexibility by allowing advantageous representations of the
data, e.g., via a set of images, instead of solely singular functions provided
by the operator. Nonlinear spectral regularizations with nonlinear dependence
on the measurement data were only considered in a few works, e.g., in the
context of CG methods [17, Cor. 7.4].

During the last decade, these research directions have been accompanied by a
highly dynamic development of learning-based methods for inverse problems (see
[4] for an earlier review). Many promising methodological directions have
arisen, such as end-to-end learned reconstructions [35, 20], learned
postprocessing [21, 40], unrolled iteration schemes [18, 1, 19], plug-and-play
priors [42, 38, 24], learned penalties in variational approaches [28, 31, 26,
2, 34], (deep) generative models [10, 13], regularization by architecture [41,
14, 3], and several more directions are developed. Besides the pure method
development, an increasing number of works investigate theoretical
justifications using the framework of regularization theory (see the recent
survey [32] and the more detailed consideration in Section 1.1).

In the present work, we follow a general supervised learning approach for a
reconstruction scheme based on the concept of invertible residual networks
[8]. We use this architecture to provide a nonlinear regularization scheme for
which we develop the general convergence theory by exploiting a local
approximation property of the underlying network. Furthermore, instead of
solving the inverse problem directly, we use a training strategy that aims to
approximate the forward operator. We also introduce an architecture type
acting on the singular function directions, and show that it provides a data-
dependent nonlinear spectral regularization which is theoretically analyzed
for specific shallow architectures and it is linked to the established linear
spectral regularization theory. In addition, we underpin our theoretical
findings with several numerical experiments.

The manuscript is structured as follows: Section 2 defines the problem setting
and the theoretical framework. We then investigate general regularization
properties in Section 3. In Section 4, we discuss specific architectures and
the outcome of a training approach aiming for approximating the forward
operator. We then relate it to the classical filter-based regularization
theory. Following the theoretical investigation, Section 5 presents a series
of numerical experiments. We conclude with a discussion and outlook in Section
6.

###  1.1 Related work

Empirically the success of learning-based methods has been established for
several applications like various imaging modalities. Theoretically, there is
still a gap, e.g., regarding the convergence properties of such reconstruction
schemes. Recently an increasing number of works have examined this aspect. We
recommend the excellent survey article [32], illustrating regularization
properties of some learned methods [32, Fig.3].

One early example providing a convergent learned postprocessing approach is
the deep null space network [40]. It utilizes established regularization
methods and turns them into learned variants. While it replaces the concept of
a minimum norm solution, it inherits convergence guarantees. [5] investigates
the convergence of regularization methods trained on a finite set of training
samples without access to the entire forward operator. Specifically, the
authors consider the projections on the subspaces spanned by the dataset.
Assumptions on the dataset ensure desired approximation properties for the
convergence analysis of the regularization method. The authors then also
consider the data-dependent projected problem in a variational Tikhonov-type
framework and derive convergence of the regularization by assuming certain
approximation properties of the projection being fulfilled for the desired
solution.

Also, Tikhonov-type/variational methods and their classical regularization
theory were used to obtain learned regularizations with convergence
guarantees. In particular, learned convex regularizers [31, 33] and the
network Tikhonov (NETT) approach [26, 34] investigated this.

[16] recently suggested convergence results for plug-and-play reconstructions
that are conceptionally motivated by the Tikhonov theory. They guarantee
convergence by assuming the required conditions on their parameterized family
of denoising operators.

The deep image prior [41] proposed exploiting a network’s training and
architecture as a regularization to achieve better image representation. Based
on this concept of unsupervised learning from one single measurement, the
authors of [14] derived an analytic deep prior framework formulated as a
bilevel optimization problem. This allowed them to verify regularization
properties for particular cases by exploiting relations to the classical
filter theory of linear regularizations. More recently, [3] investigated the
equivalence between Ivanov regularization and analytic deep prior. They
investigated the inclusion of early stopping and proved regularization
properties [3].

An earlier work [12] proposed a data-driven approach to a learned linear
spectral regularization; more recently, [7, 22] considered convergence
aspects. In [22], the authors learn a scalar for each singular function
direction in a filter-based reconstruction scheme, i.e., a linear
regularization scheme. Due to assumed properties of noise and data
distributions, their training outcome is equivalent to a standard Tikhonov
regularization with a data- and noise-dependent linear operator included in
the penalty term which is diagonal with respect to the system of singular
functions. The authors further verify convergence results concerning the data-
and noise-dependent regularization parameters.

##  2 Problem setting - theoretical framework

We consider linear inverse problems based on the operator equation

$\tilde{A}x=y,$ (2.1)

where $\tilde{A}\colon X\to Y$ is a linear operator between Hilbert spaces $X$
and $Y$. For simplification $\|\tilde{A}\|=1$ is assumed, which can be easily
obtained by a scaling of the operator. We aim to recover the unknown ground
truth $x^{\dagger}$ as good as possible by only having access to a noisy
observation $y^{\delta}\in Y$ such that
$\|y^{\delta}-\tilde{A}x^{\dagger}\|\leqslant\delta$, where $\delta>0$ is the
noise level.

Instead of solving (2.1) directly, we define $A=\tilde{A}^{*}\tilde{A}$ and
$z^{\delta}=\tilde{A}^{*}y^{\delta}$ to get the operator equation

$Ax=\tilde{A}^{\ast}\tilde{A}x=\tilde{A}^{\ast}y=z$ (2.2)

which only acts on $X$, i.e., we consider the normal equation with respect to
$\tilde{A}$. We propose a two-step data-based approach for solving (2.2),
resp. (2.1), which we refer to as the iResNet reconstruction approach:

  * (I)

Supervised training of an invertible neural network $\varphi_{\theta}\colon
X\to X$ to approximate $A$.

  * (II)

Using $\varphi_{\theta}^{-1}$ (or $\varphi_{\theta}^{-1}\circ\tilde{A}^{*}$,
respectively) to solve the inverse problem.

In general, the supervised training covers both cases, either training on
noise free data tuples $(x^{(i)},Ax^{(i)})$ or on noisy tuples
$(x^{(i)},Ax^{(i)}+\eta^{(i)})$ with noise $\eta^{(i)}$ for a given dataset
$\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}\subset X$. To guarantee the
invertibility we use a residual network of the form

$\varphi_{\theta}(x)=x-f_{\theta}(x)$ (2.3)

and restrict $f_{\theta}$ to be Lipschitz continuous with
$\mathrm{Lip}(f_{\theta})\leqslant L<1$. We refer to this architecture, which
is illustrated in Figure 1, as an invertible residual network (iResNet) [8].
Note that [8] considers concatenations of several of these invertible blocks
while we focus on a single residual block.

Figure 1: Illustration of the residual network architecture of
$\varphi_{\theta}$ with residual function $f_{\theta}$ and skip connection.

The inversion of $\varphi_{\theta}$ can be done with a fixed point iteration.
To find an $\bar{x}$ which fulfills $\varphi_{\theta}(\bar{x})=z^{\delta}$ we
consider the fixed point equation

$\bar{x}=f_{\theta}(\bar{x})+z^{\delta}.$ (2.4)

Since $\mathrm{Lip}(f_{\theta})\leqslant L<1$ holds, the iteration

$x^{k+1}=f_{\theta}(x^{k})+z^{\delta}$ (2.5)

converges to $\bar{x}$ for any initial value $x^{0}\in X$ by the fixed point
theorem of Banach.

Besides, we assume the existence of a singular value decomposition (SVD)
$(u_{j},v_{j},\sigma_{j})_{j\in\mathbb{N}}$ of the operator $\tilde{A}$, i.e.,
$u_{j}\in Y$, $v_{j}\in X$, $\sigma_{j}>0$ such that

$\tilde{A}x=\sum_{j=1}^{\infty}\sigma_{j}\langle x,v_{j}\rangle u_{j}$ (2.6)

holds111In case of $\mathrm{dim}(\mathcal{R}(\tilde{A}))<\infty$, the number
of singular values is finite.. All compact operators guarantee this. With
respect to $A$ the vectors $v_{j}$ are eigenvectors and $\sigma_{j}^{2}$ are
the corresponding eigenvalues. The system $\\{v_{j}\\}_{j}\in\mathbb{N}$
builds an orthonormal basis of
$\overline{\mathcal{R}(A)}=\overline{\mathcal{R}(\tilde{A}^{\ast})}=\mathcal{N}(\tilde{A})^{\perp}=\mathcal{N}(A)^{\perp}$.
Due to the assumption $\|\tilde{A}\|=1$, it holds $\sigma_{j}\in(0,1]$.

The essential properties of the iResNets are summarized in the following
lemma.

######  Lemma 2.1 (General properties of iResNets).

Let $\varphi_{\theta},f_{\theta}\mathrel{\mathop{\mathchar
58\relax}}X\rightarrow X$, $\varphi_{\theta}(x)=x-f_{\theta}(x)$, where
$\mathrm{Lip}(f_{\theta})\leqslant L<1$. Then it holds:

  * (i)

$\mathrm{Lip}(\varphi_{\theta})\leqslant L+1$ (2.7)

  * (ii)

$\mathrm{Lip}(\varphi_{\theta}^{-1})\leqslant\frac{1}{1-L}$ (2.8)

  * (iii)

For given $z\in X$, the unique $x=\varphi_{\theta}^{-1}(z)$ is given by
$x=\lim_{k\rightarrow\infty}x^{k}$ where

$x^{k+1}=f_{\theta}(x^{k})+z$ (2.9)

for arbitrary and given $x^{0}\in X$.

###### Proof.

See [8, Section 2]. ∎

##  3 Regularization properties of the iResNet reconstruction

In this section, we discuss regularization properties in terms of well-
definedness, stability, and convergence of the general reconstruction approach
defined by the family of operators
$T_{L}=\varphi_{\theta,L}^{-1}\circ\tilde{A}^{\ast}$. Note that $L\in[0,1)$ in
the constraint of the Lipschitz constant of the residual function undertakes
the role of the regularization parameter by $L\to 1$, as will be made apparent
in this section. To highlight the essential dependence of the trained network
on $L$, we write $\varphi_{\theta,L}$ and $f_{\theta,L}$ in the following.
Please note that the set of possible network parameters and the underlying
architecture might also depend on $L$. We thus consider a network parameter
space $\Theta(L)$ depending on $L$.

######  Lemma 3.1 (Well-definedness and stability).

Let $\varphi_{\theta,L},f_{\theta,L}\mathrel{\mathop{\mathchar
58\relax}}X\rightarrow X$, $\varphi_{\theta,L}(x)=x-f_{\theta,L}(x)$, where
$\mathrm{Lip}(f_{\theta,L})\leqslant L<1$.

  * (i)

For any $L\in[0,1)$ and arbitrary $\theta\in\Theta(L)$,
$T_{L}=\varphi_{\theta,L}^{-1}\circ\tilde{A}^{\ast}\mathrel{\mathop{\mathchar
58\relax}}Y\rightarrow X$ is well-defined.

  * (ii)

For any $L\in[0,1)$ and arbitrary $\theta\in\Theta(L)$,
$T_{L}=\varphi_{\theta,L}^{-1}\circ\tilde{A}^{\ast}\mathrel{\mathop{\mathchar
58\relax}}Y\rightarrow X$ is Lipschitz-continuous.

###### Proof.

The statements (i) and (ii) follow immediately from Lemma 2.1(iii) and (ii). ∎

While well-definedness and continuity (stability) immediately follow from the
construction of the reconstruction method $T_{L}$, the convergence property
requires a specific approximation property, which we must guarantee during
training. We do so by introducing an index function that generalizes the
concept of convergence rates (see, e.g., [39]).

######  Definition 3.1.

An index function is a mapping $\psi\colon\mathbb{R}_{\geqslant
0}\to\mathbb{R}_{\geqslant 0}$, which is continuous, strictly increasing and
it holds $\psi(0)=0$.

Now we can formulate a convergence result for the reconstruction
$T_{L}(y^{\delta})$ for $\delta\to 0$ and a suitable a priori parameter choice
$L(\delta)$.

######  Theorem 3.1 (Convergence - local approximation property).

Let $x^{\dagger}\in X$ be a solution of the problem $\tilde{A}x=y$ for $y\in
R(\tilde{A})$ and $y^{\delta}\in Y$ fulfill $\|y^{\delta}-y\|\leqslant\delta$.
Furthermore, let the network parameters $\theta(L)\in\Theta(L)$ for
$L\in[0,1)$ be obtained in a way that the local approximation property

$\|\tilde{A}^{*}\tilde{A}x^{\dagger}-\varphi_{\theta(L),L}(x^{\dagger})\|=\mathcal{O}((1-L)\psi(1-L))\quad(\text{as
}L\to 1)$ (3.1)

holds for some index function $\psi$.

If $L\mathrel{\mathop{\mathchar 58\relax}}(0,\infty)\rightarrow[0,1)$ fulfills

$L(\delta)\rightarrow 1\quad\land\quad\frac{\delta}{(1-L(\delta))}\rightarrow
0\qquad\text{for }\delta\to 0,$ (3.2)

then for $x_{L(\delta)}^{\delta}\mathrel{\mathop{\mathchar
58\relax}}=T_{L(\delta)}(y^{\delta})=(\varphi_{\theta(L(\delta)),L(\delta)}^{-1}\circ\tilde{A}^{\ast})(y^{\delta})$
it holds

$\|x_{L(\delta)}^{\delta}-x^{\dagger}\|\rightarrow 0\qquad\text{for }\delta\to
0.$ (3.3)

###### Proof.

For improved readability we write $\varphi_{\delta}\mathrel{\mathop{\mathchar
58\relax}}=\varphi_{\theta(L(\delta)),L(\delta)}$ and
$f_{\delta}\mathrel{\mathop{\mathchar
58\relax}}=f_{\theta(L(\delta)),L(\delta)}$ in the remainder of the proof.
Using Lemma 2.1, it holds

$\displaystyle\|x_{L(\delta)}^{\delta}-x^{\dagger}\|$
$\displaystyle\leqslant\|\varphi_{\delta}^{-1}(\tilde{A}^{\ast}y^{\delta})-\varphi_{\delta}^{-1}(\tilde{A}^{\ast}y)\|+\|\varphi_{\delta}^{-1}(\tilde{A}^{\ast}y)-x^{\dagger}\|$
$\displaystyle\leqslant\frac{\|\tilde{A}^{\ast}\|}{1-L(\delta)}\|y^{\delta}-y\|+\|\varphi_{\delta}^{-1}(\tilde{A}^{\ast}\tilde{A}x^{\dagger})-x^{\dagger}\|$
$\displaystyle\leqslant\frac{\delta}{1-L(\delta)}+\|\varphi_{\delta}^{-1}(\tilde{A}^{\ast}\tilde{A}x^{\dagger})-\varphi_{\delta}^{-1}(\varphi_{\delta}(x^{\dagger}))\|$
$\displaystyle\leqslant\frac{\delta}{1-L(\delta)}+\frac{1}{1-L(\delta)}\|\tilde{A}^{\ast}\tilde{A}x^{\dagger}-\varphi_{\delta}(x^{\dagger})\|.$
(3.4)

The assertion follows due to the assumptions (3.1) and (3.2) as $\lim_{L\to
1}\psi(1-L)=0$. ∎

######  Remark 3.1.

We can obtain a convergence rate of
$\|x_{L(\delta)}^{\delta}-x^{\dagger}\|=\mathcal{O}\big{(}\delta^{\,\varepsilon/(1+\varepsilon)}\big{)}$
for $\varepsilon>0$, if (3.1) is fulfilled with the index function
$\psi(\lambda)=\lambda^{\varepsilon}$, i.e.,

$\|\tilde{A}^{*}\tilde{A}x^{\dagger}-\varphi_{\theta(L),L}(x^{\dagger})\|=\mathcal{O}((1-L)^{1+\epsilon})$
(3.5)

and by choosing $1-L(\delta)\sim\delta^{1/(1+\varepsilon)}.$

######  Remark 3.2.

Note that the local approximation property (3.1) is a slightly stronger
assumption, which is derived by exploiting the properties of the iResNet to
remove the explicit dependency on the inverse within the constraint. A weaker
but sufficient condition (which is implied by (3.1)) is

$\|\varphi_{\theta(L),L}^{-1}(\tilde{A}^{\ast}\tilde{A}x^{\dagger})-x^{\dagger}\|\rightarrow
0\quad\text{as }L\to 1.$ (3.6)

The local approximation property (3.1) provides a relation between the
approximation capabilities of the network architecture and the regularization
properties of the iResNet reconstruction scheme, which also contains a certain
kind of source condition. Of particular interest for the convergence
properties is the approximation property set

$S=\\{x\in X\,|\,\exists\text{ index function }\psi\mathrel{\mathop{\mathchar
58\relax}}\,x\text{ fulfills }\eqref{eq:conv_trainingproperty_ass}\\}.$ (3.7)

There are no specific assumptions on the design of the training of the network
$\varphi_{\theta}$ made by (3.1). But since both architecture choice and
training procedure are crucial for the resulting network parameters
$\theta(L)$, they also have a strong impact on the local approximation
property and the structure of $S$. This is discussed in more detail in the
context of the specific shallow architectures in Section 4. In addition, we
report numerical experiments regarding the local approximation property in
Section 5.

A simple example illustrating the conditions under which a linear network
satisfies the local approximation property is given in the following remark.

######  Remark 3.3 (Local approximation property for linear networks).

Assume that the network $\varphi_{\theta,L}$ is linear and the prediction
error is bounded by
$\|\varphi_{\theta(L(\delta)),L(\delta)}(x^{(i)})-\tilde{A}^{*}y^{\delta,(i)}\|\leqslant\zeta(\delta)$
with $\zeta(\delta)\in\mathbb{R}$ on a dataset
$(x^{(i)},y^{\delta,(i)})_{i=1,\dots,N}$ with
$\|\tilde{A}x^{(i)}-y^{\delta,(i)}\|\leqslant\delta$. Then, for $x^{\dagger}$
given by $x^{\dagger}=\sum_{i=1}^{N}t_{i}x^{(i)}$ with $t_{i}\in\mathbb{R}$ we
have

$\displaystyle\|\tilde{A}^{*}\tilde{A}x^{\dagger}-\varphi_{\theta(L(\delta)),L(\delta)}(x^{\dagger})\|$
$\displaystyle=\Big{\|}\tilde{A}^{\ast}\tilde{A}\sum_{i=1}^{N}t_{i}x^{(i)}-\sum_{i=1}^{N}t_{i}\varphi_{\theta(L(\delta)),L(\delta)}x^{(i)}\Big{\|}$
$\displaystyle\leqslant\sum_{i=1}^{N}|t_{i}|\,\|\tilde{A}^{\ast}\tilde{A}x^{(i)}-\varphi_{\theta(L(\delta)),L(\delta)}x^{(i)}\|$
$\displaystyle\leqslant\sum_{i=1}^{N}|t_{i}|\left(\|\tilde{A}^{\ast}\tilde{A}x^{(i)}-\tilde{A}^{*}y^{\delta,(i)}\|+\|\tilde{A}^{*}y^{\delta,(i)}-\varphi_{\theta(L(\delta)),L(\delta)}x^{(i)}\|\right)$
$\displaystyle\leqslant\sum_{i=1}^{N}|t_{i}|\,(\delta+\zeta(\delta)).$ (3.8)

As as result, we obtain convergence of $x_{L(\delta)}^{\delta}$ to
$x^{\dagger}$ if

$\frac{\delta+\zeta(\delta)}{1-L(\delta)}\to 0\quad\text{for }\delta\to 0$
(3.9)

and if $L(\delta)$ satisfies the conditions of (3.2).

Similar to a source condition, (3.1) also has an influence on the type of
solution $x^{\dagger}$ one approximates in the limit $\delta\to 0$ if the
operator $\tilde{A}$ has a non-trivial kernel. For example, if
$\mathcal{R}(f_{\theta(L),L})=\mathcal{N}(\tilde{A})^{\perp}$, the local
approximation property enforces that
$P_{\mathcal{N}(\tilde{A})}x^{\dagger}=0$, i.e., together with the
approximation of $\tilde{A}^{\ast}\tilde{A}$ one would consider the minimum
norm solution. More general, we can formulate the following property of $S$.

######  Lemma 3.2.

Let $x_{1}^{\dagger},x_{2}^{\dagger}\in S$ both be solutions of $\tilde{A}x=y$
for one $y\in\mathcal{R}(\tilde{A})$. Then it holds
$x_{1}^{\dagger}=x_{2}^{\dagger}$.

###### Proof.

For abbreviation, we write $f_{L}$ and $\varphi_{L}$ instead of
$f_{\theta(L),L}$ and $\varphi_{\theta(L),L}$. Using
$\tilde{A}x_{1}^{\dagger}=\tilde{A}x_{2}^{\dagger}$, it holds

$\displaystyle\quad\,\,\|x_{1}^{\dagger}-x_{2}^{\dagger}\|$
$\displaystyle\leqslant\|x_{1}^{\dagger}-\tilde{A}^{*}\tilde{A}x_{1}^{\dagger}-f_{L}(x_{1}^{\dagger})\|+\|\tilde{A}^{*}\tilde{A}x_{1}^{\dagger}+f_{L}(x_{1}^{\dagger})-\tilde{A}^{*}\tilde{A}x_{2}^{\dagger}-f_{L}(x_{2}^{\dagger})\|+\|\tilde{A}^{*}\tilde{A}x_{2}^{\dagger}+f_{L}(x_{2}^{\dagger})-x_{2}^{\dagger}\|$
$\displaystyle=\|\varphi_{L}(x_{1}^{\dagger})-\tilde{A}^{*}\tilde{A}x_{1}^{\dagger}\|+\|f_{L}(x_{1}^{\dagger})-f_{L}(x_{2}^{\dagger})\|+\|\tilde{A}^{*}\tilde{A}x_{2}^{\dagger}-\varphi_{L}(x_{2}^{\dagger})\|$
$\displaystyle\leqslant\|\tilde{A}^{*}\tilde{A}x_{1}^{\dagger}-\varphi_{L}(x_{1}^{\dagger})\|+L\|x_{1}^{\dagger}-x_{2}^{\dagger}\|+\|\tilde{A}^{*}\tilde{A}x_{2}^{\dagger}-\varphi_{L}(x_{2}^{\dagger})\|.$
(3.10)

Subtracting $L\|x_{1}^{\dagger}-x_{2}^{\dagger}\|$, it follows

$(1-L)\|x_{1}^{\dagger}-x_{2}^{\dagger}\|\leqslant\|\tilde{A}^{*}\tilde{A}x_{1}^{\dagger}-\varphi_{L}(x_{1}^{\dagger})\|+\|\tilde{A}^{*}\tilde{A}x_{2}^{\dagger}-\varphi_{L}(x_{2}^{\dagger})\|.$
(3.11)

Since $x_{1}^{\dagger},x_{2}^{\dagger}$ both fulfill the local approximation
property (3.1), there exists an index function $\psi$ such that
$(1-L)\|x_{1}^{\dagger}-x_{2}^{\dagger}\|=\mathcal{O}((1-L)\psi(1-L))$ must
hold. This implies
$\|x_{1}^{\dagger}-x_{2}^{\dagger}\|=\mathcal{O}(\psi(1-L))$, which is only
possible for $x_{1}^{\dagger}=x_{2}^{\dagger}$. ∎

###  3.1 Diagonal architecture

We continue with the introduction of a particular network design being
exploited in the remainder of the manuscript. In order to investigate
relations to established and well-studied spectral regularization methods we
consider a specific network architecture acting on the coefficients with
respect to the singular functions $v_{j}$ of the operator $\tilde{A}$. We
assume that the residual function is given by subnetworks
$f_{\theta,j}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\rightarrow\mathbb{R}$, $j\in\mathbb{N}$, i.e.,

$f_{\theta}(x)=\sum_{j\in\mathbb{N}}f_{\theta,j}(\langle
x,v_{j}\rangle)v_{j}.$ (3.12)

We refer to this architecture as diagonal architecture which is also
illustrated in Figure 2. The naming is motivated by the similarity of the
structure of $f_{\theta}$ to a matrix multiplication with a diagonal matrix
w.r.t the basis $\\{v_{j}\\}_{j\in\mathbb{N}}$. In this case, the Lipschitz
constraint $\mathrm{Lip}(f_{\theta})\leqslant L$ is fulfilled if
$\mathrm{Lip}(f_{\theta,j})\leqslant L$ holds for all $j\in\mathbb{N}$, i.e.,
for any $x,y\in X$ it holds

$\|f_{\theta}(x)-f_{\theta}(y)\|^{2}=\sum_{j\in\mathbb{N}}|f_{\theta,j}(\langle
x,v_{j}\rangle)-f_{\theta,j}(\langle y,v_{j}\rangle)|^{2}\leqslant
L^{2}\sum_{j\in\mathbb{N}}|\langle x-y,v_{j}\rangle|^{2}\leqslant
L^{2}\|x-y\|^{2}.$ (3.13)

The local approximation property (3.1), in this case for
$x^{\dagger}\in\mathcal{N}(A)^{\perp}$, reads

$|f_{\theta(L),L,j}(x_{j}^{\dagger})-(1-\sigma_{j}^{2})x_{j}^{\dagger}|=\mathcal{O}((1-L)\psi(1-L))$
(3.14)

for any $j\in\mathbb{N}$, where $x_{j}^{\dagger}=\langle
x^{\dagger},v_{j}\rangle$.

The particular aim of using the diagonal architecture is now to find a data-
dependent filter function
$r_{L}\colon\mathbb{R}_{+}\times\mathbb{R}\to\mathbb{R}$ for self-adjoint
forward operators such that for $z\in\mathcal{R}(\tilde{A}^{\ast})$

$\displaystyle\varphi_{\theta}^{-1}(z)=\sum_{j\in\mathbb{N}}r_{L}(\sigma_{j}^{2},\langle
z,v_{j}\rangle)\langle z,v_{j}\rangle v_{j},$ (3.15)

or

$T_{L}(y)=\varphi_{\theta}^{-1}(\tilde{A}^{\ast}y)=\sum_{j\in\mathbb{N}}r_{L}(\sigma_{j}^{2},\sigma_{j}\langle
y,u_{j}\rangle)\sigma_{j}\langle y,u_{j}\rangle v_{j}.$ (3.16)

This can be seen as a nonlinear extension to the established filter theory of
linear regularization methods [37, 27], where $r_{L}$ would depend on
$\sigma_{j}^{2}$ only. The subnetworks thus play an important role in defining
the filter functions via

$(Id-f_{\theta,L,j})^{-1}(s)=r_{L}(\sigma_{j}^{2},s)s.$ (3.17)

Inputs$x_{1}$$x_{2}$$x_{j}$Outputs$f_{\theta,1}(x_{1})$ $f_{\theta,2}(x_{2})$
$f_{\theta,j}(x_{j})$ Figure 2: Diagonal structure of the residual function
including subnetworks $f_{\theta,j}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\to\mathbb{R}$ acting on $x_{j}=\langle x,v_{j}\rangle$.

######  Remark 3.4.

Some authors define filter functions in terms of the original problem with
respect to the generalized inverse $\tilde{A}^{\dagger}$ such that

$\displaystyle\varphi_{\theta}^{-1}(\tilde{A}^{\ast}y)=\sum_{j\in\mathbb{N}}F_{L}(\sigma_{j},\langle
y,u_{j}\rangle)\frac{1}{\sigma_{j}}\langle y,u_{j}\rangle
v_{j},\qquad\text{with }F_{L}(\sigma,s)=\sigma^{2}r_{L}(\sigma^{2},\sigma s).$
(3.18)

Note that the choice of the representation, either in terms of $r_{L}$ or in
terms of $F_{L}$ depends on personal preference (compare for example [37] and
[27]). In the following, we represent the filter function in terms of $r_{L}$.

One simple nonlinear dependence on $s$ will be taken into account explicitly
in the following, as it becomes relevant in later investigations. In analogy
to the bias in neural network architectures, we consider a linear filter
framework with bias, i.e.,

$T_{L}(y)=\hat{b}_{L}+\sum_{j\in\mathbb{N}}\hat{r}_{L}(\sigma_{j}^{2})\sigma_{j}\langle
y,u_{j}\rangle$ (3.19)

where $\hat{b}_{L}\in X$ is a bias term and $\hat{r}_{L}$ is a classical
filter function. With additional simple assumptions on $\hat{b}_{L}$ this
becomes a regularization method. For the sake of completeness, we provide the
following result.

######  Lemma 3.3 (Filter-based spectral regularization with bias).

Let $\hat{r}_{L}\mathrel{\mathop{\mathchar
58\relax}}[0,\|\tilde{A}\|^{2}]\to\mathbb{R}$ be a piecewise continuous
function and let $b_{L}\in X$ for $L\in[0,1)$. Furthermore, let

  * (i)

$\lim_{L\to 1}\hat{r}_{L}(\sigma_{j}^{2})=\frac{1}{\sigma_{j}^{2}}$ for any
$\sigma_{j}$, $j\in\mathbb{N}$,

  * (ii)

$\exists\,0<C<\infty\mathrel{\mathop{\mathchar 58\relax}}\
\sigma_{j}^{2}\,|\hat{r}_{L}(\sigma_{j}^{2})|\leqslant C$ for any
$\sigma_{j}$, $j\in\mathbb{N}$ and $L\in[0,1)$,

  * (iii)

$\hat{b}_{L}\in X$ for $L\in[0,1)$ and $\lim_{L\to 1}\|\hat{b}_{L}\|=0$

hold. Let $x^{\dagger}\in\mathcal{N}(\tilde{A})^{\perp}$ be a solution of the
problem $\tilde{A}x=y$ for $y\in R(\tilde{A})$, the operator
$T_{L}\mathrel{\mathop{\mathchar 58\relax}}Y\rightarrow X$ be given by (3.19)
and $y^{\delta}\in Y$ with $\|y-y^{\delta}\|\leqslant\delta$. In addition, let
$L\mathrel{\mathop{\mathchar 58\relax}}(0,\infty)\rightarrow[0,1)$ be such
that

$L(\delta)\rightarrow
1\quad\land\quad\delta\sqrt{\sup\big{\\{}|\hat{r}_{L(\delta)}(\sigma_{i})|\,\big{|}\,i\in\mathbb{N}\big{\\}}}\rightarrow
0$ (3.20)

as $\delta\to 0$. Then, $T_{L}$ is a regularization method and for
$x_{L(\delta)}^{\delta}\mathrel{\mathop{\mathchar
58\relax}}=T_{L(\delta)}(y^{\delta})$ it holds $\lim_{\delta\to
0}\|x_{L(\delta)}^{\delta}-x^{\dagger}\|=0$.

###### Proof.

By defining $\tilde{x}_{L(\delta)}^{\delta}\mathrel{\mathop{\mathchar
58\relax}}=T_{L(\delta)}(y^{\delta})-\hat{b}_{L(\delta)}$ (removing the bias),
we get an ordinary linear regularization method (see, e.g., [37, Corr.
3.3.4]). It holds

$\|x_{L(\delta)}^{\delta}-x^{\dagger}\|\leqslant\|x_{L(\delta)}^{\delta}-\tilde{x}_{L(\delta)}^{\delta}\|+\|\tilde{x}_{L(\delta)}^{\delta}-x^{\dagger}\|=\|\hat{b}_{L(\delta)}\|+\|\tilde{x}_{L(\delta)}^{\delta}-x^{\dagger}\|.$
(3.21)

Since $\|\hat{b}_{L(\delta)}\|\to 0$ by assumption and
$\tilde{x}_{L(\delta)}^{\delta}\to x^{\dagger}$ by standard theory, $T_{L}$ is
a regularization method. ∎

##  4 Approximation training of specialized architectures

Besides the architecture choice and the available data, the loss function for
network training is an important ingredient. Having regularization properties
in mind, a natural choice for training the iResNet is to approximate the
forward operator $A\colon X\to X$, $A=\tilde{A}^{*}\tilde{A}$, on a given
training dataset $\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}\subset X$. This is also
strongly motivated by the structure of the local approximation property (3.1)
to obtain convergence guarantees. The training of $\varphi_{\theta}=Id-
f_{\theta}$ then consists of solving

$\min_{\theta\in\Theta}l(\varphi_{\theta},A)=\min_{\theta\in\Theta}\frac{1}{N}\sum_{i}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}\quad\text{s.t.
}\mathrm{Lip}(f_{\theta})\leqslant L$ (4.1)

with $L<1$ as a hyperparameter, which we refer to as approximation training.
Here, we restrict ourselves to the noise-free case as the outcome of the
approximation training is independent for $N\to\infty$ (assuming stochastic
independence between data and noise).

In the following, we analyze to which extent $\varphi_{\theta}^{-1}$ acts as a
regularized inverse of $A$ for different simple diagonal iResNet architectures
trained according to (4.1). We compute the particular network parameters
$\theta$ and derive the corresponding data-dependent filter function
$r_{L}\colon\mathbb{R}\times\mathbb{R}\to\mathbb{R}$ (see Section 3.1).
Additionally, we check whether the local approximation property (3.1) is
fulfilled. For this purpose, the assumptions on the architecture as well as on
the dataset used for the training must be specified.

###  4.1 One-parameter-network – Tikhonov

We begin with a very simple linear network which only has one single scalar
learnable parameter. We show that this architecture is equivalent to Tikhonov
regularization if trained appropriately.

######  Lemma 4.1.

Let $(v_{j},\sigma_{j}^{2})_{j}$ be the eigenvectors and eigenvalues of $A$
and let $\varphi_{\theta}=Id-f_{\theta}$ be an iResNet which solves (4.1) with

  * (i)

$f_{\theta}=k(Id-A)$, where $\theta=k$ and $\Theta=\mathbb{R}$ (architecture
assumption),

  * (ii)

the training dataset $\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}$ contains at least
one $x^{(i)}$ s.t. $Ax^{(i)}\neq x^{(i)}$ (dataset assumption).

Then, the solution of (4.1) is $k=L$ and (3.15) holds with

$r_{L}(\sigma^{2},s)=\frac{1}{L}\cdot\frac{1}{\alpha+\sigma^{2}},\qquad\alpha=\frac{1-L}{L},$
(4.2)

for any $s\in\mathbb{R}$.

The proof of the lemma can be found in the Appendix A.1. The derived filter
function is illustrated in Figure 3. We plot $\sigma^{2}r_{L}(\sigma^{2},s)$
since the multiplication of the filter function with $\sigma^{2}$ corresponds
to the concatenation $\varphi_{\theta}^{-1}\circ A$ which emphasizes the
regularizing effect.

$\sigma^{2}$$\sigma^{2}r_{L}(\sigma^{2},s)$00.51$0$$0.5$$1$$L=0.9$$L=0.98$
Figure 3: The plot depicts graphs of $\sigma^{2}r_{L}(\sigma^{2},s)$ as given
in (4.2) to visualize the effect of the regularization, which is also known as
standard Tikhonov regularization.

######  Remark 4.1.

This particular choice of the residual function $f_{\theta}$ is one example
where the local approximation property is not fulfilled for any
$x^{\dagger}\in X$ except for eigenvectors corresponding to
$\sigma_{1}^{2}=1$, i.e., $S=\mathrm{span}\\{v_{1}\\}$. But the linear nature
and the specific spectral representation of the residual network allows for
the verification of the weaker constraint in (3.6) being highlighted in Remark
3.2 and being sufficient for convergence. For the residual function
$f_{\theta(L),L}=L(Id-A)$ and $x^{\dagger}\in\mathcal{N}(A)^{\bot}$ we obtain

$\displaystyle\|\varphi_{\theta(L),L}^{-1}(Ax^{\dagger})-x^{\dagger}\|^{2}$
$\displaystyle=\|((Id-L(Id-A))^{-1}A-Id)x^{\dagger}\|^{2}$
$\displaystyle=\sum_{j\in\mathbb{N}}\bigg{(}\underset{=\sigma_{j}^{2}r_{L}(\sigma_{j}^{2})}{\underbrace{\frac{\sigma_{j}^{2}}{1-L(1-\sigma_{j}^{2})}}}-1\bigg{)}^{2}|\langle
x^{\dagger},v_{j}\rangle|^{2},$ (4.3)

which converges to zero for $L\to 1$.

Alternatively, convergence can be verified by standard arguments following,
for example, the line of reasoning in the proof of [27, Thm. 3.3.3], where
properties of the filter function
$F_{L}(\sigma_{j})=\sigma_{j}^{2}r_{L}(\sigma_{j}^{2})$ are exploited. Except
for the factor $\frac{1}{L}$, $F_{L}$ defines the filter function for standard
Tikhonov regularization, for which the desired properties are fulfilled such
that we have convergence for arbitrary $x^{\dagger}\in\mathcal{N}(A)^{\perp}$.

###  4.2 Affine linear network – squared soft TSVD

We continue with a slightly more general affine linear architecture
$f_{\theta}(x)=Wx+b$. It can be observed that the linear operator
$W\mathrel{\mathop{\mathchar 58\relax}}X\to X$ only depends on $A$ and $L$
while the bias vector $b\in X$ is data dependent.

######  Lemma 4.2.

Let $(v_{j},\sigma_{j}^{2})_{j}$ be the eigenvectors and eigenvalues of $A$
and $\varphi_{\theta}=Id-f_{\theta}$ be an iResNet which solves (4.1) where

  * (i)

$f_{\theta}(x)=Wx+b$, with $\theta=(W,b)$ and

$\Theta=\left\\{(W,b)\in L(X)\times
X\,\Bigg{|}\,\exists(w_{j})_{j\in\mathbb{N}},(b_{j})_{j\in\mathbb{N}}\mathrel{\mathop{\mathchar
58\relax}}\ W=\sum_{j\in\mathbb{N}}w_{j}\langle\cdot,v_{j}\rangle
v_{j},\,b=\sum_{j\in\mathbb{N}}b_{j}v_{j}\right\\},$ (4.4)

  * (ii)

the training dataset $\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}$ and the mean
values $\mu_{j}=\frac{1}{N}\sum_{i=1}^{N}\langle x^{(i)},v_{j}\rangle$ fulfill

$\forall j\in\mathbb{N}\colon\,\exists i\in\\{1,...,N\\}\colon\quad\langle
x^{(i)},v_{j}\rangle\neq\mu_{j}.$ (4.5)

Then, the solution of the training problem (4.1) is
$w_{j}=\min\\{1-\sigma_{j}^{2},L\\}$,
$b_{j}=\max\\{0,1-L-\sigma^{2}\\}\mu_{j}$ and
$\varphi_{\theta}^{-1}\circ\tilde{A^{*}}$ is equivalent to $T_{L}$ in (3.19)
with

$\hat{r}_{L}(\sigma^{2})=\frac{1}{\max\\{\sigma^{2},1-L\\}},\qquad\hat{b}_{L}=\frac{1}{1-L}\sum_{j\in\mathbb{N}}b_{j}v_{j}=\sum_{\sigma_{j}^{2}<1-L}\frac{1-L-\sigma^{2}}{1-L}\mu_{j}v_{j}.$
(4.6)

The proof of the lemma can be found in Appendix A.2. The filter function is
illustrated in Figure 4.

######  Remark 4.2.

Note that a similar filter function has been found in the context of an
analytic deep prior approach for regularization [14]. The authors derived the
filter function $F_{L}(\sigma)=\sigma^{2}\hat{r}_{L}(\sigma^{2})$, which is
linear for small $\sigma$, and named it soft TSVD (soft truncated SVD). In
contrast, the filter function in (4.6) is a polynomial of degree two for small
$\sigma$, which we refer to as squared soft TSVD.

$\sigma^{2}$$\sigma^{2}\hat{r}_{L}(\sigma^{2})$00.51$0$$0.5$$1$$L=0.8$$L=0.9$
Figure 4: The plot depicts graphs of $\sigma^{2}\hat{r}_{L}(\sigma^{2})$ as
given in (4.6) to visualize the effect of the regularization, which is similar
to the known soft TSVD regularization.

Due to the affine linear nature of the found reconstruction scheme,
regularization properties can immediately be derived by verifying certain
properties of the filter function $\hat{r}_{L}$ and the bias $\hat{b}_{L}$
using Lemma 3.3. This is summarized in the following corollary.

######  Corollary 4.1.

The filter based method (3.19) with $\hat{r}_{L}$ and $\hat{b}_{L}$ given by
(4.6) fulfills the properties (i)-(iii) in Lemma 3.3, i.e., in the setting of
Lemma 3.3 it is a convergent regularization.

The previous result provides a link to the classical theory for linear
regularization operators. Besides the well-known classical theory, we further
want to investigate the local approximation property (3.1) to verify whether
the convergence property can also be obtained using Theorem 3.1. In a first
step, we derive a technical source condition for $x^{\dagger}$ which implies
the local approximation property.

######  Lemma 4.3.

Let the network $\varphi_{\theta}$ and the setting be that of Lemma 4.2.
Assume $x^{\dagger}\in\mathcal{N}(A)^{\bot}$ and let $\psi$ be an index
function such that

$\exists\bar{\beta}\in(0,1]\mathrel{\mathop{\mathchar
58\relax}}\forall\beta\in(0,\bar{\beta})\mathrel{\mathop{\mathchar
58\relax}}\sum_{j\mathrel{\mathop{\mathchar 58\relax}}\
\sigma_{j}^{2}\leqslant\beta}\langle
x^{\dagger},v_{j}\rangle^{2}=\mathcal{O}(\psi(\beta)^{2})$ (4.7)

holds. Then there exists another index function $\tilde{\psi}$ such that the
local approximation property (3.1) is fulfilled for $x^{\dagger}$.

###### Proof.

With the setting of Lemma 4.2 and the trained parameters $(W,b)$ of the
network $f_{\theta}$, it holds

$\|\varphi_{\theta(L)}(x^{\dagger})-Ax^{\dagger}\|=\|(Id-W)x^{\dagger}-b-Ax^{\dagger}\|\leqslant\|(Id-
W-A)x^{\dagger}\|+\|b\|.$ (4.8)

Now, we estimate the convergence order of both terms w.r.t $(1-L)\to 0$.

For the first part, we exploit the simple computations

$\displaystyle\|(Id-W-A)x^{\dagger}\|^{2}$
$\displaystyle=\left\|\sum_{j}(1-w_{j}-\sigma_{j}^{2})\langle
x^{\dagger},v_{j}\rangle v_{j}\right\|^{2}$
$\displaystyle=\sum_{j}(1-\sigma_{j}^{2}-\min\\{1-\sigma_{j}^{2},L\\})^{2}\langle
x^{\dagger},v_{j}\rangle^{2}$
$\displaystyle=\sum_{j}\max\\{1-L-\sigma_{j}^{2},0\\}^{2}\langle
x^{\dagger},v_{j}\rangle^{2}$ $\displaystyle=\sum_{\sigma_{j}^{2}\leqslant
1-L}(1-L-\sigma_{j}^{2})^{2}\langle x^{\dagger},v_{j}\rangle^{2}$
$\displaystyle\leqslant\sum_{\sigma_{j}^{2}\leqslant 1-L}(1-L)^{2}\langle
x^{\dagger},v_{j}\rangle^{2}$
$\displaystyle=(1-L)^{2}\sum_{\sigma_{j}^{2}\leqslant 1-L}\langle
x^{\dagger},v_{j}\rangle^{2}.$ (4.9)

Using assumption (4.7), we obtain

$(1-L)^{2}\sum_{\sigma_{j}^{2}\leqslant 1-L}\langle
x^{\dagger},v_{j}\rangle^{2}=\mathcal{O}\left((1-L)^{2}\psi(1-L)^{2}\right).$
(4.10)

Regarding the second part, it holds

$\displaystyle\|b\|^{2}$
$\displaystyle=\sum_{j}b_{j}^{2}=\sum_{j}\max\\{0,1-L-\sigma^{2}\\}^{2}\mu_{j}^{2}$
$\displaystyle=\sum_{\sigma_{j}^{2}\leqslant
1-L}(1-L-\sigma_{j}^{2})^{2}\mu_{j}^{2}\leqslant(1-L)^{2}\sum_{\sigma_{j}^{2}\leqslant
1-L}\mu_{j}^{2}.$ (4.11)

Since the values $\mu_{j}$ form by definition (see Lemma 4.2) an $\ell^{2}$
sequence, there exists an index function $\tilde{\psi}\geqslant\psi$ such that

$\sum_{\sigma_{j}^{2}\leqslant
1-L}\mu_{j}^{2}=\mathcal{O}\left(\tilde{\psi}(1-L)^{2}\right).$ (4.12)

Overall, we obtain

$\|(Id-
W-A)x^{\dagger}\|+\|b\|=\mathcal{O}\left((1-L)^{2}\tilde{\psi}(1-L)^{2}\right),$
(4.13)

thus, the local approximation property (3.1) is fulfilled. ∎

While (4.7) is a quite technical source condition, it is implied by standard
source conditions of the form

$\exists w\in X,\mu>0\colon\quad x^{\dagger}=A^{\mu}w.$ (4.14)

In this case, the index function $\psi$ is of the form
$\psi(\beta)=\beta^{\mu}$. The proof of the following corollary reveals the
exact relation between standard source conditions and (4.7).

######  Corollary 4.2.

We assume that the setting of Lemma 4.2 holds. Let $\tilde{A}$ be compact.
Then, for any $x^{\dagger}\in\mathcal{N}(\tilde{A})^{\bot}$ the local
approximation property (3.1) is fulfilled (i.e.,
$S=\mathcal{N}(\tilde{A})^{\bot}$).

###### Proof.

Let $x^{\dagger}\in\mathcal{N}(\tilde{A})^{\bot}$ be arbitrary. By
construction, $A=\tilde{A}^{\ast}\tilde{A}$ is self-adjoint and non-negative.
As $A$ is also compact, we can deduce from [29, Thm. 1]: For any
$\varepsilon>0$ there exists an index function $\psi\mathrel{\mathop{\mathchar
58\relax}}[0,1]\to[0,\infty)$ such that $x^{\dagger}\in\\{x\in
X\,|\,x=\psi(A)w,\|w\|\leqslant(1+\varepsilon)\|x^{\dagger}\|\\}$.

This implies that there exists a $w\in X$ such that $x^{\dagger}=\psi(A)w$. We
thus obtain for $\beta>0$

$\displaystyle\sum_{\sigma_{j}^{2}\leqslant\beta}\langle
x^{\dagger},v_{j}\rangle^{2}$
$\displaystyle=\sum_{\sigma_{j}^{2}\leqslant\beta}\langle\psi(A)w,v_{j}\rangle^{2}=\sum_{\sigma_{j}^{2}\leqslant\beta}\psi(\sigma_{j}^{2})^{2}\langle
w,v_{j}\rangle^{2}$
$\displaystyle\leqslant\sum_{\sigma_{j}^{2}\leqslant\beta}\psi(\beta)^{2}\langle
w,v_{j}\rangle^{2}\leqslant\|w\|^{2}\psi(\beta)^{2}\leqslant(1+\varepsilon)^{2}\|x^{\dagger}\|^{2}\psi(\beta)^{2}.$
(4.15)

By Lemma 4.3 we get the desired local approximation property (3.1). ∎

The verification of regularization properties by Theorem 3.1 in terms of the
local approximation property is thus in line with the existing theory
exploited in Lemma 3.3 and it further illustrates the character of the local
approximation property combining the approximation capabilities of the
residual network and a source condition on the desired solution.

###  4.3 Linear network with ReLU activation

In the following, we include selected nonlinearities as activation functions
in a shallow network which allow for an analytical investigation of the
resulting reconstruction scheme. Here, we start with a ReLU activation, which
requires a certain assumption on the training dataset depending on the chosen
nonlinearity. For simplicity, we do not include a bias in the architecture, in
contrast to the architecture from the last section.

######  Lemma 4.4.

Let $(v_{j},\sigma_{j}^{2})_{j}$ be the eigenvectors and eigenvalues of $A$
and $\varphi_{\theta}=Id-f_{\theta}$ be an iResNet which solves (4.1) with

  * (i)

$f_{\theta}(x)=\phi(Wx)$, where $\theta=W$ and $\Theta=\\{W\in
L(X)\,|\,\exists\mathrel{\mathop{\mathchar
58\relax}}(w_{j})_{j\in\mathbb{N}}\mathrel{\mathop{\mathchar
58\relax}}W=\sum_{j\in\mathbb{N}}w_{j}\langle\cdot,v_{j}\rangle v_{j}\\}$ and
$\phi$ being the ReLU function w.r.t. the eigenvectors, i.e.,
$\phi(x)=\sum_{j\in\mathbb{N}}\max(0,\langle v_{j},x\rangle)v_{j}$,

  * (ii)

for every eigenvector $v_{j}$ of $A$, the training dataset
$\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}$ contains at least one $x^{(i)}$ s.t.
$\langle x^{(i)},v_{j}\rangle>0$.

Then, the solution of (4.1) is
$W=\sum_{j\in\mathbb{N}}w_{j}\langle\cdot,v_{j}\rangle v_{j}$,
$w_{j}=\min\\{1-\sigma_{j}^{2},L\\}$ and (3.15) holds with

$r_{L}(\sigma^{2},s)=\begin{cases}\frac{1}{\max\\{\sigma^{2},1-L\\}}&\text{if
}s\geqslant 0,\\\ 1&\text{if }s<0.\end{cases}$ (4.16)

The proof of the lemma can be found in Appendix A.3. The obtained filter
function is now characterized by a varying behavior depending on the actual
measurement $y$. This is expressed via the variable $s$ which represents the
coefficients $\langle z,v_{i}\rangle=\sigma_{i}\langle y,u_{i}\rangle$ (see
(3.15), (3.16)). Whenever the coefficient is positive, the reconstruction
scheme behaves like the squared soft TSVD without bias discussed in Section
4.2, i.e., they share the same filter function for those
$x^{\dagger}\in\mathcal{N}(A)$. In all other cases, the reconstruction method
does not change the coefficients of the data, i.e., it behaves like the
identity. Due to the relationship to the squared soft TSVD, we can immediately
specify those $x^{\dagger}$ fulfilling the local approximation property, i.e.,
for the ReLU network we have

$S=\\{x\in\mathcal{N}(\tilde{A})^{\bot}\,|\,\forall
j\in\mathbb{N}\mathrel{\mathop{\mathchar 58\relax}}\langle
v_{j},x\rangle\geqslant 0\\}.$ (4.17)

Thus, the nonlinearity in the network architecture introduces restricted
approximation capabilities as well as convergence guarantees on a limited
subset of $X$ only.

###  4.4 Linear network with soft thresholding activation

At last, we want to analyze a network with soft thresholding activation, e.g.,
known from LISTA [18]. This function is known to promote sparsity since it
shrinks all coefficients and sets those under a certain threshold to zero.
That is why only training data with sufficiently large coefficients matters
for the result of the training and the condition $|\langle
x,v_{j}\rangle|>\frac{\alpha_{j}}{L}$ is crucial in the following lemma.

######  Lemma 4.5.

Let $(v_{j},\sigma_{j}^{2})_{j}$ be the eigenvectors and eigenvalues of $A$,
$\alpha=(\alpha_{j})_{j},\alpha_{j}\geqslant 0$ and $\varphi_{\theta}=Id-
f_{\theta}$ be an iResNet which solves (4.1) with

  * (i)

$f_{\theta}(x)=\phi_{\alpha}(Wx)$, where $\theta=W$, $\Theta=\\{W\in
L(X)\,|\,\exists(w_{j})_{j\in\mathbb{N}}\mathrel{\mathop{\mathchar
58\relax}}W=\sum_{j\in\mathbb{N}}w_{j}\langle\cdot,v_{j}\rangle v_{j}\\}$ and
$\phi_{\alpha}$ is the soft thresholding function w.r.t. the eigenvectors,
i.e., $\phi_{\alpha}(x)=\sum_{j\in\mathbb{N}}\mathrm{sign}(\langle
x,v_{j}\rangle)\max(0,|\langle x,v_{j}\rangle|-\alpha_{j})v_{j}$,

  * (ii)

the training data set is $\\{x^{(i)}\\}_{i\in\\{1,\ldots,N\\}}$, where for any
$j\in\mathbb{N}\mathrel{\mathop{\mathchar 58\relax}}\exists
i\mathrel{\mathop{\mathchar 58\relax}}\ |\langle
x^{(i)},v_{j}\rangle|>\frac{\alpha_{j}}{L}$.

Then, a solution of (4.1) is
$W=\sum_{j\in\mathbb{N}}w_{j}\langle\cdot,v_{j}\rangle v_{j}$,
$w_{j}=\min\\{\frac{\alpha_{j}}{p_{L,j}}+1-\sigma_{j}^{2},L\\}$,
$p_{L,j}=\frac{\sum_{i\in I_{j}(L)}|\langle
x^{(i)},v_{j}\rangle|^{2}}{\sum_{i\in I_{j}(L)}|\langle
x^{(i)},v_{j}\rangle|}$, $I_{j}(L)=\\{i\in\\{1,...,N\\}|\ |\langle
x^{(i)},v_{j}\rangle|>\frac{\alpha_{j}}{L}\\}$ and (3.15) holds with

$r_{L}(\sigma_{j}^{2},s)=\begin{cases}\frac{1}{\max\\{\sigma_{j}^{2}-\frac{\alpha_{j}}{p_{L,j}},1-L\\}}\frac{|s|-\alpha_{j}}{|s|}&\text{if
}|s|>\frac{\alpha_{j}}{w_{j}},\\\ 1&\text{if
}|s|\leqslant\frac{\alpha_{j}}{w_{j}}.\end{cases}$ (4.18)

For singular values $\sigma_{j}^{2}=1$, $w_{j}$ is not uniquely determined.

The proof of the lemma can be found in Appendix A.4. It follows the same line
of reasoning as in the previous sections but is more technical due to the
effects of the nonlinearity.

So far, the filter function $r_{L}$ in Lemma 4.5 is only defined on the
discrete values $\sigma_{j}^{2}$ (and not continuous for $\sigma^{2}\in[0,1]$)
since it depends on the coefficients $p_{L,j}$, $\alpha_{j}$ and $w_{j}$.
However, if we assume continuous extensions $p_{L}(\sigma^{2})$ with
$p_{L}(\sigma_{j}^{2})=p_{L,j}$, $w_{L}(\sigma^{2})$ with
$w_{L}(\sigma_{j}^{2})=w_{j}$, and $\alpha(\sigma^{2})$ with
$\alpha(\sigma_{j}^{2})=\alpha_{j}$, the function $r_{L}=r_{L}(\sigma^{2},s)$
also becomes continuous. The continuity at the point
$|s|=\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}$ is assured by

$\frac{1}{\max\\{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})},1-L\\}}\frac{\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}-\alpha(\sigma^{2})}{\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}}=\frac{1-w_{L}(\sigma^{2})}{\max\\{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})},1-L\\}}=1=r_{L}(\sigma^{2},s).$
(4.19)

To be able to interpret the filter function, suitable values for $s$
representing the coefficients $\langle z,v_{i}\rangle=\sigma_{i}\langle
y,u_{i}\rangle$ (see (3.15), (3.16)) need to be considered. One reasonable
option is to choose $s$ according to the data on which $\varphi_{\theta}$ has
been trained. We thus consider an eigenvector $v_{j}$ scaled with the
coefficient $p_{L,j}$. Since $\varphi_{\theta}$ minimizes (4.1), we expect
$\varphi_{\theta}(p_{L,j}v_{j})\approx
p_{L,j}Av_{j}=p_{L,j}\sigma_{j}^{2}v_{j}$ and
$\varphi_{\theta}^{-1}(p_{L,j}\sigma_{j}^{2}v_{j})\approx p_{L,j}v_{j}$,
respectively. Accordingly, we choose $|s|$ proportional to
$p_{L}(\sigma^{2})\sigma^{2}$, i.e., $|s|=\gamma\,p_{L}(\sigma^{2})\sigma^{2}$
for different values $\gamma>0$. Hence, the case $\gamma=1$ corresponds to a
test vector $z$ with coefficients $|\langle
z,v_{i}\rangle|=|s|=p_{L}(\sigma_{j}^{2})\sigma_{j}^{2}$, which perfectly fits
to the training data. Analogously, the cases $\gamma<1$ (and $\gamma>1$,
respectively) correspond to test data whose coefficients are smaller (or
bigger, respectively) than the average coefficients of the training data.

For $\gamma=1$ the filter function $r_{L}$ can be written in a form, which
allows for an easier interpretation. It holds

$r_{L}\left(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2}\right)=\begin{cases}1&\text{if
}\,\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})},\\\
\frac{1}{1-L}\left(1-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}\right)&\text{if
}\,\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})}<\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L,\\\
\frac{1}{\sigma^{2}}&\text{if
}\,\sigma^{2}>\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L.\end{cases}$
(4.20)

The derivation can be found in Appendix A.5. Note that the filter function
depends especially on the quotient of $\alpha(\sigma^{2})$ (soft thresholding
parameter) and $p_{L}(\sigma^{2})$ (property of training data). To visualize
the influence, we depicted the graph in Figure 5 for two different (constant)
values of $\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}$. As can be seen, the
graph of
$\sigma^{2}\mapsto\sigma^{2}r(\sigma^{2},p_{L}(\sigma^{2})\sigma^{2})$ has two
kinks. The first one depends mainly on the choice of
$\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}$, the second one mainly on the
choice of $L$.

For $\gamma\neq 1$, the filter functions cannot be written in a similar
compact and easy form as in (4.20). Instead, we illustrate them in Figure 6.
In contrast to the case of $\gamma=1$, the graph of
$\sigma^{2}\mapsto\sigma^{2}r(\sigma^{2},p_{L}(\sigma^{2})\sigma^{2}$ is not
equal to one for the larger values of $\sigma^{2}$.

$\sigma^{2}$$\sigma^{2}r_{L}(\sigma,\pm\,p_{L}\sigma^{2})$00.51$0$$0.5$$1$$L=0.8$$L=0.9$$\nicefrac{{\alpha}}{{p_{L}}}=0.1$

$\sigma^{2}$$\sigma^{2}r_{L}(\sigma,\pm\,p_{L}\sigma^{2})$00.51$0$$0.5$$1$$L=0.8$$L=0.9$$\nicefrac{{\alpha}}{{p_{L}}}=0.05$

Figure 5: The plots depict graphs of
$\sigma^{2}r_{L}(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2})$ as given in
(4.18) and (4.20). The choice $s=p_{L}(\sigma^{2})\sigma^{2}$ corresponds to
the case that the test data is similar to the average training data. In this
case, $r_{L}$ shows a quite similar behavior as the squared soft TSVD
function.

$\sigma^{2}$$\sigma^{2}r_{L}(\sigma,\pm\,\gamma
p_{L}\sigma^{2})$00.51$0$$0.5$$1$$L=0.8$$L=0.9$$\alpha=0.1,p_{L}=1$$\gamma=0.85$

$\sigma^{2}$$\sigma^{2}r_{L}(\sigma,\pm\,\gamma
p_{L}\sigma^{2})$00.51$0$$0.5$$1$$L=0.8$$L=0.9$$\alpha=0.1,p_{L}=1$$\gamma=1.15$

Figure 6: The plots depict graphs of $\sigma^{2}r(\sigma^{2},\pm\,\gamma
p_{L}(\sigma^{2})\sigma^{2})$ as given in (4.18). The case $\gamma\neq 1$
corresponds to test data which differs from the average training data.
Especially for larger singular values $\sigma^{2}$ (to the right of the two
kinks), the filter function shows a suboptimal behavior, compared to the
squared soft TSVD filter function.

Finally, we want to analyze to which extend $\varphi_{\theta}$ fulfills the
local approximation property (3.1). The illustration of the filter function in
Figure 5 implies that convergence can only be possible if both kinks tend to
zero. So, we need $L\to 1$ and $\alpha\to 0$. Thus, the structure of the
nonlinearity (soft thresholding) has a severe influence on the regularization
properties of the reconstruction scheme. But in contrast to the ReLU
architecture, the soft thresholding operator provides the opportunity to
control its similarity to a linear operator which can be controlled via the
coefficients $\alpha$. In the following lemma, we discuss in more detail how
$\alpha$ can be chosen depending on the regularization parameter $L$, to
obtain the desired regularization properties.

######  Lemma 4.6.

Let all assumptions of Lemma 4.5 hold. Further, assume that $\tilde{A}$ is
compact and let $p^{\dagger}=\sum_{j\in\mathbb{N}}p_{j}^{\dagger}v_{j}$ be
given by $p_{j}^{\dagger}=\frac{\sum_{i=1}^{N}|\langle
x^{(i)},v_{j}\rangle|^{2}}{\sum_{i=1}^{N}|\langle x^{(i)},v_{j}\rangle|}$. In
addition, consider $x^{\dagger}\in X$ as well as strictly monotonic and
continuous architecture parameter choices
$\alpha_{j},\beta\mathrel{\mathop{\mathchar
58\relax}}(0,1)\rightarrow[0,\infty)$ with

$\alpha_{j}(L)\leqslant p_{j}^{\dagger}\beta(L),\text{ with
}\beta(L)=\mathcal{O}((1-L)\psi(1-L))$ (4.21)

for an index function $\psi\mathrel{\mathop{\mathchar
58\relax}}[0,1]\to[0,\infty)$.

Then, the local approximation property (3.1) holds for any
$x^{\dagger}\in\mathcal{N}(A)^{\perp}$.

###### Proof.

We first verify that $p_{L},p^{\dagger}\in X$
($p_{L}=\sum_{j\in\mathbb{N}}p_{L,j}v_{j}$). For this, we exploit for each
$j\in\mathbb{N}$ that

$\frac{|\langle x^{(i)},v_{j}\rangle|}{\sum_{i=1}^{N}|\langle
x^{(i)},v_{j}\rangle|}\leqslant 1$ (4.22)

for any $i$. We thus obtain via Young’s inequality

$\|p^{\dagger}\|^{2}=\sum_{j\in\mathbb{N}}(p_{j}^{\dagger})^{2}\leqslant\sum_{j\in\mathbb{N}}\bigg{(}\sum_{i=1}^{N}|\langle
x^{(i)},v_{j}\rangle|\bigg{)}^{2}\leqslant
2^{N-1}\sum_{i=1}^{N}\|x^{(i)}\|^{2}<\infty.$ (4.23)

Analogously for $p_{L}$. Due to the finite nature of the data set and the
properties of $\alpha$ we immediately have $p_{L}\rightarrow p^{\dagger}$ for
$L\to 1$ and $p_{L}=p^{\dagger}$ for $L$ being sufficiently close to $1$.

We now continue with the approximation property making use of the notation
$x_{j}\mathrel{\mathop{\mathchar 58\relax}}=\langle x^{\dagger},v_{j}\rangle$
and a sufficiently large $L<1$ where
$w_{j}=\min\\{\frac{\alpha_{j}(L)}{p^{\dagger}_{j}}+1-\sigma_{j}^{2},L\\}\geqslant
0$:

$\displaystyle\|\varphi_{\theta(L)}(x^{\dagger})-Ax^{\dagger}\|^{2}$
$\displaystyle=\sum_{j\in\mathbb{N}}\left(\max(0,w_{j}|x_{j}|-\alpha_{j}(L))-(1-\sigma_{j}^{2})|x_{j}|\right)^{2}$
$\displaystyle\leqslant
2\sum_{j\in\mathbb{N}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}+(w_{j}|x_{j}|-\max(0,w_{j}|x_{j}|-\alpha_{j}(L)))^{2}$
$\displaystyle=2\sum_{j\in\mathbb{N}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}+(\max(0,w_{j}|x_{j}|)-\max(0,w_{j}|x_{j}|-\alpha_{j}(L)))^{2}$
$\displaystyle\leqslant
2\sum_{j\in\mathbb{N}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}+\alpha_{j}(L)^{2}$
$\displaystyle\leqslant
2\sum_{j\in\mathbb{N}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}+2\beta(L)^{2}\|p^{\dagger}\|^{2}$
(4.24)

due to Young’s inequality, the Lipshitz continuity of the ReLU function, and
the assumption on $\alpha_{j}$. We further obtain

$\displaystyle\sum_{j\in\mathbb{N}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}$
$\displaystyle=\underset{=\mathrel{\mathop{\mathchar
58\relax}}(I)}{\underbrace{\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L+\frac{\alpha_{j}(L)}{p_{j}^{\dagger}}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}}}$
$\displaystyle+\underset{=\mathrel{\mathop{\mathchar
58\relax}}(II)}{\underbrace{\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}>1-L+\frac{\alpha_{j}(L)}{p_{j}^{\dagger}}}(1-\sigma_{j}^{2}-w_{j})^{2}|x_{j}|^{2}}}.$
(4.25)

We thus need to derive estimates for $(I)$ and $(II)$.

$\displaystyle(I)$ $\displaystyle=\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L+\frac{\alpha_{j}(L)}{p_{j}^{\dagger}}}(1-\sigma_{j}^{2}-L)^{2}|x_{j}|^{2}$
$\displaystyle\leqslant\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L+\beta(L)}(1-\sigma_{j}^{2}-L)^{2}|x_{j}|^{2}$
$\displaystyle\leqslant\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L}(1-\sigma_{j}^{2}-L)^{2}|x_{j}|^{2}+\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}>1-L\land\sigma_{j}^{2}\leqslant
1-L+\beta(L)}(1-\sigma_{j}^{2}-L)^{2}|x_{j}|^{2}$
$\displaystyle\leqslant\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L}(1-L)^{2}|x_{j}|^{2}+\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}>1-L\land\sigma_{j}^{2}\leqslant
1-L+\beta(L)}\beta(L)^{2}|x_{j}|^{2}$
$\displaystyle\leqslant(1-L)^{2}\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant
1-L}|x_{j}|^{2}+\beta(L)^{2}\|x^{\dagger}\|^{2},$ (4.26)

where we again exploit Young’s inequality and the properties of $a_{j}$. For
$(II)$ we immediately obtain

$(II)=\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}>1-L+\frac{\alpha_{j}(L)}{p_{j}^{\dagger}}}\left(\frac{\alpha_{j}(L)}{p^{\dagger}_{j}}\right)^{2}|x_{j}|^{2}\leqslant\beta(L)^{2}\|x^{\dagger}\|^{2}.$
(4.27)

Due to the properties of $\beta$, $(II)$ already has the desired
characteristics for $L\to 1$. In contrast, $(I)$ requires some further
calculations. Following the same line of reasoning as in the proof of
Corollary 4.2 we can find an index function $\psi^{\prime}$ for any
$\varepsilon>0$ such that

$\displaystyle\sum_{j\mathrel{\mathop{\mathchar
58\relax}}\sigma_{j}^{2}\leqslant 1-L}\langle x^{\dagger},v_{j}\rangle^{2}$
$\displaystyle\leqslant(1+\varepsilon)^{2}\|x^{\dagger}\|^{2}\psi^{\prime}\left(1-L\right)^{2}$
$\displaystyle=\mathcal{O}(\psi^{\prime}(1-L)^{2}).$ (4.28)

Combining this with (4.24), (4.26) and (4.27) we obtain the desired result. ∎

######  Remark 4.3.

Analogous to the line of reasoning in the proof of Lemma 4.3 we split the
series in two sums, (I) and (II). (I) takes care of small singular values and
needs to be related to a property of the element $x^{\dagger}$ in terms of a
certain kind of source condition. The second sum (II) is somehow related to
the approximation properties of the underlying architecture. In the proof of
Lemma 4.3 it is immediately zero for any $x^{\dagger}$ which is due to the
linear network assumption therein. In contrast, in the previous proof we had
to carefully control the behavior of this term since it is strongly influenced
by the structure of the nonlinearity.

In general, the previous lemma serves as an example for the interplay between
approximation capabilities of the network and the resulting regularization
properties.

##  5 Numerical experiments

In this section, we present numerical results in order to compare our
theoretical findings from Section 3 and Section 4 to its corresponding
numerical implementations and extend these by experiments in a more general
setting by learning from noisy measurements. To this end, we use a Radon
transform with $30$ angles and $41$ detector pixels as our forward operator.
We discretize the Radon transform for $28\times 28$ px images. This results in
a linear forward operator $\tilde{A}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{28\times 28}\to\mathbb{R}^{30\times 41}$. For the
training of the different iResNet architectures we use the well-known MNIST
dataset of handwritten digits [25]. This dataset is split into $60\,000$
images for training and $10\,000$ images for testing. For our experiments we
treat the $28\times 28$ images as column vectors of length $784$ and use fully
connected layers in all network architectures. We optimize the networks using
Adam [23] and a learning rate scheduling scheme.

There are multiple ways to satisfy the Lipschitz constraint during training.
Bungert et al. [11] use an additional penalty term in the loss function.
However, this method does not strictly enforce the constraint. Behrmann et al.
[8] use contractive nonlinearities and only constrain the Lipschitz constant
of each individual linear mapping. They compute the Lipschitz constant using a
power iteration and normalize the weights after each training step. We observe
an improvement of convergence when directly parameterizing the weights to
fulfill a specific Lipschitz constant by extending the approach in [30].

In Section 5.1, we show experiments on the local approximation property in
Theorem 3.1. To this end, we use the diagonal architecture proposed in Section
3.1. The resulting data-dependent filter functions are visualized in Section
5.2 and the convergence property is considered in Section 5.3.

We implement each subnetwork $f_{\theta,j},\ j=1,\dots,n$, (see (3.12)) as a
small fully-connected neural network with independent weights for each $j$,
where $n$ denotes the total number of singular values in the discrete setting.
Each subnetwork consists of three layers, where the first two layers
$f_{\theta,j}^{k}$, $k=1,2$ each contain $35$ hidden neurons and its final
layer $f_{\theta,j}^{3}$ contains $1$ neuron. Every layer is equipped with a
linear matrix multiplication and corresponding additive bias and a ReLU
activation function ($k=1,2$). Accordingly, each subnetwork has $1366$
trainable parameters. An illustration of the architecture is provided in
Figure 7. Altogether, the parameters of the subnetworks determine the
parameters $\theta\in\Theta(L)$ of $\varphi_{\theta}=\mathrm{Id}-f_{\theta}$,
where $\Theta(L)$ includes the network parameters as well as the Lipschitz
constraint being realized by constraints on the network parameter. Here, we
enforce the Lipschitz constants $\mathrm{Lip}(f_{\theta,j})\leqslant L$,
$j=1,\ldots,n$, by constraining $\mathrm{Lip}(f_{\theta,j}^{k})\leqslant 1$
for $k=1,2$ and $\mathrm{Lip}(f_{\theta,j}^{3})\leqslant L$. In the following,
we thus write $\theta(L)$ in order to give emphasis to the regularization
parameter. Our training objective is the minimization of the _approximation
loss_ (4.1), i.e., minimize the loss

$l(\varphi_{\theta},A)=\frac{1}{N}\sum_{i}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}$
(5.1)

subject to $\mathrm{Lip}(f_{\theta})\leqslant L$.

$x_{j}$$f_{\theta,j}(x_{j})$
$f_{\theta,j}^{1}$$f_{\theta,j}^{2}$$f_{\theta,j}^{3}$ReLUReLU$35$ neurons$35$
neurons$\mathrm{Lip}(f_{\theta,j}^{1})\leqslant
1$$\mathrm{Lip}(f_{\theta,j}^{2})\leqslant
1$$\mathrm{Lip}(f_{\theta,j}^{3})\leqslant L$ Figure 7: Illustration of the
subnetworks $f_{\theta,j}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\to\mathbb{R},\ j=1,\ldots,n$. Each having $1366$
trainable parameters.

###  5.1 Local approximation property

According to Theorem 3.1, we expect $\varphi_{\theta(L)}^{-1}$ to act as a
regularizer as soon as the local approximation property (3.1) is fulfilled.
Note that the architecture does not change for varying $L$ in our experiments
such that we write $\varphi_{\theta(L)}$ instead of $\varphi_{\theta(L),L}$.
To be able to observe this behavior of our trained network, the local
approximation property needs to be verified with respect to the trained
network parameters. Therefore, we evaluate the approximation error with
respect to fixed data samples

$\mathcal{E}_{\mathrm{x^{(i)}}}(\varphi_{\theta(L)},A)=\|\varphi_{\theta(L)}(x^{(i)})-Ax^{(i)}\|$
(5.2)

as well as the mean approximation error over the test data set

$\mathcal{E}_{\mathrm{mean}}(\varphi_{\theta(L)},A)=\frac{1}{N}\sum_{i}\|\varphi_{\theta(L)}(x^{(i)})-Ax^{(i)}\|.$
(5.3)

Figure 8 indicates superlinear convergence of $\mathcal{E}_{\mathrm{mean}}$
for $L\to 1$, i.e., the existence of an index function $\psi$ such that the
local approximation property (3.1) is satisfied within our test data set on
average. This underpins the capability of the chosen network architecture and
training to locally approximate the operator $A$ in terms of (3.1).
Furthermore, the evaluation of $\mathcal{E}_{x^{(1)}}$ shows that the chosen
sample $x^{(1)}$ fulfills the superlinear convergence and behaves very similar
to the mean error over the test data set. However, from the selection of data
sample $x^{(2)}$ and corresponding error $\mathcal{E}_{x^{(2)}}$ we notice
that the local approximation property does not hold for all samples of our
test data set. Figure 8 additionally shows that some coefficients
$x^{(2)}_{j}$ of $x^{(2)}$, corresponding to large singular values, severely
deviate from the corresponding mean values with respect to the data
distribution. This effect is reduced for $x^{(1)}$. Therefore, the slow,
respectively non-, convergence of $\mathcal{E}_{x^{(2)}}$ could possibly be
explained by the fact that structures outside of the core of the distribution
of chosen singular values have not been learned properly during network
training.

(a) $x^{(1)}$

(b) $x^{(2)}$

$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{-3}$$10^{-2}$$10^{-1}$$1-L$$\mathcal{O}\big{(}(1-L)\big{)}$$\mathcal{O}\big{(}(1-L)^{2}\big{)}$$\mathcal{E}_{\mathrm{mean}}(\varphi_{\theta(L)},A)$$\mathcal{E}_{x^{(1)}}(\varphi_{\theta(L)},A)$$\mathcal{E}_{x^{(2)}}(\varphi_{\theta(L)},A)$$0$$250$$500$$750$$0$$5$$10$$15$$20$$25$$30$$j$$|x^{(2)}_{j}-\frac{1}{N}\sum_{i}x^{(i)}_{j}|^{2}$$|x^{(1)}_{j}-\frac{1}{N}\sum_{i}x^{(i)}_{j}|^{2}$

Figure 8: Test samples $x^{(1)}$ and $x^{(2)}$ (left). Evaluations of
$\mathcal{E}_{\mathrm{mean}}(\varphi_{\theta(L_{m})},A),\
\mathcal{E}_{x^{(1)}}(\varphi_{\theta(L_{m})},A)$ and
$\mathcal{E}_{x^{(2)}}(\varphi_{\theta(L_{m})},A)$ for
$L_{m}=1-\nicefrac{{1}}{{2^{m}}}$ with $m=1,\ldots,8$ on the MNIST test data
set (middle). Squared absolute differences of the coefficients
$x_{j}^{(k)}=\langle x^{(k)},v_{j}\rangle$ for $k=1,2$ to the mean value of
corresponding coefficients in the test data set (right).

###  5.2 Data-dependent filter functions for diagonal architecture

For the experiments in this and the subsequent section we train networks
$\varphi_{\theta(L)}$ with Lipschitz bounds
$L\in\\{L_{m}=1-\nicefrac{{1}}{{3^{m}}}\,|\,m=1,\ldots 5\\}$. We also include
additive Gaussian noise in the network training via minimizing the
approximation loss (4.1). More precisely, for each training sample $x^{(i)}$
we generate $Ax^{(i)}+\eta^{(i)}$ with
$\eta^{(i)}\sim\mathcal{N}(0,\delta_{\ell}Id)$ and relative noise level
${\delta_{\ell}=\hat{\delta}_{\ell}\cdot\mathrm{std}_{\mathrm{MNIST}}}$, where
${\mathrm{std}_{\mathrm{MNIST}}}$ denotes the averaged standard deviation of
the coefficients $\langle x^{(i)},v_{j}\rangle$ of the MNIST data set
(standard deviation with respect to $i$, mean with respect to $j$) and

$\hat{\delta}_{\ell}=\begin{cases}\left(\frac{1}{3}\right)^{\ell}&\text{ for
}\ell>0\\\ 0&\text{ for }\ell=0.\end{cases}$ (5.4)

The loss function $l$ then amounts to

$\min_{\theta\in\Theta(L)}l(\varphi_{\theta},A)=\min_{\theta\in\Theta(L)}\frac{1}{N}\sum_{i=1}^{N}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}-\eta^{(i)}\|^{2}.$
(5.5)

Note that this includes (4.1) in the noise-free case. Trained networks on
noisy samples with noise level $\delta$ with a particular Lipschitz bound $L$
are denoted by $\varphi_{\theta(L,\delta)}$. The noise free training outcome
is denoted by $\varphi_{\theta(L)}$.

Utilizing identical network architectures as in the previous subsection, we
evaluate the learned filter functions of chosen networks. Analogously to
(3.19) the reconstruction can be written as

$T_{L}(y)=\varphi_{\theta}^{-1}(\tilde{A}^{\ast}y)=\hat{b}_{L}+\sum_{j\in\mathbb{N}}\hat{r}_{L}(\sigma_{j}^{2},\sigma_{j}\langle
y,u_{j}\rangle)\sigma_{j}\langle y,u_{j}\rangle v_{j}$ (5.6)

for $\hat{b}_{L}\in X=\mathbb{R}^{n}$. As stated in Section 3.1, the learned
filter function of a trained network with diagonal architecture follows
immediately from its subnetworks. Due to the additional bias, we make a small
adaption to (3.17), which gives

$(\mathrm{Id}-f_{\theta,j})^{-1}(s)-\hat{b}_{L,j}=r_{L}(\sigma_{j}^{2},s)s\qquad\text{for}\
s\in\mathbb{R},$ (5.7)

where each entry $\hat{b}_{L,j}=\langle b_{L},v_{j}\rangle$ corresponds to the
axis intercept of a subnetwork $\varphi_{\theta,j},\ j=1,\ldots,n$, i.e.,
$\hat{b}_{L,j}=(\mathrm{Id}-f_{\theta,j})^{-1}(0)$. Since $f_{\theta,j}$
corresponds to $\sigma_{j}$ or $\sigma_{j}^{2}$, respectively, we also write
$\hat{b}_{L,j}=\hat{b}_{L}(\sigma_{j}^{2})$.

Adapting to the mean values with respect to the distribution of each
coefficient, we compute

$\mu_{j}\mathrel{\mathop{\mathchar 58\relax}}=\frac{1}{N}\sum_{i=1}^{N}\langle
x^{(i)},v_{j}\rangle\quad\text{for}\ j=1,\ldots,n$ (5.8)

and evaluate the filter function with respect to its second variable at the
continuous extension $\sigma^{2}\mu(\sigma^{2})$ with
$\mu(\sigma_{j}^{2})=\mu_{j}$, $j=1,\ldots,n$.

From Figure 9 we notice that the regularization of smaller singular values
increases with decreasing Lipschitz bound $L$ for the mean MNIST data sample.
This is in line with the regularization theory, as $L$ serves as our
regularization parameter. Thus, the filter plots in Figure 9 visualize how $L$
acts as a regularization parameter. The trained diagonal networks
$\varphi_{\theta(L_{m})}$ for $L_{2}$ and $L_{5}$ show a slightly different
behavior at $\sigma_{j}^{2}\approx 0.27$. The distributions of the
coefficients $\langle x,v_{j}\rangle_{j=1,\ldots,n}$ are treated independently
and in this particular singular value we observed a wider distribution of
coefficients $\langle x,v_{j}\rangle$ in the dataset. The inherent structure
within the dataset might be one possible explanation for this outlier. In
general, when neglecting the previously mentioned outlier, for the mean MNIST
sample one can observe a similarity to the squared soft TSVD filter function
(see Section 4.2) to some extent. In additon, the observed decay of
$\|\hat{b}_{L}\|$ with increased $L$ is also in line with the necessary
condition for a classical filter function with bias to become a convergent
regularization (see Lemma 3.3). The observed increase for the largest $L$ is
most likely caused by a numerical instability when evaluating
$(\mathrm{Id}-f_{\theta,j})^{-1}(0)$ via the fixed point iteration performed
for 30 iterations.

$0$$0.25$$0.5$$0.75$$1$$0$$0.5$$1$$\sigma^{2}$$\sigma^{2}\hat{r}_{L}(\sigma^{2},\sigma^{2}\mu(\sigma^{2}))$$L_{5}=0.996$$L_{2}=0.889$$L_{1}=0.667$$0.7$$0.8$$0.9$$0$$0.5$$1$$L$$\|\hat{b}_{L}\|$
Figure 9: Learned filter functions
$r_{L}(\sigma^{2},\sigma^{2}\mu(\sigma^{2}))$ of the iResNet with diagonal
architecture and zero noise level $\delta_{0}$ for different Lipschitz
constraints $L$ (left) where the eigenvalues $\sigma_{j}^{2}$ are highlighted
as black bullets on the x-axis. $\|\hat{b}_{L}\|$ corresponding to the axis
intersection $\hat{b}_{L}=\sum_{j}\hat{b}_{L,j}v_{j}$, evaluated for $L_{m}$,
$m=1,\ldots,5$ (right).

Figure 10 includes reconstructions of a fixed sample from the test data set.
This example illustrates the effect of a regularization induced by the
Lipschitz bound even if the training takes noise into account. It becomes
stronger for small values of $L$, which coincides with our previous
observations. Reconstructions resulting from increased noise levels require
stronger regularizations in order to improve reconstruction quality.
Therefore, the best reconstructions in case of $\hat{\delta}_{0}$ and
$\hat{\delta}_{6}$ result from $\varphi_{\theta(L_{3},\delta_{0})}^{-1}$ and
$\varphi_{\theta(L_{3},\delta_{6})}^{-1}$. In comparison,
$\varphi_{\theta(L_{2},\delta_{3})}^{-1}$ and
$\varphi_{\theta(L_{2},\delta_{4})}^{-1}$ provide improved reconstructions for
the cases $\hat{\delta}_{3}$ and $\hat{\delta}_{4}$. Moreover, at $L_{1}$ we
notice similar blurred structures in the background of the reconstructed digit
for all noise levels. One might argue that its structure compares to a
handwritten digit itself, making the learned background pattern being encoded
in the bias $\hat{b}_{L}$ suitable to the data set. These additional
observations from Figure 10 indicate a dependency of the regularization on
learned data structures. The corresponding filter functions being illustrated
in the right column of Figure 10 show a similar behavior for all training
noise levels which underpins that the outcome of the approximation training is
independent of the noise for a sufficiently large number of training samples.
The outlier in the filter function for the mean sample $(\mu_{j})_{j}$ can
also be observed in Figure 10. In addition, this test sample has a slightly
different behavior with respect to the second singular value. Note that the
seemingly linear behavior for $\sigma_{j}^{2}>0.4$ is only due to the fact
that this is the linear continuous extension between the first and second
singular value. In summary, the resulting filter behaves similar to the
squared soft TSVD filter independent of the training noise and with exceptions
at two singular values.

Figure 10: Reconstructions of an MNIST sample $\tilde{x}$ from the test data
set by computing
$\varphi_{\theta(L_{m},\delta_{\ell})}^{-1}(A\tilde{x}+\tilde{\eta})$ with
$\tilde{\eta}\sim\mathcal{N}(0,\delta_{\ell}Id)$ for Lipschitz bounds $L_{m}$
(columns) and noise levels
$\delta_{\ell}=\hat{\delta}_{\ell}\cdot\mathrm{std}_{\mathrm{MNIST}}$ with
$\ell=0,6,4,3$ (rows). The last column depicts the filter functions at
$L_{1}$, $L_{2}$, and $L_{3}$ for each noise level with respect to the sample
$\tilde{x}$.

###  5.3 Convergence property

After investigating approximation capabilities of trained networks with
respect to the operator $A$ in Section 5.1 in the noise-free case and
extending the training to the noisy case in Section 5.2, we continue verifying
the convergence property with respect to different noise levels. We analyze
the convergence property by computing the averaged _reconstruction error_

$\mathrm{MSE}_{\mathrm{reco}}^{\delta_{\ell}}(\varphi_{\theta(L,\delta)},A)=\frac{1}{N}\sum_{i=1}^{N}\|x^{(i)}-\varphi_{\theta(L,\delta)}^{-1}(Ax^{(i)}+\eta^{(i)})\|^{2},$
(5.9)

including the noise levels $\delta_{\ell}$ and noise samples
$\eta^{(i)}\sim\mathcal{N}(0,\delta_{\ell}Id)$ in the reconstruction process
where the network has been trained on noise level $\delta$. We thus can
evaluate the noise-free training case, which solely aims to impart data-
dependency, on reconstructions from noisy data and the noisy training case
where training and test noise share the same level. Reconstructions, i.e., the
inverse of the iResNet, are obtained by using 30 fixed point iterations.

Figure 11 shows that $\varphi_{\theta(L,\delta_{0})}^{-1}$ as well as
$\varphi_{\theta(L,\delta_{\ell})}^{-1}$ provides more accurate
reconstructions with respect to (5.9) at large $L$ and low noise levels,
whereas this behavior is reversed for decreasing Lipschitz bound and
increasing noise levels. This is consistent with the regularization theory and
the visualized reconstructions in Figure 10, as high noise levels require
strong regularizations and vice versa. The behavior of
$\mathrm{MSE}_{\mathrm{reco}}^{\delta_{\ell}}(\varphi_{\theta(L,\delta_{\ell})},A)$
for small $L$ and small noise levels $\delta_{\ell}$ is rather intuitive,
since its approximation capability is limited as a consequence of strong
regularization. In addition, from Figure 11 one can also extract a suitable
candidate for the parameter choice $L(\delta)$ to obtain convergence. The
similarity in the behavior of $\varphi_{\theta(L,\delta_{0})}^{-1}$ and
$\varphi_{\theta(L,\delta_{\ell})}^{-1}$ underpins that the outcome of the
approximation training is independent of noise if data and noise are
independent.

$10^{-2}$$10^{-1}$$10^{-2}$$10^{-1}$$\hat{\delta}_{\ell}$$\mathrm{MSE}_{\mathrm{reco}}^{\delta\ell}(\varphi_{\theta(L,\delta_{0})},A)$$10^{-2}$$10^{-1}$$10^{-2}$$10^{-1}$$\hat{\delta}_{\ell}$$\mathrm{MSE}_{\mathrm{reco}}^{\delta\ell}(\varphi_{\theta(L,\delta_{\ell})},A)$$L_{1}=0.667$$L_{2}=0.889$$L_{3}=0.963$$L_{4}=0.988$$L_{5}=0.996$
Figure 11: Reconstruction errors
$\mathrm{MSE}_{\mathrm{reco}}^{\delta_{\ell}}(\varphi_{\theta(L,\delta_{0})},A)$
with training on noise free samples and reconstruction from noisy samples
(left) and
$\mathrm{MSE}_{\mathrm{reco}}^{\delta_{\ell}}(\varphi_{\theta(L,\delta_{\ell})},A)$
with training and reconstruction on the same noise level $\delta_{\ell}$
(right) for $\varphi_{\theta(L_{m},\delta_{\ell})}$ over different noise
levels $\delta_{\ell}=\hat{\delta}_{\ell}\cdot\mathrm{std}_{\mathrm{MNIST}}$,
$\ell=0,\ldots,6$, and Lipschitz bounds $L_{m}$, $m=1,\ldots,5$.

##  6 Discussion and Outlook

In the present work, we developed and investigated the regularization theory
for the proposed iResNet reconstruction approach providing a learned method
from data samples. The network’s local approximation property is fundamental
to delivering a convergent regularization scheme. It comprises approximation
properties of the architecture and training, a definition of solution type,
and a source condition. The approximation loss used for training is motivated
by this property. In the most general version, the framework can be understood
as a fully-learned end-to-end reconstruction scheme with minor limitations as
it relies on the concatenation with $\tilde{A}^{\ast}$, i.e., some a priori
knowledge on the forward operator is at least included in a hand-crafted way
in the reconstruction scheme. Introducing a diagonal architecture type relying
on the SVD of the forward operator $\tilde{A}$ allowed for an analytical
investigation of the resulting learned nonlinear spectral reconstruction
method, which becomes a convergent regularization when fulfilling the local
approximation property. The analysis of trained shallow architectures revealed
the link between the classical linear filter-based regularization theory and
the concept of the local approximation property, and it illustrated the
interplay between the approximation capability of the nonlinear network and
the source condition. In addition, we validated and extended the theoretical
findings by a series of numerical experiments on the MNIST data set and
provided further insights into the learned nonlinear spectral regularization,
such as similarity to the analytically determined linear regularization
(squared soft TSVD) and data-dependency of the learned regularization.

Our iResNet method using the diagonal architecture can be seen as a
generalization of the learned linear spectral regularization considered in
[22]. Using a different loss, which measures the reconstruction error, the
authors of [22] obtain learned filter functions corresponding to Tikhonov
regularization with data- and noise-dependent regularization parameters for
each singular vector. An extension of the iResNet approach to this kind of
training loss is thus desirable for future comparison.

The present work also serves as a starting point for various future research
directions:

As we obtain the inverse of the residual network via a fixed point iteration,
findings in the context of the iResNet reconstruction may be related to
learned optimization methods [6], maximally monotone operators [36], or
potentially to plug-and-play methods [16]. Investigations at the interface to
these methods are far beyond the present work and will remain future work.

The concept of a local criterion to guarantee convergence of the iResNet
reconstruction method provides further opportunities to investigate
regularization properties of learned reconstruction schemes. In a different
context, the authors of [5] also consider a local property. They consider a
learned projection onto subspaces concatenated with the forward operator,
i.e., they consider a projected linear problem and combine it with Tikhonov
regularization. There, the authors guarantee convergence of the regularization
method by the Tikhonov theory, including a different but also localized
assumption. The findings in the present work and the literature thus reveal
the necessity to explicitly take the approximation capabilities of the learned
schemes and the training procedure into account. Potential future
investigations in this context are twofold. On the one hand, one may exploit
universal approximation properties of the underlying architectures to
guarantee convergence for a sufficiently large subset of $X$, e.g., by
verifying the local approximation property. On the other hand, one may adapt
the definition of a regularization to the approximate and data-based nature of
the learning-based setting.

Besides the more general open research questions at the interface to other
methods and data-based concepts, the iResNet reconstruction approach itself
provides various potential future works. The observed similarity of the
learned nonlinear filter functions in the numerical experiments to the
analytical solution of the affine network with bias (Section 4.2) immediately
raises the question: how much does the data’s structure influence the
resulting reconstruction scheme? Phrased differently, one needs to test the
limits of the proposed approximation training. Here, it would be desirable to
change the perspective to a Bayesian one considering data and noise
distributions and further investigate the training outcome, including
additional variations of the training loss and its influence on the local
approximation property. Natural extension and generalizations, such as
including learned solution types, e.g., by combining our approach with null
space networks [40] or relaxing the diagonal assumption of the architecture,
remain future work.

In summary, the present first work on learned iResNet regularization schemes
builds the basis for various future works of theoretical and numerical nature.

## Acknowledgments

A. Denker, N. Heilenkötter, and M. Iske acknowledge the support by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project
number 281474342/GRK2224/2. T. Kluth acknowledges support from the DELETO
project funded by the Federal Ministry of Education and Research (BMBF,
project number 05M20LBB).

##  Appendix A Proofs of Section 4

###  A.1 Proof of Lemma 4.1

Equivalent filter function to the one-parameter-network.

###### Proof.

At first, we observe that the Lipschitz constraint in (4.1) is fulfilled for
$|k|\leqslant L$. Thus, (4.1) is equivalent to

$\displaystyle\min_{|k|\leqslant
L}\sum_{i}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{|k|\leqslant
L}\sum_{i}\|x^{(i)}-kx^{(i)}+kAx^{(i)}-Ax^{(i)}\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{|k|\leqslant
L}\sum_{i}(1-k)^{2}\|x^{(i)}-Ax^{(i)}\|^{2}.$ (A.1)

Since there is one $x^{(i)}$ s.t. $Ax_{i}\neq x^{(i)}$ the solution is
obviously $k=L$ and we have $\varphi_{\theta}(x)=x-L(x-Ax)$.

Now we use the SVD of $\tilde{A}$ to solve $\varphi_{\theta}(\bar{x})=z$ for
$z\in\mathcal{R}(A)$ component wise. For all $j\in\mathbb{N}$ it holds

$\displaystyle\langle\bar{x}-L(\bar{x}-A\bar{x}),v_{j}\rangle$
$\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle(1-L)\langle\bar{x},v_{j}\rangle+L\langle\bar{x},Av_{j}\rangle$
$\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle(1-L+L\sigma_{j}^{2})\langle\bar{x},v_{j}\rangle$
$\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle\langle\bar{x},v_{j}\rangle$
$\displaystyle=\frac{1}{1-L+L\sigma_{j}^{2}}\langle z,v_{j}\rangle$ (A.2)

Thus, the filter function

$r_{L}(\sigma^{2},s)=\frac{1}{1-L+L\sigma^{2}}=\frac{1}{L}\cdot\frac{1}{\alpha+\sigma^{2}},\qquad\alpha=\frac{1-L}{L}$
(A.3)

fulfills (3.15). ∎

###  A.2 Proof of Lemma 4.2

Equivalent filter function to the linear network.

###### Proof.

At first, we observe that the Lipschitz constraint in (4.1) is fulfilled if
the eigenvalues $(w_{j})_{j}$ of $W$ are restricted by $|w_{j}|\leqslant L$.
Defining $\Theta_{L}=\left\\{(W,b)\in\Theta\,\big{|}\,|w_{j}|\leqslant
L\right\\}$, (4.1) is equivalent to

$\displaystyle\min_{\theta\in\Theta_{L}}\,\sum_{i=1}^{N}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}$
$\displaystyle\Leftrightarrow$
$\displaystyle\min_{(W,b)\in\Theta_{L}}\,\sum_{i=1}^{N}\left\|x^{(i)}-Wx_{i}-b-Ax^{(i)}\right\|^{2}$
$\displaystyle\Leftrightarrow$
$\displaystyle\min_{(W,b)\in\Theta_{L}}\,\sum_{i=1}^{N}\|P_{\mathcal{N}(A)}(x^{(i)})\|^{2}+\left\|P_{\mathcal{N}(A)^{\bot}}(x^{(i)})-Wx_{i}-b-Ax^{(i)}\right\|^{2}.$
(A.4)

The part with $P_{\mathcal{N}(A)}(x^{(i)})$ is independent of $W$ and $b$ and
therefore not relevant. Furthermore, it holds

$\displaystyle P_{\mathcal{N}(A)^{\bot}}(x^{(i)})-Wx_{i}-b-Ax^{(i)}$
$\displaystyle=\sum_{j\in\mathbb{N}}\langle x,v_{j}\rangle
v_{j}-\sum_{j\in\mathbb{N}}\left(w_{j}\langle
x,v_{j}\rangle+b_{j}\right)v_{j}-\sum_{j\in\mathbb{N}}\sigma_{j}^{2}\langle
x,v_{j}\rangle v_{j}$
$\displaystyle=\sum_{j\in\mathbb{N}}\left((1-w_{j}-\sigma_{j}^{2})\langle
x,v_{j}\rangle-b_{j}\right)v_{j}.$ (A.5)

Thus, the minimizing problem can be written as

$\min_{|w_{j}|\leqslant
L,\,(b_{j})\in\ell^{2}}\,\sum_{i=1}^{N}\sum_{j\in\mathbb{N}}\left((1-w_{j}-\sigma_{j}^{2})\langle
x,v_{j}\rangle-b_{j}\right)^{2}.$ (A.6)

We can solve this problem for each $j$ separately. For determining $b_{j}$ we
set the derivative to zero and get

$\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle!}}{{=}}-2\sum_{i=^{1}}^{N}\left((1-w_{j}-\sigma_{j}^{2})\langle
x^{(i)},v_{j}\rangle-b_{j}\right)$ $\displaystyle\Leftrightarrow$
$\displaystyle b_{j}$
$\displaystyle\stackrel{{\scriptstyle!}}{{=}}\frac{1}{N}(1-w_{j}-\sigma_{j}^{2})\sum_{i=1}^{N}\langle
x^{(i)},v_{j}\rangle=(1-w_{j}-\sigma_{j}^{2})\mu_{j}.$ (A.7)

Since for every $i$, the sequence $(\langle
x^{(i)},v_{j}\rangle)_{j\in\mathbb{N}}$ is in $\ell^{2}$, $(b_{j})\in\ell^{2}$
is also fulfilled. It remains to solve

$\min_{|w_{j}|\leqslant
L,}\,\sum_{i=1}^{N}\left(\left(1-w_{j}-\sigma_{j}^{2}\right)\left(\langle
x,v_{j}\rangle-\mu_{j}\right)\right)^{2}.$ (A.8)

By assumption, for every $v_{j}$ there is at least one $x^{(i)}$ s.t. $\langle
x^{(i)},v_{j}\rangle\neq\mu_{j}$, thus, the problem can be simplified to

$\min_{|w_{j}|\leqslant L}\,(1-w_{j}-\sigma_{j}^{2})^{2}\qquad\forall
j\in\mathbb{N}.$ (A.9)

The obvious solution is $w_{j}=\min\\{1-\sigma_{j}^{2},L\\}$, since
$1-\sigma_{j}^{2}$ is by assumption always positive.

Now we use the SVD of $\tilde{A}$ to solve $\varphi_{\theta}(\bar{x})=z$ for
some $z\in\mathcal{R}(A)$ component-wise. For all $j\in\mathbb{N}$ it holds

$\displaystyle\langle\bar{x},v_{j}\rangle-w_{j}\langle\bar{x},v_{j}\rangle-
b_{j}$ $\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle(1-w_{j})\langle\bar{x},v_{j}\rangle$ $\displaystyle=\langle
z,v_{j}\rangle+b_{j}$ $\displaystyle\Leftrightarrow$
$\displaystyle\langle\bar{x},v_{j}\rangle$
$\displaystyle=\frac{1}{1-w_{j}}\langle z,v_{j}\rangle+\frac{b_{j}}{1-w_{j}}$
$\displaystyle\Leftrightarrow$ $\displaystyle\langle\bar{x},v_{j}\rangle$
$\displaystyle=\frac{1}{\max\\{\sigma_{j}^{2},1-L\\}}\langle
z,v_{j}\rangle+\frac{\max\\{0,1-L-\sigma_{j}^{2}\\}}{\max\\{\sigma_{j}^{2},1-L\\}}\mu_{j}.$
(A.10)

Thus, using the filter function and bias

$\hat{r}_{L}(\sigma^{2})=\frac{1}{\max\\{\sigma^{2},1-L\\}},\qquad\hat{b}_{L}=\sum_{\sigma_{j}^{2}<1-L}\frac{1-L-\sigma^{2}}{1-L}\mu_{j}v_{j},$
(A.11)

the reconstruction schemes $\varphi_{\theta}^{-1}\circ\tilde{A^{*}}$ and
$T_{L}$ from (3.19) are equivalent. ∎

###  A.3 Proof of Lemma 4.4

Equivalent filter function to the ReLU-network.

###### Proof.

At first, we observe that the Lipschitz constraint in (4.1) is fulfilled if
the eigenvalues $(w_{j})_{j}$ of $W$ are restricted by $|w_{j}|\leqslant L$.
Thus, (4.1) is equivalent to

$\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\left\|x^{(i)}-\phi(Wx_{i})-Ax^{(i)}\right\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\left\|P_{\mathcal{N}(A)}(x^{(i)})+\sum_{j\in\mathbb{N}}\langle
x^{(i)},v_{j}\rangle v_{j}-\sum_{j\in\mathbb{N}}\max\\{0,w_{j}\langle
x^{(i)},v_{j}\rangle\\}v_{j}-\sum_{j\in\mathbb{N}}\sigma_{j}^{2}\langle
x^{(i)},v_{j}\rangle v_{j}\right\|^{2}$ $\displaystyle\Leftrightarrow$
$\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\left\|P_{\mathcal{N}(A)}(x^{(i)})\right\|^{2}+\left\|\sum_{j\in\mathbb{N}}\left((1-\sigma_{j}^{2})\langle
x^{(i)},v_{j}\rangle-\max\\{0,w_{j}\langle
x^{(i)},v_{j}\rangle\\}\right)v_{j}\right\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\sum_{j\in\mathbb{N}}\left((1-\sigma_{j}^{2})\langle
x^{(i)},v_{j}\rangle-\max\\{0,w_{j}\langle
x^{(i)},v_{j}\rangle\\}\right)^{2}.$ (A.12)

At first we focus on the components $\langle x^{(i)},v_{j}\rangle$ which are
negative or zero. They do not influence the solution of the problem since

$\min_{|w_{j}|\leqslant
L}\,\left(-(1-\sigma_{j}^{2})-\max\\{0,-w_{j}\\}\right)^{2}\qquad\forall
j\in\mathbb{N}$ (A.13)

is solved by any $w_{j}\in[0,L]$, as $1-\sigma_{j}^{2}$ is by assumption
always positive. Thus, it suffices to consider the cases $\langle
x^{(i)},v_{j}\rangle>0$, which lead to

$\min_{|w_{j}|\leqslant
L}\,(1-\sigma_{j}^{2}-\max\\{0,w_{j}\\})^{2}\qquad\forall j\in\mathbb{N}.$
(A.14)

The obvious solution is $w_{j}=\min\\{1-\sigma_{j}^{2},L\\}$.

Now we use the SVD of $\tilde{A}$ to solve $\varphi_{\theta}(\bar{x})=z$ for
$z\in\mathcal{R}(A)$ or

$\bar{x}-\sum_{j\in\mathbb{N}}\max\\{0,w_{j}\langle\bar{x},v_{j}\rangle\\}v_{j}=\sum_{j\in\mathbb{N}}\langle
z,v_{j}\rangle v_{j},$ (A.15)

respectively, component-wisely. For all $j\in\mathbb{N}$ it holds

$\displaystyle\langle\bar{x},v_{j}\rangle-\max\\{0,w_{j}\langle\bar{x},v_{j}\rangle\\}$
$\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle\min\\{\langle\bar{x},v_{j}\rangle,(1-w_{j})\langle\bar{x},v_{j}\rangle\\}$
$\displaystyle=\langle z,v_{j}\rangle.$ (A.16)

Since $1-w_{j}>0$, the signs of $\langle\bar{x},v_{j}\rangle$ and $\langle
z,v_{j}\rangle$ must coincide. In case of positive signs, the minimum on the
left hand side takes the value $(1-w_{j})\langle\bar{x},v_{j}\rangle$, in case
of negative signs, it is $\langle\bar{x},v_{j}\rangle$. Thus, we get

$\langle\bar{x},v_{j}\rangle=\begin{cases}\frac{1}{1-w_{j}}\langle
z,v_{j}\rangle&\text{if }\langle z,v_{j}\rangle\geqslant 0,\\\ \langle
z,v_{j}\rangle&\text{if }\langle z,v_{j}\rangle<0.\end{cases}$ (A.17)

Using $1-w_{j}=\max\\{\sigma_{j}^{2},1-L\\}$, it can be seen that

$r_{L}(\sigma^{2},s)=\begin{cases}\frac{1}{\max\\{\sigma^{2},1-L\\}}&\text{if
}s\geqslant 0,\\\ 1&\text{if }s<0\end{cases}$ (A.18)

fulfills (3.15). ∎

###  A.4 Proof of Lemma 4.5

Equivalent filter function to the soft-thresholding-network.

###### Proof.

At first, we observe that the Lipschitz constraint in (4.1) is fulfilled if
the eigenvalues $(w_{j})_{j}$ of $W$ are restricted by $|w_{j}|\leqslant L$.
Besides, $f_{\theta}(x)$ can be written as

$\displaystyle\phi_{\alpha}(Wx)$
$\displaystyle=\sum_{j\in\mathbb{N}}\mathrm{sign}(\langle
Wx,v_{j}\rangle)\max(0,|\langle Wx,v_{j}\rangle|-\alpha_{j})v_{j}$
$\displaystyle=\sum_{j\in\mathbb{N}}\mathrm{sign}(w_{j}\langle
x,v_{j}\rangle)\max(0,|w_{j}\langle x,v_{j}\rangle|-\alpha_{j})v_{j}.$ (A.19)

For any $x^{(i)}$ from the training data set, it follows

$\displaystyle\quad\,\,x^{(i)}-\phi_{\alpha}(Wx^{(i)})-Ax^{(i)}$
$\displaystyle=P_{\mathcal{N}(A)}(x^{(i)})+\sum_{j\in\mathbb{N}}\langle
x^{(i)},v_{j}\rangle v_{j}-\sum_{j\in\mathbb{N}}\mathrm{sign}(w_{j}\langle
x^{(i)},v_{j}\rangle)\max(0,|w_{j}\langle
x^{(i)},v_{j}\rangle|-\alpha_{j})v_{j}-\sum_{j\in\mathbb{N}}\sigma_{j}^{2}\langle
x^{(i)},v_{j}\rangle v_{j}$
$\displaystyle=P_{\mathcal{N}(A)}(x^{(i)})+\sum_{j\in\mathbb{N}}\left[x_{j}^{(i)}-\mathrm{sign}(w_{j}x_{j}^{(i)})\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})-\sigma_{j}^{2}x_{j}^{(i)}\right]v_{j}$
(A.20)

where $x_{j}^{(i)}\mathrel{\mathop{\mathchar 58\relax}}=\langle
x^{(i)},v_{j}\rangle$. Therefore, (4.1) is equivalent to

$\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\|\varphi_{\theta}(x^{(i)})-Ax^{(i)}\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{i}\left\|x^{(i)}-\phi_{\alpha}(Wx^{(i)})-Ax^{(i)}\right\|^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{(w_{j}),\,|w_{j}|\leqslant
L}\,\sum_{j}\sum_{i}\left(x_{j}^{(i)}-\mathrm{sign}(w_{j}x_{j}^{(i)})\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})-\sigma_{j}^{2}x_{j}^{(i)}\right)^{2}.$
(A.21)

For each $j$, we have to find the solution of

$\displaystyle\min_{|w_{j}|\leqslant
L}\sum_{i}\left(x_{j}^{(i)}-\mathrm{sign}(w_{j}x_{j}^{(i)})\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})-\sigma_{j}^{2}x_{j}^{(i)}\right)^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{|w_{j}|\leqslant
L}\sum_{i}\left((1-\sigma_{j}^{2})\mathrm{sign}(x_{j}^{(i)})|x_{j}^{(i)}|-\mathrm{sign}(w_{j})\mathrm{sign}(x_{j}^{(i)})\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})\right)^{2}$
$\displaystyle\Leftrightarrow$ $\displaystyle\min_{|w_{j}|\leqslant
L}\sum_{i}\left((1-\sigma_{j}^{2})|x_{j}^{(i)}|-\mathrm{sign}(w_{j})\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})\right)^{2}.$
(A.22)

Note that for $\sigma_{j}^{2}=1$, any
$w_{j}\in\left[-\frac{\alpha_{j}}{\max^{(i)}|x_{j}^{(i)}|},\frac{\alpha_{j}}{\max_{i}|x_{j}^{(i)}|}\right]$
solves the problem and hence the solution is not unique. In all other cases,
we have $1-\sigma_{j}^{2}>0$, thus, only a positive sign of $w_{j}$ makes
sense and we restrict ourselves to this case. Due to $w_{j}\geqslant 0$, we
have

$\quad\,\,\left((1-\sigma_{j}^{2})|x_{j}^{(i)}|-\max(0,|w_{j}x_{j}^{(i)}|-\alpha_{j})\right)^{2}=\begin{cases}\left((1-\sigma^{2})x_{j}^{(i)}\right)^{2}&\text{if
}|w_{j}x_{j}^{(i)}|\leqslant\alpha_{j},\\\
\left((1-\sigma_{j}^{2}-w_{j})|x_{j}^{(i)}|+\alpha_{j}\right)^{2}&\text{if
}|w_{j}x_{j}^{(i)}|>\alpha_{j}.\end{cases}$ (A.23)

As the first case is an upper bound of the second case, it is highly desirable
to choose $w_{j}$ large enough such that all training samples are contained in
the second case. But we also need to take into account the upper bound
$w_{j}\leqslant L$. Due to the data set assumption (ii), the minimization
problem (A.22) becomes

$\min_{w_{j}\in[0,L]}\sum_{i\in
I_{j}(L)}\left((1-\sigma_{j}^{2}-w_{j})|x_{j}^{(i)}|+\alpha_{j}\right)^{2}.$
(A.24)

The minimizer is

$w_{j}=\begin{cases}\frac{\alpha_{j}}{p_{j}}+1-\sigma_{j}^{2}&\text{if
}\frac{\alpha_{j}}{p_{j}}+1-\sigma_{j}^{2}\leqslant L,\\\
L&\text{else},\end{cases}$ (A.25)

where $p_{L,j}=p_{j}=\frac{\sum_{i\in I_{j}(L)}|x_{j}^{(i)}|^{2}}{\sum_{i\in
I_{j}(L)}|x_{j}^{(i)}|}$. Equivalently, we can write
$w_{j}=\min\\{\frac{\alpha_{j}}{p_{j}}+1-\sigma_{j}^{2},L\\}$.

Now we use the SVD of $\tilde{A}$ to solve $\varphi_{\theta}(\bar{x})=z$ for
$z\in\mathcal{R}(A)$ or

$\bar{x}-\sum_{j\in\mathbb{N}}\mathrm{sign}(w_{j}\langle\bar{x},v_{j}\rangle)\max(0,|w_{j}\langle\bar{x},v_{j}\rangle|-\alpha_{j})v_{j}=\sum_{j\in\mathbb{N}}\langle
z,v_{j}\rangle v_{j},$ (A.26)

respectively, componentwise. For all $j\in\mathbb{N}$ it holds

$\displaystyle\langle\bar{x},v_{j}\rangle-\mathrm{sign}(w_{j}\langle\bar{x},v_{j}\rangle)\max(0,|w_{j}\langle\bar{x},v_{j}\rangle|-\alpha_{j})$
$\displaystyle=\langle z,v_{j}\rangle$ $\displaystyle\Leftrightarrow$
$\displaystyle\mathrm{sign}(\langle\bar{x},v_{j}\rangle)\left(|\langle\bar{x},v_{j}\rangle|-\max(0,w_{j}|\langle\bar{x},v_{j}\rangle|-\alpha_{j})\right)$
$\displaystyle=\langle z,v_{j}\rangle.$ (A.27)

Thus, the sign of $\langle\bar{x},v_{j}\rangle$ is the same as the one of
$\langle z,v_{j}\rangle$. Assuming
$w_{j}|\langle\bar{x},v_{j}\rangle|\leqslant\alpha_{j}$, we obtain

$|\langle\bar{x},v_{j}\rangle|=|\langle z,v_{j}\rangle|\qquad\text{if
}|\langle z,v_{j}\rangle|\leqslant\frac{\alpha_{j}}{w_{j}}$ (A.28)

and assuming $w_{j}|\langle\bar{x},v_{j}\rangle|>\alpha_{j}$, we obtain

$|\langle\bar{x},v_{j}\rangle|=\frac{|\langle
z,v_{j}\rangle|-\alpha_{j}}{1-w_{j}}\qquad\text{if }\frac{|\langle
z,v_{j}\rangle|}{1-w_{j}}-\frac{\alpha_{j}}{1-w_{j}}>\frac{\alpha_{j}}{w_{j}}.$
(A.29)

In total, this leads to

$\langle\bar{x},v_{j}\rangle=\begin{cases}\langle z,v_{j}\rangle&\text{if
}|\langle z,v_{j}\rangle|\leqslant\frac{\alpha_{j}}{w_{j}},\\\
\mathrm{sign}(\langle z,v_{j}\rangle)\frac{|\langle
z,v_{j}\rangle|-\alpha_{j}}{1-w_{j}}&\text{if }|\langle
z,v_{j}\rangle|>\frac{\alpha_{j}}{w_{j}}.\end{cases}$ (A.30)

Using $1-w_{j}=\max\\{\sigma_{j}^{2}-\frac{\alpha_{j}}{p_{j}},1-L\\}$, we see
that the filter function

$r_{L}(\sigma_{j}^{2},s)=\begin{cases}1&\text{if
}|s|\leqslant\frac{\alpha_{j}}{w_{j}},\\\
\frac{1}{\max\\{\sigma_{j}^{2}-\frac{\alpha_{j}}{p_{j}},1-L\\}}\frac{|s|-\alpha_{j}}{|s|}&\text{if
}|s|>\frac{\alpha_{j}}{w_{j}}\end{cases}$ (A.31)

fulfills (3.15). $F_{L}$ is then derived by exploiting (3.18). ∎

###  A.5 Computation of filter function (4.20)

Assuming continuous extensions $p_{L}(\sigma^{2})$ s.t.
$p_{L}(\sigma_{j}^{2})=p_{L,j}$, $w_{L}(\sigma^{2})$ s.t.
$w_{L}(\sigma_{j}^{2})=w_{j}$, and $\alpha(\sigma^{2})$ s.t.
$\alpha(\sigma_{j}^{2})=\alpha_{j}$, we can write the filter function from
Lemma 4.5 as

$r_{L}(\sigma^{2},s)=\begin{cases}\frac{1}{\max\\{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})},1-L\\}}\frac{|s|-\alpha(\sigma^{2})}{|s|}&\text{if
}|s|>\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})},\\\ 1&\text{if
}|s|\leqslant\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}.\end{cases}$ (A.32)

Besides, it holds
$w_{L}(\sigma^{2})=\min\left\\{\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-\sigma^{2},L\right\\}\leqslant
L$, also by Lemma 4.5. Now, we choose $|s|=p_{L}(\sigma^{2})\sigma^{2}$ and
distinguish between three different cases (small, medium and large values of
$\sigma^{2}$).

At first, we consider the case
$\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})}$. Then, it
holds

$|s|=p_{L}(\sigma^{2})\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{L}\leqslant\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}.$
(A.33)

Thus, it follows

$r_{L}\left(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2}\right)=1.$ (A.34)

Secondly, we consider the case
$\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})}<\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L$.
We shortly make sure, that the lower bound is in fact smaller than the upper
one. The assumptions of Lemma 4.5 imply that
$p_{L}(\sigma^{2})>\frac{\alpha(\sigma^{2})}{L}$ holds. From this, we can
follow

$\displaystyle p_{L}(\sigma^{2})(1-L)$
$\displaystyle>\frac{\alpha(\sigma^{2})}{L}-\frac{L\alpha(\sigma^{2})}{L}$
$\displaystyle\Rightarrow$
$\displaystyle\alpha(\sigma^{2})+p_{L}(\sigma^{2})(1-L)$
$\displaystyle>\frac{\alpha(\sigma^{2})}{L}$ $\displaystyle\Rightarrow$
$\displaystyle\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L$
$\displaystyle>\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})}.$ (A.35)

Besides, note that
$\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L$ implies
$w_{L}(\sigma^{2})=L$. Thus, it holds

$|s|=p_{L}(\sigma^{2})\sigma^{2}>\frac{\alpha(\sigma^{2})}{L}=\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}.$
(A.36)

Accordingly, we obtain

$\displaystyle r_{L}\left(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2}\right)$
$\displaystyle=\frac{1}{\max\\{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})},1-L\\}}\frac{p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}$
$\displaystyle=\frac{1}{1-L}\frac{p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}=\frac{1}{1-L}\left(1-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}\right).$
(A.37)

At last, we consider the case
$\sigma^{2}>\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L$. Note that this
implies
$w_{L}(\sigma^{2})=\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-\sigma^{2}$.
Besides it holds $\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}>0$,
which implies

$\displaystyle\left(\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}\right)(1-\sigma^{2})$
$\displaystyle>0$ $\displaystyle\Rightarrow$
$\displaystyle\sigma^{2}\left(\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-\sigma^{2}\right)$
$\displaystyle>\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}$
$\displaystyle\Rightarrow$ $\displaystyle|s|=p_{L}(\sigma^{2})\sigma^{2}$
$\displaystyle>\frac{\alpha(\sigma^{2})}{\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-\sigma^{2}}=\frac{\alpha(\sigma^{2})}{w_{L}(\sigma^{2})}.$
(A.38)

Accordingly, we obtain

$\displaystyle r_{L}\left(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2}\right)$
$\displaystyle=\frac{1}{\max\\{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})},1-L\\}}\frac{p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}$
$\displaystyle=\frac{1}{\sigma^{2}-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}}\frac{p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}$
$\displaystyle=\frac{p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})}{\left(p_{L}(\sigma^{2})\sigma^{2}-\alpha(\sigma^{2})\right)\sigma^{2}}=\frac{1}{\sigma^{2}}.$
(A.39)

Thus, in total, it holds

$r_{L}\left(\sigma^{2},\pm\,p_{L}(\sigma^{2})\sigma^{2}\right)=\begin{cases}1&\text{if
}\,\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})},\\\
\frac{1}{1-L}\left(1-\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})\sigma^{2}}\right)&\text{if
}\,\frac{\alpha(\sigma^{2})}{Lp_{L}(\sigma^{2})}<\sigma^{2}\leqslant\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L,\\\
\frac{1}{\sigma^{2}}&\text{if
}\,\sigma^{2}>\frac{\alpha(\sigma^{2})}{p_{L}(\sigma^{2})}+1-L.\end{cases}$
(A.40)

## References

  * [1] J. Adler and O. Öktem.  Learned primal-dual reconstruction.  IEEE transactions on medical imaging, 37(6):1322–1332, 2018.
  * [2] G. S. Alberti, E. De Vito, M. Lassas, L. Ratti, and M. Santacesaria.  Learning the optimal tikhonov regularizer for inverse problems.  Advances in Neural Information Processing Systems, 34:25205–25216, 2021.
  * [3] C. Arndt.  Regularization theory of the analytic deep prior approach.  Inverse Problems, 38(11):115005, 2022.
  * [4] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [5] A. Aspri, Y. Korolev, and O. Scherzer.  Data driven regularization by projection.  Inverse Problems, 36(12):125009, 2020.
  * [6] S. Banert, A. Ringh, J. Adler, J. Karlsson, and O. Oktem.  Data-driven nonsmooth optimization.  SIAM Journal on Optimization, 30(1):102–131, 2020.
  * [7] H. Bauermeister, M. Burger, and M. Moeller.  Learning spectral regularizations for linear inverse problems.  In NeurIPS 2020 Workshop on Deep Learning and Inverse Problems.
  * [8] J. Behrmann, W. Grathwohl, R. T. Chen, D. Duvenaud, and J.-H. Jacobsen.  Invertible residual networks.  In International Conference on Machine Learning, pages 573–582. PMLR, 2019.
  * [9] M. Benning and M. Burger.  Modern regularization methods for inverse problems.  Acta Numerica, 27:1–111, 2018.
  * [10] A. Bora, A. Jalal, E. Price, and A. G. Dimakis.  Compressed sensing using generative models.  In International Conference on Machine Learning, pages 537–546. PMLR, 2017.
  * [11] L. Bungert, R. Raab, T. Roith, L. Schwinn, and D. Tenbrinck.  Clip: Cheap lipschitz training of neural networks.  In Scale Space and Variational Methods in Computer Vision: 8th International Conference, SSVM 2021, Virtual Event, May 16–20, 2021, Proceedings, pages 307–319. Springer, 2021.
  * [12] J. Chung, M. Chung, and D. P. O’Leary.  Designing optimal spectral filters for inverse problems.  SIAM Journal on Scientific Computing, 33(6):3132–3152, 2011.
  * [13] A. Denker, M. Schmidt, J. Leuschner, and P. Maass.  Conditional invertible neural networks for medical imaging.  Journal of Imaging, 7(11):243, 2021.
  * [14] S. Dittmer, T. Kluth, P. Maass, and D. Otero Baguer.  Regularization by Architecture: A Deep Prior Approach for Inverse Problems.  Journal of Mathematical Imaging and Vision, 62:456–470, 2020.
  * [15] A. Ebner, J. Frikel, D. Lorenz, J. Schwab, and M. Haltmeier.  Regularization of inverse problems by filtered diagonal frame decomposition.  Applied and Computational Harmonic Analysis, 62:66–83, 2023.
  * [16] A. Ebner and M. Haltmeier.  Plug-and-play image reconstruction is a convergent regularization method.  arXiv preprint arXiv:2212.06881, 2022.
  * [17] H. W. Engl, M. Hanke, and A. Neubauer.  Regularization of inverse problems, volume 375.  Springer Science & Business Media, 1996.
  * [18] K. Gregor and Y. LeCun.  Learning fast approximations of sparse coding.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML’10, pages 399––406, Madison, WI, USA, 2010. Omnipress.
  * [19] A. Hauptmann, F. Lucka, M. Betcke, N. Huynh, J. Adler, B. Cox, P. Beard, S. Ourselin, and S. Arridge.  Model-based learning for accelerated, limited-view 3-d photoacoustic tomography.  IEEE transactions on medical imaging, 37(6):1382–1393, 2018.
  * [20] J. He, Y. Wang, and J. Ma.  Radon inversion via deep learning.  IEEE transactions on medical imaging, 39(6):2076–2087, 2020.
  * [21] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser.  Deep convolutional neural network for inverse problems in imaging.  IEEE Transactions on Image Processing, 26(9):4509–4522, 2017.
  * [22] S. Kabri, A. Auras, D. Riccio, H. Bauermeister, M. Benning, M. Moeller, and M. Burger.  Convergent data-driven regularizations for ct reconstruction.  arXiv preprint arXiv:2212.07786, 2022.
  * [23] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [24] R. Laumont, V. D. Bortoli, A. Almansa, J. Delon, A. Durmus, and M. Pereyra.  Bayesian imaging using plug & play priors: when langevin meets tweedie.  SIAM Journal on Imaging Sciences, 15(2):701–737, 2022.
  * [25] Y. LeCun.  The mnist database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * [26] H. Li, J. Schwab, S. Antholzer, and M. Haltmeier.  Nett: Solving inverse problems with deep neural networks.  Inverse Problems, 36(6):065005, 2020.
  * [27] A. K. Louis.  Inverse und schlecht gestellte Probleme.  Teubner, 1989.
  * [28] S. Lunz, O. Öktem, and C.-B. Schönlieb.  Adversarial regularizers in inverse problems.  Advances in neural information processing systems, 31, 2018.
  * [29] P. Mathé and B. Hofmann.  How general are general source conditions?  Inverse Problems, 24(1):015009, 2008.
  * [30] T. Miyato, T. Kataoka, M. Koyama, and Y. Yoshida.  Spectral normalization for generative adversarial networks.  arXiv preprint arXiv:1802.05957, 2018.
  * [31] S. Mukherjee, S. Dittmer, Z. Shumaylov, S. Lunz, O. Öktem, and C.-B. Schönlieb.  Learned convex regularizers for inverse problems.  arXiv preprint arXiv:2008.02839, 2020.
  * [32] S. Mukherjee, A. Hauptmann, O. Öktem, M. Pereyra, and C.-B. Schönlieb.  Learned reconstruction methods with convergence guarantees: A survey of concepts and applications.  IEEE Signal Processing Magazine, 40(1):164–182, 2023.
  * [33] S. Mukherjee, C.-B. Schönlieb, and M. Burger.  Learning convex regularizers satisfying the variational source condition for inverse problems.  In NeurIPS 2021 Workshop on Deep Learning and Inverse Problems, 2021\.
  * [34] D. Obmann, L. Nguyen, J. Schwab, and M. Haltmeier.  Augmented nett regularization of inverse problems.  Journal of Physics Communications, 5(10):105002, 2021.
  * [35] C. Oh, D. Kim, J.-Y. Chung, Y. Han, and H. Park.  Eter-net: End to end mr image reconstruction using recurrent neural network.  In Machine Learning for Medical Image Reconstruction: First International Workshop, MLMIR 2018, Held in Conjunction with MICCAI 2018, Granada, Spain, September 16, 2018, Proceedings 1, pages 12–20. Springer, 2018\.
  * [36] J.-C. Pesquet, A. Repetti, M. Terris, and Y. Wiaux.  Learning maximally monotone operators for image recovery.  SIAM Journal on Imaging Sciences, 14(3):1206–1237, 2021.
  * [37] A. Rieder.  Keine Probleme mit inversen Problemen: eine Einführung in ihre stabile Lösung.  Springer-Verlag, 2013.
  * [38] E. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, and W. Yin.  Plug-and-play methods provably converge with properly trained denoisers.  In International Conference on Machine Learning, pages 5546–5557. PMLR, 2019.
  * [39] T. Schuster, B. Kaltenbacher, B. Hofmann, and K. S. Kazimierski.  Regularization Methods in Banach Spaces.  De Gruyter, Berlin, Boston, 2012.
  * [40] J. Schwab, S. Antholzer, and M. Haltmeier.  Deep null space learning for inverse problems: convergence analysis and rates.  Inverse Problems, 35(2):025008, 2019.
  * [41] D. Ulyanov, A. Vedaldi, and V. Lempitsky.  Deep image prior.  International Journal of Computer Vision, 128:1867–1888, 2020.
  * [42] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg.  Plug-and-play priors for model based reconstruction.  In 2013 IEEE Global Conference on Signal and Information Processing, pages 945–948. IEEE, 2013.