Regularization by architecture: A deep prior approach for inverse problems

TL;DR

This paper explores deep image prior (DIP) methods for inverse problems, providing theoretical insights by interpreting DIP as Tikhonov functional optimization and validating findings with numerical experiments.

Contribution

It introduces a novel interpretation of DIP as Tikhonov functional optimization, offering analytical results and bridging theory with practical inverse problem solutions.

Findings

DIP can be viewed as Tikhonov regularization.

Theoretical analysis of specific network designs and operators.

Numerical experiments support the analytical results.

Abstract

The present paper studies so-called deep image prior (DIP) techniques in the context of ill-posed inverse problems. DIP networks have been recently introduced for applications in image processing; also first experimental results for applying DIP to inverse problems have been reported. This paper aims at discussing different interpretations of DIP and to obtain analytic results for specific network designs and linear operators. The main contribution is to introduce the idea of viewing these approaches as the optimization of Tikhonov functionals rather than optimizing networks. Besides theoretical results, we present numerical verifications.

Tables

Table 1. Table 1: Values of ν 𝜈 \nu for which TSVD, Tikhonov and the Soft TSVD are order optimal. For more details see the Proof of Theorem 34 in Appendix II.

Method	$F_{\alpha }\left(\sigma \right)$	$\gamma$	$\nu$
Tikhonov	$\frac{{\sigma }^2}{{\sigma }^2+\alpha }$	$1/2$	$2>\nu >0$
TSVD		$1$	$\nu >0$
Soft TSVD		$1/2$	$\nu >0$

Equations

$$\underset{\Theta }{\min }\Vert A{\phi }_{\Theta }(z)-y^{\delta }{\Vert }^2.$$

$$x(B)=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert Bx-y^{\delta }{\Vert }^2+\alpha R(x),$$

$$x^{k+1}=\varphi \left(W_k{x}^k+b_k\right)$$

$${\phi }_{\Theta }(z)=x^L.$$

$$\Theta =\left\{W_0,..,W_{L-1},b_0,..,b_{L-1}\right\}$$

$$\underset{\Theta }{\min }\sum _{i=1}^N\Vert {\phi }_{\Theta }(y_i^{\delta })-x_i{\Vert }^2.$$

$$y^{\delta }=A{x}^{†}+\tau ,$$

$${\phi }_{\Theta }(z)=x^L,$$

$$x^{k+1}=\varphi \left(W{x}^k+b\right)$$

$$J_B(x)=\frac{1}{2}\Vert Bx-y^{\delta }{\Vert }^2+\alpha R(x),$$

$$\underset{\Theta }{\min }\Vert A{\phi }_{\Theta }(z)-y^{\delta }{\Vert }^2$$

$$=\underset{x\in \mathcal{R}({\phi }_{\cdot }(z))}{\min }\sum _{i=1}^m(x^{\ast }{a}_i-y_i^{\delta }{)}^2,$$

$$\underset{x\in X,\Theta }{\min }\Vert Ax-y^{\delta }{\Vert }^2\text{ s.t. }\Vert x-{\phi }_{\Theta }(z){\Vert }^2=0.$$

$$\mathcal{L}(\Theta ,x,\lambda )=\Vert Ax-y^{\delta }{\Vert }^2+\lambda \Vert x-{\phi }_{\Theta }(z){\Vert }^2,$$

$$\hat{x}=\underset{x\in X}{\mathrm{argmax}}p(x|y^{\delta }),$$

$$\hat{x}=\underset{x\in X}{\mathrm{argmin}}\left\{-\log [p(y^{\delta }|x)]-\log [p(x)]\right\}.$$

$$\hat{x}$$

$$=\underset{x\in X}{\mathrm{argmin}}\Vert A{x}_{\perp }-y^{\delta }{\Vert }^2-\lambda \log p(x_{\perp })-\lambda \log p(x_{\mathcal{N}}|x_{\perp }).$$

$$\displaystyle=\operatorname*{arg\,min}_{x\in X}\ \hbox to0.0pt{$\overbrace{\phantom{\|Ax_{\perp}-y^{\delta}\|^{2}-\lambda\log\ p(x_{\perp})}}^{(I)}$\hss}\|Ax_{\perp}-y^{\delta}\|^{2}\underbrace{-\lambda\log\ p(x_{\perp})-\lambda\log\ p(x_{\mathcal{N}}|x_{\perp})}_{(II)}.$$

$$\underset{B}{\min }F(B)=\underset{B}{\min }\frac{1}{2}\Vert Ax(B)-y^{\delta }{\Vert }^2,$$

$$x(B)=\underset{x}{\mathrm{argmin}}{J}_B(x)=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert Bx-y^{\delta }{\Vert }^2+\alpha R(x).$$

$$x^{k+1}=\underset{\lambda \alpha R}{\mathrm{Prox}}\left(x^k-\lambda B^{\ast }(B{x}^k-y^{\delta })\right).$$

$$x^L=x(B).$$

$$B^{\ell +1}=B^{\ell }-\eta \partial F(B^{\ell }).$$

$$\psi (x,B)=\underset{\lambda \alpha R}{\mathrm{Prox}}\left(x-\lambda B^{\ast }(Bx-y^{\delta })\right)-x.$$

$$\partial F(B)=\partial x(B{)}^{\ast }{A}^{\ast }(Ax(B)-y^{\delta })$$

$$\partial x(B)=-{\psi }_x(x(B),B{)}^{-1}{\psi }_B(x(B),B),$$

$$B^{\ell +1}=B^{\ell }-\eta \partial F(B^{\ell }).$$

$$x(A)=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert Ax-y^{\delta }{\Vert }^2+\frac{\alpha }{2}\Vert x{\Vert }^2,$$

$$R(x)=\frac{1}{2}\Vert x{\Vert }^2.$$

$$x(B)=(B^{\ast }B+\alpha I{)}^{-1}{B}^{\ast }{y}^{\delta },$$

Full text

††journal: Journal of Mathematical Imaging and Vision (JMIV)

1 Introduction

Deep image priors (DIP) were recently introduced in deep learning for some tasks in image processing ulyanov2017dip . Usually, deep learning approaches to inverse problems proceed in two steps. In a first step (training) the parameters $\Theta$ of the deep neural network $\varphi_{\Theta}$ are optimized by minimizing a suitable loss function using large sets of training data. In a second step (application), new data is fed into the network for solving the desired task.

DIP approaches are radically different; they are based on unsupervised training using only a single data point $y^{\delta}$. More precisely, in the context of inverse problems, where we aim at solving ill-posed operator equations $Ax\sim y^{\delta}$, the task of DIP is to train a network $\varphi_{\Theta}(z)$ with parameters $\Theta$ by minimizing the simple loss function

[TABLE]

The minimization is with respect to $\Theta,$ the random input $z$ is kept fixed. After training the solution to the inverse problem is approximated directly by $\hat{x}=\varphi_{\Theta}(z).$

In image processing, common choices for $A$ are the identity operator (denoising) or a projection operator to a subset of the image domain (inpainting). For these applications, it has been observed, that minimizing the functional iteratively by gradient descent methods in combination with a suitable stopping criterion leads to amazing results ulyanov2017dip .

Training with a single data point is the most striking property, which separates DIP from other neural network concepts. One might argue that the astonishing results ulyanov2017dip ; mataev2019deepred ; Cheng_2019_CVPR ; van2018compressed are only possible if the network architecture is fine-tuned to the specific task. This is true for obtaining optimal performance; nevertheless, the presented numerical results perform well even with somewhat generic network architectures such as autoencoders.

We are interested in analyzing DIP approaches for solving ill-posed inverse problems. As a side remark, we note that the applications (denoising, inpainting) mentioned above are modeled by either identity or projection operators, which are not ill-posed in the functional analytical setting louis ; engl ; rieder . Typical examples of ill-posed inverse problems correspond to compact linear operators such as a large variety of tomographic measurement operators or parameter-to-state mappings for partial differential equations.

We aim at analyzing a specific network architecture $\varphi_{\Theta}$ and at interpreting the resulting DIP approach as a regularization technique in the functional analytical setting, and also at proving convergence properties for the minimizers of (1). In particular, we are interested in network architectures, which themselves can be interpreted as a minimization algorithm that solves a regularized inverse problem of the form

[TABLE]

where $R$ is a given convex function and $B$ a learned operator.

In general, deep learning approaches for inverse problems have their own characteristics, and naive applications of neural networks can fail for even the most simple inverse problems, see maass2018trivial . However, there is a growing number of compelling numerical experiments using suitable network designs for some of the toughest inverse problems such as photo-acoustic tomography hauptmann2018model or X-ray tomography with very few measurements jin2017deep ; adler2018learned . Concerning networks based on deep prior approaches for inverse problems, first experimental investigations have been reported, see ulyanov2017dip ; van2018compressed ; mataev2019deepred . Similar as for the above-mentioned tasks in image processing, DIPs for inverse problems rely on two ingredients:

A suitable network design, which leads to our phrase “regularization by architecture”. 2. 2.

Training algorithms for iteratively minimizing (1) with respect to $\Theta$ in combination with a suitable stopping criterion.

In this paper, we present different mathematical interpretations of DIP approaches, and we analyze two network designs in the context of inverse problems in more detail. It is organized as follows: In Section 2, we discuss some relations to existing results and make a short survey of the related literature. In Section 3, we then state different interpretations of DIP approaches and the network architectures that we use, as a basis for the subsequent analysis. We start with a first mathematical result for a trivial network design, which yields a connection to Landweber iterations. We then consider a fully connected feedforward network with $L$ identical layers, which generates a proximal gradient descent for a modified Tikhonov functional. In Section 4, we use this last connection to define the notion of analytic deep prior networks, for which one can strictly analyze its regularization and convergence properties. The key to the theoretical findings is a change of view, which allows for the interpretation of DIP approaches as optimizing families of Tikhonov functionals. Finally, we exemplify our theoretical findings with numerical examples for the standard linear integration operator.

2 Deep prior and related research

We start with a description of general deep prior concepts. Afterwards, we address similarities and differences to other approaches, such as LISTA lista , in more detail.

2.1 The deep prior approach

Present results on deep prior networks utilize feedforward architectures. In general, a feedforward neural network is an algorithm that starts with input $x^{0}=z$, computes iteratively

[TABLE]

$\mbox{for}\ \ k=0,..,L-1$ and outputs

[TABLE]

The parameters of this system are denoted by

[TABLE]

and $\phi$ denotes a non-linear activation function.

In order to highlight one of the unique features of deep image priors, let us first refer to classical generative networks that require training on large data sets.

In this classical setting we are given an operator $A:X\to Y$ between Hilbert spaces $X,Y$, as well as a set of training data $(x_{i},y_{i}^{\delta})$, where $y_{i}^{\delta}$ is a noisy version of $Ax_{i}$ satisfying $\|y_{i}^{\delta}-Ax_{i}\|\leq\delta$. Here the usual deep learning approach is to use a network for direct inversion and the parameters $\Theta$ of the network are obtained by minimizing the loss function

[TABLE]

After training $\Theta$ is fixed and the network is used to approximate the solution of the inverse problem with new data $y^{\delta}$ by computing $x=\varphi_{\Theta}(y^{\delta})$. For a recent survey on this approach and more general deep learning concepts for inverse problems see arridge2019solving .

In general, this approach relies on the underlying assumption, that complex distributions of suitable solutions $x$, e.g., the distribution of natural images, can be approximated by neural networks fista ; unser2007fista ; bruna2013invariant . The parameters $\Theta$ are trained for the specific distribution of training data and are fixed after training. One then expects, that choosing a new data set as input, i.e., $z=y^{\delta}$ will generate a suitable solution to $Ax\sim y^{\delta}$  bora2017compressed . Hence, after training the distribution of solutions is parametrized by the inputs $z$.

In contrast, DIP is an unsupervised approach using only a single data point for training. That means, for given data $y^{\delta}$ and fixed $z$, the parameters $\Theta$ of the network $\varphi_{\Theta}$ are obtained by minimizing the loss function (1). The solution to the inverse problem is then denoted by $\hat{x}=\varphi_{\Theta}(z)$. Hence, deep image priors keep $z$ fixed and aim at parameterizing the solution with $\Theta$. It has been observed in several works ulyanov2017dip ; van2018compressed ; mataev2019deepred ; Cheng_2019_CVPR that this approach indeed leads to remarkable results for problems such as inpainting or denoising.

To some extent, the success of deep image priors is rooted in the careful design of network architectures. For example,  ulyanov2017dip uses a U-Net-like “hourglass” architecture with skip connections, and the amazing results show that such an architecture implicitly captures some statistics of natural images. However, in general, the DIP learning process may converge towards noisy images or undesirable reconstructions. The whole success relies on a combination of the architecture with a suitable optimization method and stopping criterion. Nevertheless, the authors claim the architecture has a positive impact on the exploration of the solution space during the iterative optimization of $\Theta$. They show that the training process descends quickly to “natural-looking” images but requires much more steps to produce noisy images. This is also supported by the theoretical results of saxe2013exact and the observations of  zhang2016understanding , which shows that deep networks can fit noise very well but need more training time to do so. Another paper that hints in this direction is anonymous2019on , which analyzes whether neural networks could have a bias towards approximating low frequencies.

There are already quite a few works that deal with deep prior approaches. Following, we mention the most relevant ones to our work. The original deep image prior article ulyanov2017dip introduces the DIP concept and presents experimental evidence that today’s network architectures are in and of themselves conducive to image reconstruction. Another work van2018compressed explores the applicability of DIP to problems in compressed sensing. Also, mataev2019deepred discusses how to combine DIP with the regularization by denoising approach and Cheng_2019_CVPR explores DIP in the context of stationary Gaussian processes. All of these introduce and discuss variants of DIP concepts; however, neither of them addresses the intrinsic regularizing properties of the network concerning ill-posed inverse problems.

2.2 Deep prior and unrolled proximal gradient architectures

A major part of this paper is devoted to analyzing the DIP approach in combination with an unrolled proximal gradient network $\varphi_{\Theta}$. Hence, there is a natural connection to the well-established analysis of LISTA schemes. Before we sketch the state of research in this field, we highlight the two major differences (loss function, training data) to the present approach. LISTA is based on a supervised training using multiple data points $(x_{i},y_{i}^{\delta})$, $i=1,..N$ where $y_{i}^{\delta}$ is a noisy representation of $Ax_{i}$. The loss function is (3). DIP, however, is based on unsupervised learning using the loss function (1) and a single data point $y^{\delta}$. Hence, DIP with the unrolled proximal gradient network shares the architecture with LISTA, but its concept, as well as its analytic properties, are different. Nevertheless, the analysis we will present in Section 4 will exhibit structures similar to the ones appearing in the LISTA-related literature. Hence we shortly review the major contributions in this field.

Similarities are most visible when considering algorithms and convergence analysis for sparse coding applications xin2016maximal ; moreau2016understanding ; sulam2019multi ; liu2018alista ; sprechmann2015learning . The field of sparse coding makes heavy use of proximal splitting algorithms and, since the advent of LISTA, of trained architectures inspired by truncated versions of these algorithms. In the broadest sense, all of these methods are expressions of “Learning to learn by gradient descent” andrychowicz2016learning . Once more, we would like to emphasize that these results utilize multiple data points while DIP does not require any training data but only one measurement. Another key difference is that we approach the topic from an ill-posed inverse problem perspective, which (a) grounds our approach in the functional analytic realm and (b) considers ill-posed (not only ill-conditioned) problems in the Nashed sense, i.e., allows the treatment of unstable inverses engl . These two points fundamentally differentiate the present approach from traditional compressed sensing considerations which usually deal with (a) finite dimensional formulations and (b) forward operators given by well-conditioned, carefully hand-crafted settings or dictionaries, which are optimized using large sets of training data sprechmann2015learning .

Coming back to LISTA for sparse coding applications, there are many excellent papers moreau2016understanding ; giryes2018tradeoffs ; meinhardt2017learning which are devoted to a strict mathematical analysis of different aspects of LISTA-like approaches. In moreau2016understanding , the authors show under which conditions sparse coding can benefit from LISTA-like trained structures and asks how good trained sparsity estimators can be, given a computational budget. The article giryes2018tradeoffs deals with a similar trade-off proposing the quite exciting, “inexact proximal gradient descent”. The paper chen2018theoretical proposes, based on theoretically founded considerations, a sibling architecture to LISTA. Moreover, papyan2018theoretical argues that deep learning architectures, in general, can be interpreted as multi-stage proximal splitting algorithms.

Finally, we want to point at publications, which address deep learning with only a few data points for training, see, e.g., joergluecke2018 and the references therein. However, they do not address the architectures relevant for our publication, and they do not refer to the specific complications of inverse problems.

3 Deep prior architectures and interpretations

In this section, we discuss different perspectives on deep prior networks, which open the path to provable mathematical results. The first two subsections are devoted to special network architectures, and the last two subsections deal with more general points of view.

3.1 A trivial architecture

We aim at solving ill-posed inverse problems. For a given operator $A,$ the general task in inverse problems is to recover an approximation for $x^{\dagger}$ from measured noisy data

[TABLE]

where $\tau$, with $\|\tau\|\leq\delta,$ describes the noise in the measurement.

The deep image prior approach to inverse problems asks to train a network $\varphi_{\Theta}(z)$ with parameters $\Theta$ and fixed input $z$ by minimizing  $\|A\varphi_{\Theta}(z)-y^{\delta}\|^{2}$ with an optimization method such as gradient descent with early stopping. After training, a final run of the network computes $\hat{x}=\varphi_{\Theta}(z)$ as an approximation to $x^{\dagger}$.

We consider a trivial single-layer network without activation function, see Figure 1. This network simply outputs $\Theta,$ i.e., $\varphi_{\Theta}(z)=\Theta$. In this case, the network parameter $\Theta$ is a vector, which is chosen to have the same dimension as $x$. That means, that training the network by gradient descent of $\|A\varphi_{\Theta}(z)-y^{\delta}\|^{2}=\|A\Theta-y^{\delta}\|^{2}$ with respect to $\Theta$ is equivalent to the classical Landweber iteration, which is a gradient descent method for $\|Ax-y^{\delta}\|^{2}$ with respect to $x$.

Landweber iterations are slowly converging. However, in combination with a suitable stopping rule, they are optimal regularization schemes for diminishing noise level $\delta\rightarrow 0$, louis ; engl ; rieder . Despite the apparent trivialization of the neural network approach, this shows that there is potential in training such networks with a single data point for solving ill-posed inverse problems.

3.2 Unrolled proximal gradient architecture

In this section, we aim at rephrasing DIP, i.e., the minimization of (1) with respect to $\Theta$, as an approach for learning optimized Tikhonov functionals for inverse problems. This change of view, i.e., regarding deep inverse priors as optimization of functionals rather than networks, opens the way for analytic investigations in Section 4.

We use the particular architecture, which was introduced in lista , i.e. a fully connected feedforward network with $L$ layers of identical size,

[TABLE]

where

[TABLE]

The affine linear map $\Theta=(W,b)$ is the same for all layers. The matrix $W$ is restricted to obey $I-W=\lambda B^{*}B$ ($I$ denotes the identity operator) for some $B$ and the bias is determined via $b=\lambda B^{*}y^{\delta}$, see Figure 2. If the activation function of the network is chosen as the proximal mapping of a regularizing functional $\lambda\alpha R$, then $\varphi_{\Theta}(z)$ is identical to the $L$-th iterate of a proximal gradient descent method for minimizing

[TABLE]

see daubechies2004surrogate or Appendix I.

Remark 1

Restricting activation functions to be proximal mappings is not as severe as it might look at first glance. E.g.,  ReLU is the proximal mapping for the indicator function of positive real numbers, and soft shrinkage is the proximal mapping for the modulus function.

This allows the interpretation that every weight update, i.e., every gradient step for minimizing (1) with respect to $\Theta$ or $B$, changes the functional $J_{B}$. Hence, DIP can be regarded as optimizing a functional, which in-turn is minimized by the network. This view is the starting point for investigating convergence properties in Section 4.

3.3 Two perspectives based on regression

The following subsections address more general concepts, which open the way to further analytic investigations, which, however, are not considered further in this paper. The reader interested in the regularization properties for DIP approaches for inverse problems only may jump directly to Section 4.

In this subsection we present two different perspectives on solving inverse problems with the DIP via the minimization of a functional as discussed in the subsection above. The first perspective is based on a reinterpretation of the minimization of the functional (1) in the finite, real setting, i.e. $A\in\mathbb{R}^{m\times n}$. This setting allows us to write

[TABLE]

where $\mathcal{R}(\varphi_{\cdot}(z))$ denotes the range of the network with regard to $\Theta$ for a fixed $z$ and $a_{i}$ the rows of the matrix $A$ as well as $y^{\delta}_{i}$ the entries of the vector $y^{\delta}$. This setting allows for the interpretation that we are solving a linear regression, parameterized by $x$, which is constrained by a deep learning hypothesis space and given by data pairs of the form $(a_{i},y^{\delta}_{i})$.

The second perspective is based on the rewriting of the optimization problem via the method of Lagrange multipliers. We start by considering the constrained optimization problem

[TABLE]

If we now assume that $\varphi$ has continuous first partial derivatives with regard to $\Theta$, the Lagrange functional

[TABLE]

with the correct Lagrange multiplier $\lambda=\lambda_{0}$, has a stationary point at each minimum of the original constraint optimization problem. This gives us a direct connection to unconstrained variational approaches like Tikhonov functionals.

3.4 The Bayesian point of view

The Bayesian approach to inverse problems focuses on computing MAP (maximum a posteriori probability) estimators, i.e. one aims for

[TABLE]

where $p:X\times Y\to\mathbb{R}_{+}\cup\{0\}$ is a conditional PDF. From standard Bayesian theory we obtain

[TABLE]

The setting for inverse problems, i.e. $Ax+\tau=y^{\delta}$ with $\tau\sim\mbox{Normal}(0,\sigma^{2}\mathbb{1}_{Y})$, yields ($\lambda=2\sigma^{2}$)

[TABLE]

We now decompose $x$ into $x_{\perp}:=P_{\mathcal{N}(A)^{\perp}}(x)$,  $and$  $x_{\mathcal{N}}:=P_{\mathcal{N}(A)}(x)$, where $\mathcal{N}(A)$ denotes the nullspace of $A$ and where $P_{\mathcal{N}(A)}(x)$, resp. $P_{\mathcal{N}(A)^{\perp}}(x)$, denotes the orthogonal projection onto $\mathcal{N}(A)$, resp. $\mathcal{N}(A)^{\perp}$. Setting $\hat{x}=(x_{\mathcal{N}},x_{\perp})$ yields

[TABLE]

The data $y^{\delta}$ only contains information about $x_{\perp}$, which in classical regularization is exploited by restricting any reconstruction to ${\cal N}(A)^{\perp}$.

However, if available, $p(x_{\mathcal{N}}|x_{\perp})$ is a measure on how to extend $x_{\perp}$ with an $x_{\perp}\in{\cal N}(A)^{\perp}$ to a suitable $x=(x_{\mathcal{N}},x_{\perp})$. The classical regularization of inverse problems uses the trivial extension by zero, i.e.,  $x=(0,x_{\perp})$, which is not necessarily optimal. If we accept the interpretation that a network can be a meaningful parametrization of the set of suitable solutions $x$, then $p(x)\equiv 0$ for all $x$ not in the range of the network and optimizing the network will indeed yield a non-trivial completion $x=(x_{\mathcal{N}},x_{\perp})$. More precisely (I) can be interpreted to be a deep prior on the measurement and (II) to be a deep prior on the nullspace part of the problem.

4 Deep priors and Tikhonov functionals

In this section, we consider the particular network architecture given by unrolled proximal gradient schemes, see Section 3.2. We aim at embedding this approach into the classical regularization theory for inverse problems. For a strict mathematical analysis, we will introduce the notion of an analytic deep prior network, which then allows interpreting the training of the deep prior network as an optimization of a Tikhonov functional. The main result of this section is Theorem 34, which states that analytic deep priors in combination with a suitable stopping rule are indeed order optimal regularization schemes. Numerical experiments in Section 4.2 demonstrate that such deep prior approaches lead to smaller reconstruction errors when compared with standard Tikhonov reconstructions. The superiority of this approach can be proved, however, only for the rather unrealistic case, that the solution coincides with a singular function of $A$.

4.1 Unrolled proximal gradient networks as deep priors for inverse problems

In this section, we consider linear operators $A$ and aim at rephrasing DIP, i.e., the minimization of (1) with respect to $\Theta$, as a constrained optimization problem. This change of view, i.e., regarding deep inverse priors as an optimization of a simple but constrained functional, rather than networks, opens the way for analytic investigations. We will use an unrolled proximal gradient architecture for the network $\varphi_{\Theta}(z)$ in (1). The starting point for our investigation is the common observation, see combettes2005splitting ; lista or Appendix I, that an unrolled proximal gradient scheme as defined in Section 3.2 approximates a minimizer $x(B)$ of (6). Assuming that a unique minimizer $x(B)$ exists as well as neglecting the difference between $x(B)$ and the approximation $\varphi_{\Theta}(z)$ achieved by the unrolled proximal gradient motivates the following definition of analytic deep priors.

Definition 1

Let us assume that measured data $y^{\delta}\in Y$, a fixed $\alpha>0$, a convex penalty functional $R:X\to\mathbb{R}$, and a measurement operator $A\in\mathcal{L}(X,Y)$ are given. We consider the minimization problem

[TABLE]

subject to the constraint

[TABLE]

We assume that for every $B\in\mathcal{L}(X,Y)$ there is a unique minimizer $x(B)$. We call this constrained minimization problem an analytic deep prior and denote by $x(B)$ the resulting solution to the inverse problems posed by $A$ and $y^{\delta}$.

We can also use this technical definition as the starting point of our consideration and retrieve the neural network architecture by considering the following approach for solving the minimization problem stated in the above definition. Assuming that $R$ has a proximal operator, we can compute $x(B)$, given $B$, via proximal gradient method. I.e., via the (for a suitable choice of $\lambda>0$ and an arbitrary $x^{0}=z\in X$) converging iteration

[TABLE]

Following this iteration for $L$ steps can be seen as the forward pass of a particular architecture of a fully connected feed-forward network with $L$ layers of identical size as described in (4) and (5). The affine linear map given by $\Theta=(W,b)$ is the same for all layers. Moreover, the activation function of the network is given by the proximal mapping of $\lambda\alpha R$, the matrix $W$ is given via $I-W=\lambda B^{*}B$ ($I$ denotes the identity operator), and the bias is determined by $b=\lambda B^{*}y^{\delta}$.

From now on we will assume that the difference between $x^{L}$ and $x(B)$ is negligible, i.e.,

[TABLE]

Remark 2

The task in the DIP approach is to find $\Theta$ (network parameters). Analogously, in the analytic deep prior, we try to find the operator $B$.

We now examine the analytic deep image prior utilizing the proximal gradient descent approach to compute $x(B)$. Therefore we will focus on the minimization of (13) with respect to $B$ for given data $y^{\delta}$ by means of gradient descent.

The stationary points are characterized by $\partial F(B)=0$ and gradient descent iterations with stepsize $\eta$ are given by

[TABLE]

Hence we need to compute the derivative of $F$ with respect to $B$.

Lemma 1

Consider an analytic deep prior with the proximal gradient descent approach as described above. We define

[TABLE]

Then

[TABLE]

with

[TABLE]

which leads to the gradient descent

[TABLE]

This lemma allows to obtain an explicit description of the gradient descent for $B$, which in turn leads to an iteration of functionals $J_{B}$ and minimizers $x(B)$. We will now exemplify this derivation for a rather academic example, which however highlights in particular the differences between a classical Tikhonov minimizer, i.e.

[TABLE]

and the solution of the DIP approach.

4.1.1 Example

In this example we examine analytic deep priors for linear inverse problems $A:X\rightarrow Y$, i.e. $A,B\in{\cal L}(X,Y)$, and

[TABLE]

The rather abstract characterization of the previous section can be made explicit for this setting. Since $J_{B}(x)$ is the classical Tikhonov regularization, which can be solved by

[TABLE]

we can rewrite the analytic deep prior reconstruction as $x(B)$, where $B$ is minimizing

[TABLE]

Lemma 2

Following Lemma 21, assuming $B^{0}=A$ and computing one step of gradient descent to minimize the functional with respect to $B$, yields

[TABLE]

with

[TABLE]

This expression nicely collapses if ${y^{\delta}}({y^{\delta}})^{*}$ commutes with $AA^{*}$. For illustration we assume the rather unrealistic case that $x^{+}=u$, where $u$ is a singular function for $A$ with singular value $\sigma$. The dual singular function is denoted by $v$, i.e. $Au=\sigma v$ and $A^{*}v=\sigma u$ and we further assume, that the measurement noise in $y^{\delta}$ is in the direction of this singular function, i.e., $y^{\delta}=(\sigma+\delta)v$, see Figure 3. In this case, the problem is indeed one-dimensional and we obtain an iteration restricted to the span of $u$, resp. the span of $v$.

Lemma 3

The setting described above yields the following gradient step for the functional in (24):

[TABLE]

with

[TABLE]

and the iteration (28) in-turn results in the sequence $x(B^{\ell})$ with the unique attractive stationary point

[TABLE]

For comparison, the classical Tikhonov regularization would yield $\frac{\sigma}{\sigma^{2}+\alpha}(\sigma+\delta)u$. This is depicted in Figure 4.

4.1.2 Constrained system of singular functions

In the previous example we showed that if we do gradient descent starting from $B_{0}=A$ and assume the rather simple case $y^{\delta}=(\sigma+\delta)v$, we obtain the iteration $B^{\ell+1}=B^{\ell}-c_{\ell}vu^{*}$, i.e. $B^{\ell+1}$ has the same singular functions as $A$ and only one of the singular values is different.

We now analyze the optimization from a different perspective. Namely, we focus on finding directly a minimizer of (13) for a general $y^{\delta}\in Y$, however, we restrict $B$ to be an operator such that $B^{\ast}B$ commutes with $A^{\ast}A$, i.e. $A$ and $B$ share a common system of singular functions. Hence, $B$ has the following representation.

[TABLE]

where $\{u_{i},\sigma_{i},v_{i}\}$ is the singular value decomposition of $A$. That means we restrict the problem to finding optimal singular values $\beta_{i}$ for $B$. In this case we show that a global minimizer exists and that it has interesting properties.

Theorem 4.1

For any $y^{\delta}\in Y$ there exist a global minimizer (in the constrained singular functions setting) of (13) given by $B_{\alpha}=\sum\beta_{i}^{\alpha}v_{i}u_{i}^{*}$ with

[TABLE]

Remark 3

The singular values obtained in Theorem 31 match the ones obtained in the previous section for general $B$ but simple $y^{\delta}=(\sigma+\delta)v$.

Remark 4

The minimizer from Theorem 31 does not depend on $y^{\delta}$, i.e. $\forall:y^{\delta}\in Y$ it holds that $B_{\alpha}$ is a minimizer of (13). The solution to the inverse problem does still depend on $y^{\delta}$ since

[TABLE]

Remark 5

In the original DIP approach, some of the parameters of the network may be similar for different $y^{\delta}$, for example, the parameters of the first layers of the encoder part of the UNet. Other parameters may strongly depend on $y^{\delta}$. In this particular case of the analytic deep prior (constrained system of singular functions) we have a explicit separation of which parameters ($b=\lambda B^{*}y^{\delta}$) depend on $y^{\delta}$ and which do not ($W=I-\lambda B^{*}B$).

From now on we consider the notation $x(B,\,y^{\delta})$ to incorporate the dependency of $x(B)$ on $y^{\delta}$. Following the classical filter theory for order optimal regularization schemes, louis ; rieder ; engl , we obtain the following theorem.

Theorem 4.2

The pseudo inverse $K_{\alpha}:Y\to X$ defined as

[TABLE]

is an order optimal regularization method given by the filter functions

[TABLE]

The regularized pseudoinverse $K_{\alpha}$ is quite similar to the Truncated Singular Value Decomposition (TSVD) but is a softer version because it does not have a jump (see Fig. 5). We call this method Soft TSVD.

The disadvantage of Tikhonov, in this case, is that it damps all singular values, and the disadvantage of TSVD is that it throws away all the information related to small singular values. On the other hand, the Soft TSVD does not damp the higher singular values (similar to TSVD) and does not throw away the information related to smaller singular values but does damp it (similar to Tikhonov). For a comparison of the filter functions, see Table 1. Moreover, what is interesting is how this method comes out from Def. 14, which is stated in terms of the Tikhonov pseudoinverse, and that the optimal singular values do not depend on $y^{\delta}$.

At this point the relation to the original DIP approach becomes more abstract. We considered a simplified network architecture where all layers share the same weights that comes from an iterative algorithm for solving inverse problems. That means, we let the solution to the original inverse problem be the solution of another problem with different operator $B$. The DIP approach in this case is transformed to finding an optimal $B$ and allows us to do the analysis in the functional analysis setting. What we learn from the previous results is that we can establish interesting connections between the DIP approach and the classical Inverse Problems theory. This is important because it shows that deep inverse priors can be used to solve really ill-posed inverse problems.

Remark 6

In the original DIP the input $z$ to the network is chosen arbitrarily and is of minor importance. However, once the weights have been trained for a given $y^{\delta}$, $z$ cannot be changed because it would affect the output of the network, i.e. it would change the obtained reconstruction. In the analytic deep prior the input to the unrolled proximal gradient method is completely irrelevant (assuming an infinite number of layers). After finding the “weights” $B$ a different input will still produce the same solution $\hat{x}=x(B)=\varphi_{\Theta}(z)$.

Remark 6 tells us that there is still a gap between the original DIP and the analytic one. This was expected because of the obvious trivialization of the network architecture but serves as motivation for further research.

4.2 Numerical experiments

We now use the analytic deep inverse prior approach for solving an inverse problem with the following integration operator $A:~{}L^{2}\left(\left[0,1\right]\right)~{}\rightarrow~{}L^{2}\left(\left[0,1\right]\right)$

[TABLE]

$A$ is a linear and compact operator, hence the inverse problem is ill-posed. Let $A_{n}\in\mathbb{R}^{n\times n}$ be a discretization of $A$ and $x^{\dagger}\in\mathbb{R}^{n}$ to be one of its discretized singular vectors $u$. We set the noisy data ${y^{\delta}=A_{n}x^{\dagger}+\delta\tau}$ with ${\tau\sim\mbox{Normal}(0,\mathbb{1}_{n})}$, see Figure 6. A more general example, i.e where $x^{\dagger}$ is not restricted to be a singular function, is also included (Fig. 7).

We aim at recovering $x^{\dagger}$ from $y^{\delta}$ considering the setting established in Def. 14 for ${R(\cdot)=\frac{1}{2}\|\cdot\|^{2}}$. That means that the solution $x$ is parametrized by the operator $B$. Solving the inverse problem is now equivalent to finding optimal $B$ that minimizes the loss function (1) for the single data point $(z,y^{\delta})$.

To find such a $B$, we go back to the DIP and the neural network approach. We write $x(B)$ as the output of the network $\varphi_{\Theta}$ defined in (4) with some randomly initialized input $z$. We optimize with respect to $B$, which is a matrix in the discretized setting, and obtain a minimizer $B_{\text{opt}}$ of (1). For more details, please refer to Appendix III.

In Figure 8 we show some reconstruction results. The first plot of each row contains the true solution $x^{\dagger},$ the standard Tikhonov solution $x(A)$ and the reconstruction obtained with the analytic deep inverse approach $x(B_{\text{opt}})$ after $B$ converged. For each case we provide additional plots depicting:

• •

The true error of the network’s output $x(B)$ after each update of $B$ in a logarithmic scale.

• •

The squared Frobenius norm of $B_{k}-B_{k+1}$ after each update of $B$.

• •

The matrix $B_{\text{opt}}$.

For all choices of $\alpha$ the training of $B$ converges to a matrix $B_{\text{opt}}$, such that $x(B_{\text{opt}})$ has a smaller true error than $x(A)$. In the third plot of each row, one can check that $B$ indeed converges to some matrix $B_{\text{opt}}$, which is shown in the last plot. The networks were trained using gradient descent with $0.05$ as learning rate.

The theoretical findings of the previous subsections allow us to compute, either the exact update (28) for $B$ in the rather unrealistic case that $y^{\delta}=(\sigma+\delta)v$ , or the exact solution $x(B_{\alpha},y^{\delta})$ if we restrict $B$ to have the same system of singular functions as $A$ (Theorem 31). In the numerical experiments we do not consider any of these restrictions and therefore we cannot directly apply our theoretical results. Instead we implement the network approach (see Appendix III) to be able to find $B_{\text{opt}}$ in a more general scenario. Nevertheless, as it can be observed in the last plot of each row in Figure 8, $B_{\text{opt}}$ contains some patterns that reflect, to some extent, that $B$ keeps the same singular system but with different singular values. Namely, $B$ is updated in a similar way as in (28). With the current implementation we could also use more complex regularization functionals $R$, in order to reduce the gap between our analytic approach and the original DIP. This is also a motivation for further research.

5 Summary and conclusion

In this paper, we investigated the concept of deep inverse priors / regularization by architecture. This approach neither requires massive amounts of ground truth / surrogate data, nor pretrained models / transfer learning. The method is based on a single measurement. We started by giving different qualitative interpretations of what regularization is and specifically how regularization by architecture fits into this context.

We followed up with the introduction of the analytic deep prior by explicitly showing how unrolled proximal gradient architectures, allow for a somewhat transparent regularization by architecture. Specifically, we showed that their results can be interpreted as solutions of optimized Tikhonov functionals and proved precise equivalences to regularization techniques. We further investigated this point of view with an academic example, where we implemented the analytic deep inverse prior and tested its numerical applicability. The results confirmed our theoretical findings and showed promising results.

There is obviously, like in deep learning in general, much work to be done in order to have a good understanding of deep inverse priors, but we see much potential in the idea of using deep architectures to regularize inverse problems; especially since an enormous part of the deep learning community is already concerned with the understanding of deep architectures.

Acknowledgements.

S. Dittmer and T. Kluth acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG) within the framework of GRK 2224/1 “Pi 3 : Parameter Identification—Analysis, Algorithms, Applications”. D. Otero Baguer acknowledges the financial support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) –- Projektnummer 276397488 –- SFB 1232, sub-project “P02-Heuristic, Statistical and Analytical Experimental Design”. Peter Maass acknowledges funding by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 765374, sub-project ’Data driven model adaptations of coil sensitivities in MR systems’. The authors would also like to acknowledge the anonymous reviewers for their helpful comments and suggestions.

Appendix I: A reminder on minimization of Tikhonov functionals and the LISTA approach

In this section we consider only linear operators $A$ and we review the well known theory for the Iterative Soft Shrinkage Algorithm (ISTA) as well as the slightly more general Proximal Gradient (PG) combettes2005splitting ; nesterov04convex method for minimizing Tikhonov functionals of the type

[TABLE]

We recapitulate the main steps in deriving ISTA and PG, as far as we need it for our motivation. The necessary first-order condition for a minimizer is given by

[TABLE]

Multiplying with an arbitrary real positive number $\lambda$ and adding $x$ plus rearranging yields

[TABLE]

For convex $R$, the term of the right hand side is inverted by the (single valued) proximal mapping of $\lambda\alpha R$, which yields

[TABLE]

Hence this is a fixed point condition, which is a necessary condition for all minimizers of $J$. Turning the fixed point condition into an iteration scheme yields the PG method

[TABLE]

This structure is also the motivation for LISTA lista approaches where fully connected networks with $L$ internal layers of identical size are used. Moreover, in some versions of LISTA, the affine maps between the layers are assumed to be identical. The values at the $k$-th layer are denoted by $x^{k}$, hence,

[TABLE]

LISTA then trains $(W,b)$ on some given training data. More precisely, it trains two matrices $W=I-\lambda A^{*}A$ and $S=\lambda A^{*}$ such that

[TABLE]

This derivation can be rephrased as follows.

Lemma 4

Let $\varphi_{\Theta}$, $\Theta=(W,b)$, denote a fully connected network with input $x^{0}$ and $L$-internal layers. Further assume, that the activation function is identical to a proximal mapping for a convex functional $\lambda\alpha R:X\rightarrow I\!\!R$. Assume $W$ is restricted, such that $I-W$ is positive definite, i.e., there exists a matrix $B$ such that

[TABLE]

Furthermore, we assume that the bias term is fixed as $b=\lambda B^{*}y^{\delta}$. Then $\varphi_{\Theta}(z)$ is the $L$-th iterate of an ISTA scheme with starting value $x^{0}=z$ for minimizing

[TABLE]

Proof

Follows directly from equation (40).

Appendix II: Proofs

Proof of Lemma 21

$F$ is a functional which maps operators $B$ to real numbers, hence, its derivative is given by

[TABLE]

which follows from classical variational calculus, see, e.g. engl . The derivative of $x(B)$ with respect to $B$ can be computed using the fix point condition for a minimizer of $J_{B}$, namely

[TABLE]

which is equivalent to

[TABLE]

We apply the implicit function theorem and obtain the derivative

[TABLE]

Combining $\partial F(B)$ with $\partial x(B)$ yields the required result.$\square$

Proof of Lemma 2

We start with the the explicit description of the iteration

[TABLE]

with

[TABLE]

The derivative of $x(B)$ with respect to $B$ is a linear map $\partial x(B):{\cal L}(X,Y)\rightarrow X.$ For $\delta B\in{\cal L}(X,Y)$ we obtain

[TABLE]

The adjoint operator is a mapping from $X$ to ${\cal L}(X,Y)$, which can be derived from the defining relation

[TABLE]

Hence,

[TABLE]

Here, ${y^{\delta}}z^{*}\in{\cal L}(X,Y)$ denotes a linear map, which maps an $x\in X$ to $\langle z,x\rangle_{X}\ y^{\delta}$.

First of all, we now aim at determining explicitly $\partial F(B)$ at the starting point of our iteration, i.e., with $B^{0}=A$.

From this follows the rather lengthy expression

[TABLE]

This enables us to compute the update

[TABLE]

as well as the output of the analytic deep prior approach after one iteration of updating $B$ (assuming a suitably chosen $\eta$)

[TABLE]

Proof of Lemma 3

A lengthy computation exploiting $B^{0}=A$ and $\beta_{0}=\sigma$ shows that the singular value $\beta_{\ell}$ of $u$ in the spectral decomposition of $B^{\ell}$ obeys the iteration

[TABLE]

i.e.,

[TABLE]

with

[TABLE]

We will now consider the stability of the fixed points of the sequence $x(B^{\ell})$, i.e., we will analyze the fixed points of the iteration described in (56), that is,

[TABLE]

where

[TABLE]

This iteration in-turn gives you via the Tikhonov filter function, the sequence

[TABLE]

of reconstructions. To find the fixed points of the iteration, we analyze the real roots of $c$, which are

• •

$\beta^{(1)}=\sqrt{\alpha}$,

• •

$\beta^{(2)}=-\sqrt{\alpha}$,

• •

$\beta^{(3)}=\frac{\sigma}{2}+\sqrt{\frac{\sigma^{2}}{4}-\alpha}$, for $\sigma\geq 2\sqrt{\alpha}$ and

• •

$\beta^{(4)}=\frac{\sigma}{2}-\sqrt{\frac{\sigma^{2}}{4}-\alpha}$, for $\sigma\geq 2\sqrt{\alpha}$.

Simple calculations show that

• •

$\partial_{\beta}c(\beta^{(1)})\begin{cases}>0,&\sigma<2\sqrt{\alpha}\\ \leq 0,&\text{otherwise.}\end{cases}$

• •

$\partial_{\beta}c(\beta^{(2)})>0$

• •

$\partial_{\beta}c(\beta^{(3)})>0$, for $\sigma\geq 2\sqrt{\alpha}$ and

• •

$\partial_{\beta}c(\beta^{(4)})>0$ for $\sigma\geq 2\sqrt{\alpha}$.

This leads to the single attractive fixed point $\beta^{(1)}$ for $\sigma<2\sqrt{\alpha}$ and the two attractive fixed points $\beta^{(3)}$ and $\beta^{(4)}$ otherwise. Since,

[TABLE]

we therefore have a unique reconstruction, namely

[TABLE]

$\square$

Proof of Theorem 31

Let $B=\sum\beta_{i}v_{i}u_{i}^{*}$. We want to find $\{\beta_{i}\}$ to minimize

[TABLE]

The Tikhonov solution is given by

[TABLE]

the result of applying the operator $A$ to $x(B)$ is

[TABLE]

and

[TABLE]

Inserting (65) and (66) in (63) yields

[TABLE]

In order to minimize $F(B)$, we should set $\frac{\sigma\beta_{i}}{\beta_{i}^{2}+\alpha}=1$ which implies $\beta_{i}^{2}-\sigma\beta_{i}+\alpha=0$. The roots of the previous equation are $\beta_{i}=\frac{\sigma_{i}}{2}\pm\sqrt{\frac{\sigma_{i}^{2}}{4}-\alpha}$ and they are real only if $\frac{\sigma_{i}^{2}}{4}\geq\alpha$. If it does not hold then $\frac{\alpha\beta_{i}}{\beta_{i}^{2}+\alpha}<1$ and the optimal choice is to find its maximum value which is attained at $\beta_{i}=\sqrt{\alpha}$.

Therefore we set

[TABLE]

and we minimize every term in the sum (67), which means we have found singular values $\{\beta_{i}\}$ that minimize $F(B)$.$\square$

Proof of Theorem 34

In order to prove that $K_{\alpha}$ is a proper order optimal regularization method we need to check if the corresponding filters $F_{\alpha}$ from (34) satisfy the three conditions of optimality louis ; rieder .

These conditions state that a filter $F_{\alpha}:\mathbb{R}\to\mathbb{R}$ is an order optimal regularization filter if $\exists\,\gamma,\,c_{1},\,c_{2},\,c_{3}>0$ such that

$\sup_{\sigma}{\left|F_{\alpha}(\sigma)\sigma^{-1}\right|}\leq c_{1}\alpha^{-\gamma}$ 2. 2.

$\sup_{\sigma}{\left|1-F_{\alpha}(\sigma)\right|\sigma^{\nu}}<c_{2}\alpha^{\gamma\nu}$ 3. 3.

$\forall\alpha>0,\sigma>0:\left|F_{\alpha}(\sigma)\right|\leq c_{3}$

In the following we show that they hold $\forall\nu>0$ with $\gamma=\frac{1}{2},\,c_{1}=\frac{1}{2},\,c_{2}=2^{\nu},\,c_{3}=1$:

1. i

If $\sigma\geq 2\sqrt{\alpha}$

$\sup_{\sigma}{\left|F_{\alpha}(\sigma)\sigma^{-1}\right|}=\sup_{\sigma}{\left|\sigma^{-1}\right|}\leq\frac{1}{2}\alpha^{-\frac{1}{2}}$ 2. 2.

$\sup_{\sigma}{\left|1-F_{\alpha}(\sigma)\right|\sigma^{\nu}}=0\leq\alpha^{\nu}$ 3. 3.

$\forall\alpha>0,\sigma>0:\left|F_{\alpha}(\sigma)\right|=1$ 2. ii

If $\sigma<2\sqrt{\alpha}$

$\sup_{\sigma}{\left|F_{\alpha}(\sigma)\sigma^{-1}\right|}=\frac{1}{2}\alpha^{-\frac{1}{2}}$ 2. 2.

$\begin{aligned} {\sup}_{\sigma}{\left|1-F_{\alpha}(\sigma)\right|\sigma^{\nu}}&=\sup_{\sigma}{\left|\frac{2\sqrt{\alpha}-\sigma}{2\sqrt{\alpha}}\right|\sigma^{\nu}}\\ &\leq 2^{\nu}\alpha^{\frac{\nu}{2}}\end{aligned}$ 3. 3.

$\forall\alpha>0,\sigma>0:\left|F_{\alpha}(\sigma)\right|=\frac{\sigma}{2\sqrt{\alpha}}\leq 1$

$\square$

Appendix III: Numerical experiments

In this section, we provide details about the implementation of the analytic deep inverse prior and the academic example. We start by discretizing the integration operator, which yields the matrix ${A_{n}\in\mathbb{R}^{n\times n}}$, that has $\frac{h}{2}$ on the main diagonal, $h$ everywhere under the main diagonal and [math] above (here $h=\frac{1}{n}$). In our experiments we use $n=200$.

The analytic deep inverse prior network $\varphi^{L}_{\Theta}$ is implemented using Python and Tensorflow tensorflow2015-whitepaper . Initially, we create the matrix $B\in\mathbb{R}^{n\times n}$ and add $L$ fully connected layers to the network, all having the same parameters $\Theta=(W,b)$, with weight matrix $W=I-\lambda B^{T}B$, bias $b=\lambda B^{T}y^{\delta}$ and activation function given by the $\ell_{2}$ proximal operator. That means the network contains in total $4\times 10^{4}$ parameters (the number of components in $B$). For the experiments shown in the paper, the input $z$ is randomly initialized with a small norm and $\lambda$ is $\frac{1}{\mu}$, where $\mu$ is the biggest eigenvalue of $A^{T}A$.

We follow the DIP approach and minimize (1) using gradient descent. To guarantee that $\varphi^{L}_{\Theta}(z)=x(B)$ holds, the network should have thousands of layers, because of the slow convergence of the PG method. This is prohibitive from the implementation point of view. Therefore, we consider only a reduced network with a small number of layers, $L=10,$ and at each iteration we set the input of the network to be the network’s output after the previous iteration. This is equivalent to adding $L$ new identical layers after each update of $B$, with

[TABLE]

and

[TABLE]

where $B_{i}$ refers to the value of $B$ at the $i$-th iteration. After $k$ iterations, we implicitly create a network that has $(k+1)L$ layers, however, each time we update $B$, we back-propagate only through the last $L$ layers.