A Lifted Bregman Formulation for the Inversion of Deep Neural Networks
We propose a novel framework for the regularised inversion of deep neural
networks. The framework is based on the authors' recent work on training
feed-forward neural networks without the differentiation of activation
functions. The framework lifts the parameter space into a higher dimensional
space by introducing auxiliary variables, and penalises these variables with
tailored Bregman distances. We propose a family of variational regularisations
based on these Bregman distances, present theoretical results and support their
practical application with numerical examples. In particular, we present the
first convergence result (to the best of our knowledge) for the regularised
inversion of a single-layer perceptron that only assumes that the solution of
the inverse problem is in the range of the regularisation operator, and that
shows that the regularised inverse provably converges to the true inverse if
measurement errors converge to zero.
Xiaoyu Wang, Martin Benning
# A Lifted Bregman Formulation for the Inversion of Deep Neural Networks

###### Abstract

We propose a novel framework for the regularised inversion of deep neural networks. The framework is based on the authors' recent work on training feed-forward neural networks without the differentiation of activation functions. The framework lifts the parameter space into a higher dimensional space by introducing auxiliary variables, and penalises these variables with tailored Bregman distances. We propose a family of variational regularisations based on these Bregman distances, present theoretical results and support their practical application with numerical examples. In particular, we present the first convergence result (to the best of our knowledge) for the regularised inversion of a single-layer perceptron that only assumes that the solution of the inverse problem is in the range of the regularisation operator, and that shows that the regularised inverse provably converges to the true inverse if measurement errors converge to zero.

Inverse problems, regularisation theory, lifted network training, Bregman distance, perceptron, multi-layer perceptron, variational regularisation, total variation regularisation Article Xiaoyu Wang, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning, Martin Benning
modelling that include (but are not limited to) (variational) Autoencoders (Kingma and Welling, 2013), Normalising Flows (Rezende and Mohamed, 2015; Dinh et al., 2015), Cycle-Consistent Generative Adversarial Networks (Zhu et al., 2017), and Probabilistic Diffusion Models (Sohl-Dickstein et al., 2015; Ho et al., 2020).

While several approaches for the inversion of neural networks have been proposed especially in the context of generative modelling (see for example (Behrmann et al., 2019, 2021) in the context of normalising flows, (Xia et al., 2022) in the context of generative adversarial networks and (Gal et al., 2022) in the context of probabilistic diffusion models), an important aspect, which is often overlooked, is that invertible operations alone are not automatically stable with respect to small variations in the data. For example, computing the solution of the heat equation after a fixed termination time is stable with respect to variations in the initial condition, but estimating the initial condition from the terminal condition of the heat equation is not stable with respect to perturbations in the terminal condition. This issue cannot be resolved without approximation of the inverse with a family of continuous operators, also known as _regularisation_. The research field of _Inverse and Ill-posed Problems_ and its branch _Regularisation Theory_ focus strongly on the stable approximation of ill-posed and ill-conditioned inverses via _regularisations_(Engl et al., 1996) and so-called _variational regularisations_(Scherzer et al., 2009; Benning and Burger, 2018) that are a special class of (nonlinear) regularisations. The optimisation model proposed in (Mahendran and Vedaldi, 2015) can be considered as a variational regularisation method with total variation regularisation; however, the work in (Mahendran and Vedaldi, 2015) is purely empirical, and to the best of our knowledge no works exist that rigorously prove that the proposed approach is a variational regularisation.

In this work, we propose a novel regularisation framework based on lifting with tailored Bregman distances and prove that the proposed framework is a convergent variational regularisation for the inverse problem of estimating the inputs from single-layer perceptrons or the inverse problem of estimating hidden variables in a multi-layer perceptron sequentially. While there has been substantial work in previous years that focuses on utilising neural networks as nonlinear operators in variational regularisation methods (Lunz et al., 2018; Arridge et al., 2019; Schwab et al., 2019; Li et al., 2020; Mukherjee et al., 2021), this is the first work that provides theoretical guarantees for the stable, model-based inversion of neural networks to the best of our knowledge.

Our contributions are three-fold. 1) We propose a novel framework for the regularised inversion of multi-layer perceptrons, respectively feed-forward neural networks, that is based on the lifted Bregman framework recently proposed by the authors in Wang and Benning (2022). 2) We show that for the single-layer perceptron case, the proposed variational regularisation approach is a provably convergent regularisation under very mild assumptions. To our knowledge, this is the first time that an inversion method has been proposed that does not just allow to perform inversion empirically, but for which we can prove that the proposed method is a convergent regularisation method without overly restrictive assumptions such as differentiability of the activation function and the presence of a tangential cone condition. 3) We propose a proximal first-order optimisation strategy to solve the proposed variational regularisation method and present several numerical examples that support the effectiveness of the proposed model-based regularisation approach.

The paper is structured as follows. In Section 2 we introduce the lifted Bregman formulation for the model-based inversion of feed-forward neural networks. In Section 3 we prove that for the single-layer perceptron case the proposed model is a convergent variational regularisation method and provide general error estimates as well as error estimates for a concrete example of a perceptron with ReLU activation function. In Section 4 we discuss how to implement the proposed variational regularisation computationally
for both the single-layer and multi-layer perceptron setting with a generalisation of the primal-dual hybrid gradient method and coordinate descent. Subsequently, we present numerical results that demonstrate empirically that the proposed approach is a model-based regularisation in Section 5, before we conclude this work with a brief section on conclusions and outlook in Section 6.

## 2 Model-based inversion of Feed-forward Networks

Suppose we are given an \(L\)-layer feed-forward neural network \(\mathcal{N}:\mathbb{R}^{n}\times\mathcal{P}\rightarrow\mathbb{R}^{m}\) of the form

\[\mathcal{N}(x,\mathbf{\Theta})=\sigma_{L}(f(\sigma_{L-1}(f(\ldots\sigma_{1}(f(x, \Theta_{1}))\ldots)),\Theta_{L})), \tag{1}\]

for input data \(x\in\mathbb{R}^{n}\) and **pre-trained** parameters \(\mathbf{\Theta}\in\mathcal{P}\). Here, \(\{\sigma_{l}\}_{l=1}^{L}\) denotes the collection of nonlinear activation functions and \(f\) denotes a generic function parametrised by parameters \(\{\Theta_{l}\}_{l=1}^{L}\). For ease of notation, we use \(\mathbf{\Theta}\) to refer to all parameters \(\{\Theta_{l}\}_{l=1}^{L}\). For a given network output \(y\in\mathbb{R}^{m}\), our goal is to solve the inverse problem

\[\mathcal{N}(x,\mathbf{\Theta})=y\]

for the unknown input \(x\in\mathbb{R}^{n}\). We propose to approximate the inverse of this nonlinear, potentially ill-posed inverse problem via the minimisation of a lifted Bregman formulation of the form

\[\begin{pmatrix}x^{\alpha}\\ x_{1}^{\alpha}\\ \vdots\\ x_{L-1}^{\alpha}\end{pmatrix}\in\operatorname*{arg\,min}_{x,x_{1},\ldots,x_{L -1}}\left\{\sum_{l=1}^{L}B_{\Psi_{l}}(x_{l},f(x_{l-1},\Theta_{l}))+\alpha R(x )\right\}\,, \tag{2}\]

where we assume \(x_{0}=x\) and \(x_{L}=y^{\delta}\) for simplicity of notation. The data \(y^{\delta}\) is a perturbed version of \(y\), for which we assume \(B_{\Psi_{L}}(y^{\delta},f(x_{L-1}^{\dagger},\Theta_{L}))\leq\delta^{2}\), for some constant \(\delta\geq 0\) and \(y=\sigma_{L}(f(x_{L-1}^{\dagger},\Theta_{L}))\). The functions \(B_{\Psi_{l}}\) for \(l=1,\ldots,L\) are defined as

\[B_{\Psi_{l}}(x,z)=\frac{1}{2}\|x\|^{2}+\Psi_{l}(x)+\left(\frac{1}{2}\|\cdot\| ^{2}+\Psi_{l}\right)^{*}(z)-\left\langle x,z\right\rangle, \tag{3}\]

for a proper, convex and lower semi-continuous function \(\Psi_{l}\colon\mathbb{R}^{n_{l}}\rightarrow\mathbb{R}\cup\{\infty\}\). The notation \(\left(\frac{1}{2}\|\cdot\|^{2}+\Psi_{l}\right)^{*}\) refers to the convex or Fenchel conjugate of \(\frac{1}{2}\|\cdot\|^{2}+\Psi_{l}\), i.e. \(\left(\frac{1}{2}\|\cdot\|^{2}+\Psi_{l}\right)^{*}(z)=\sup_{y}\langle z,y \rangle-\frac{1}{2}\|y\|^{2}-\Psi_{l}(y)\). Last but not least, the function \(R:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{\infty\}\) is a proper, convex, and lower semi-continuous function that enables us to incorporate a-priori information into the inversion process. The impact of this is controlled by the parameter \(\alpha>0\).

Please note that the functions \(B_{\Psi_{l}}\) have some useful properties and are directly connected to the chosen activation functions \(\{\sigma_{l}\}_{l=1}^{L}\). Following (Wang and Benning, 2022), we observe

\[B_{\Psi_{l}}(x,z)\geq\frac{1}{2}\|\sigma_{l}(z)-x\|^{2}\,,\]
where \(\sigma_{l}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l}}\) is the proximal map with respect to \(\Psi_{l}\), i.e.

\[\sigma_{l}(z)=\operatorname*{arg\,min}_{y\in\mathbb{R}^{n_{l}}}\left\{\frac{1}{ 2}\|y-z\|^{2}+\Psi_{l}(y)\right\}\,,\]

for all \(l\in\{1,\ldots,L\}\). This means that we will solely focus on feed-forward neural networks with **nonlinear activation functions that are proximal maps**.

Another useful property is that the functions \(B_{\Psi_{l}}\) are continuously differentiable with respect to their second argument. If we define \(F_{x}^{l}(z):=B_{\Psi_{l}}(x,z)\), we observe

\[\nabla F_{x}^{l}(z)=\sigma_{l}(z)-x\,. \tag{4}\]

Please note that the family of objective functions \(B_{\Psi_{l}}\) satisfies several other interesting properties; we refer the interested reader to (Wang and Benning, 2022, Theorem 10).

For the remainder of this work, we assume that the parametrised functions \(f\) are affine-linear in the first argument, with parameters \(\Theta_{l}\). A concrete example is the affine-linear transformation \(f(x,\Theta_{l})=W_{l}x+b_{l}\), for a (weight) matrix \(W_{l}\in\mathbb{R}^{n_{l}\times n_{l-1}}\), a (bias) vector \(b_{l}\in\mathbb{R}^{n_{l}}\) and the collection of parameters \(\Theta_{l}=(W_{l},b_{l})\).

In the next section we show that (2) is a variational regularisation method for \(L=1\) and prove a convergence rate with which the solution of (2) converges towards the true input of a perceptron when \(\delta\) converges to zero.

## 3 Convergence analysis and error estimates

In this section we show that the proposed model (2) is a convergent variational regularisation for the specific choice \(L=1\) and the assumption \(f(x,\Theta)=Wx+b\) for \(\Theta=(W,b)\), which reduces (2) to a variational regularisation model for the perceptron case studied in Wang and Benning (2020). In contrast to Wang and Benning (2020) we are not interested in estimating the perceptron parameters \(W\) and \(b\) but assume that these are fixed, and that we study the regularisation operator

\[\mathcal{R}_{\alpha}\colon\,\text{dom}(\Psi)\rightrightarrows\mathbb{R}^{n} \,,\qquad\mathcal{R}_{\alpha}\colon\,y^{\delta}\rightrightarrows x_{\alpha} \in\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\left\{B_{\Psi}\left(y^{ \delta},Wx+b\right)+\alpha R(x)\right\}\,, \tag{5}\]

where \(\text{dom}(\Psi)\) is defined as \(\text{dom}(\Psi):=\{y\in\mathbb{R}^{m}\,|\,\Psi(y)<\infty\}\). We first want to establish under which assumptions (5) is well-defined for all \(y^{\delta}\).

### Well-definedness

For simplicity, we focus on the finite-dimensional setting with network inputs in \(\mathbb{R}^{n}\) and outputs in \(\text{dom}(\Psi)\). However, the following analysis also extends to more general Banach space settings with additional assumptions on the operator \(W\), see for instance (Benning and Burger, 2018, Section 5.1). Following (Benning and Burger, 2018), we assume that \(R\) is non-negative and the polar of a proper function, i.e. \(R=H^{*}\) for a proper function \(H:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}\). Note that this automatically implies convexity of \(R\). Moreover, we assume that \(\Psi\) is a proper, non-negative and convex function that is continuous on \(\text{dom}(\Psi)\), which implies that \(B_{\Psi}\) is proper, non-negative, convex in its second argument and
continuous in its first argument for every \(y^{\delta}\in\text{dom}(\Psi)\). Then, for every \(g\in\text{dom}(\Psi)\) there exists \(x\) with

\[B_{\Psi}(g,Wx+b)+\alpha R(x)<\infty\,.\]

Last but not least, we assume that \(R\) and \(\Psi\) are chosen such that for each \(g\in\text{dom}(\Psi)\) and \(\alpha>0\) we have

\[\|x\|\leq c(a,b,\|g\|),\qquad\text{if}\quad B_{\Psi}(g,Wx+b)\leq a\quad\text{ and}\quad\alpha R(x)\leq d\,,\]

for constants \(a,d\) and a constant \(c\) that depends monotonically non-decreasing on all arguments. With these assumptions, we can then verify the following lemma.

**Theorem 1**.: _Let the assumptions outlined in the previous paragraph be satisfied._

1. _Then, for every_ \(y\in\left\{g\in\text{dom}(\Psi)\,\,\Big{|}\arg\min_{x,\in\mathbb{R}^{n},R(x)< \infty}B_{\Psi}(g,Wx+b)\neq\emptyset\right\}\) _the selection operator_ \[\mathcal{S}(y)=\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\left\{R(x)\,\, \bigg{|}\,x\in\operatorname*{arg\,min}_{\tilde{x}\in\mathbb{R}^{n}}B_{\Psi}(y,W\tilde{x}+b)\right\}\] _is well-defined._
2. _The regularisation operator_ \(\mathcal{R}_{\alpha}\) _as defined in (_2_) is well-defined in the sense that for every_ \(y\in\text{dom}(\Psi)\) _there exists_ \(x_{\alpha}\in\mathbb{R}^{n}\) _with_ \(x_{\alpha}\in\mathcal{R}_{\alpha}(y)\)_. Moreover, the set_ \(\mathcal{R}_{\alpha}(y)\) _is a convex set._
3. _For every sequence_ \(y_{n}\to y\in\text{dom}(\Psi)\) _there exists a subsequence_ \(x_{n_{k}}\in\mathcal{R}_{\alpha}(y_{n_{k}})\) _converging to an element_ \(x^{*}\in\mathcal{R}_{\alpha}(y)\)_._

Proof.: The results follow directly from (Banning and Burger, 2018), Lemma 5.5, Theorem 5.6 and Theorem 5.7. The latter statement originally only implies convergence in the weak-star topology; however, since we are in a finite-dimensional Hilbert space, this automatically implies strong convergence here.

### Error estimates

Having established that (2) is a regularisation operator, we now want to prove that it is also a convergent regularisation operator in the sense of the estimate

\[D_{R}(x^{\dagger},x^{\alpha})\leq C\delta\,, \tag{6}\]

such that

\[\lim_{\delta\to 0}\sup\left\{D_{R}(x^{\dagger},x^{\alpha})\,\,\Big{|}\,x^{ \alpha}\in\mathcal{R}_{\alpha}(y^{\delta}),\,y^{\delta}\in\text{dom}(\Psi),\, B_{\Psi}(y^{\delta},y)\leq\delta^{2}\right\}=0\,.\]

Here, the term \(D_{R}\) denotes the (generalised) Bregman distance (or divergence) (cf. (Bregman, 1967; Kiwiel, 1997)) with respect to \(R\), i.e.

\[D_{R}(x,\tilde{x})=R(x)-R(\tilde{x})-\langle\tilde{q},x-\tilde{x}\rangle\,,\]

for two arguments \(x,\tilde{x}\in\text{dom}(R)\) and a subgradient \(\tilde{q}\in\partial R(\tilde{x})=\{q\in\mathbb{R}^{n}\,|\,R(x)\geq R(\tilde {x})+\langle q,x-\tilde{x}\rangle,\,\forall\,x\in\text{dom}(R)\}\). The vector \(x^{\alpha}\) is a solution of (2) with data \(y^{\delta}\) for which we assume \(B_{\Psi}(y^{\delta},y)\leq\delta^{2}\), and \(C\geq 0\) is a constant. The vector \(x^{\dagger}\) is an element of the selection operator as specified in Lemma 1.1, i.e. \(x^{\dagger}\in\mathcal{S}(y)\) for \(y\in\text{dom}(\Psi)\). Note that \(x^{\dagger}\in\mathcal{S}(y)\) is equivalent to \(x^{\dagger}\) being a \(R\)-minimising vector amongst
all vectors that satisfy \(0=W^{*}\left(\sigma(Wx^{\dagger}+b)-y\right)\), where \(\sigma\) denotes the proximal map with respect to \(\Psi\). This is due to the fact that \(x^{\dagger}\in\arg\min_{\tilde{x}\in\mathbb{R}^{n}}B_{\Psi}(y,W\tilde{x}+b)\) is equivalent to \(0=\nabla B_{\Psi}(y,Wx^{\dagger}+b)=W^{*}\left(\sigma(Wx^{\dagger}+b)-y\right)\). Assuming that \(\sigma(Wx^{\dagger}+b)-y\) does not lie in the nullspace of \(W^{*}\), this further implies \(y=\sigma(Wx^{\dagger}+b)\).

In order to be able to derive error estimates of the form (6), we restrict ourselves to solutions \(x^{\dagger}\) that are in the range of \(\mathcal{R}_{\alpha}\). This means that there exists \(y^{\dagger}\) such that \(x^{\dagger}\in\mathcal{R}_{\alpha}(y^{\dagger})\). Considering the optimality condition of (2) for \(y^{\dagger}\), this implies

\[W^{*}\left(\frac{y^{\dagger}-\sigma(Wx^{\dagger}+b)}{\alpha}\right)\in\partial R (x^{\dagger})\,,\]

which for \(v^{\dagger}:=(y^{\dagger}-\sigma(Wx^{\dagger}+b))/\alpha=(y^{\dagger}-y)/\alpha\) is equivalent to the existence of a source condition element \(v^{\dagger}\) that satisfies the source condition (cf. Engl et al. (1996); Benning and Burger (2018))

\[W^{*}v^{\dagger}\in\partial R(x^{\dagger})\,,\] (SC)

In the following, we verify that the symmetric Bregman distance with respect to \(R\) between a solution of the regularisation operator and the solution of the inverse problem is converging to zero if the error in the data is converging to zero. The symmetric Bregman distance or Jeffreys distance between two vectors \(x\) and \(\tilde{x}\) simply is the sum of two Bregman distances with interchanged arguments, i.e.

\[D_{R}^{\text{symm}}(x,\tilde{x}):=D_{R}(x,\tilde{x})+D_{R}(\tilde{x},x)=\langle x -\tilde{x},q-\tilde{q}\rangle\,,\]

for \(q\in\partial R(x)\) and \(\tilde{q}\in\partial R(\tilde{x})\); hence, an error estimate in the symmetric Bregman distance also implies an error estimate in the classical Bregman distance.

Before we begin our analysis, we recall the concept of the _Jensen-Shannon divergence_(Lin, 1991), which for general proper, convex and lower semi-continuous functions \(F:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}\) generalises to so-called _Burbea-Rao divergences_(Burbea and Rao, 1982a,b; Nielsen and Boltz, 2011) and are defined as follows.

[Burbea-Rao divergence] Suppose \(F:\mathbb{R}^{n}\to\mathbb{R}\cup\{\infty\}\) is a proper, convex and lower semi-continuous function. The corresponding Burbea-Rao divergence is defined as

\[J_{F}(x,\tilde{x}):=\frac{1}{2}\left(F(x)+F(\tilde{x})-2F\left(\frac{x+\tilde {x}}{2}\right)\right)\,, \tag{7}\]

for all \(x,\tilde{x}\in\text{dom}(F)\).

Another important concept that we need in order to establish error estimates is that of Fenchel conjugates (cf. Beck (2017)).

[Fenchel conjugate] The Fenchel (or convex) conjugate \(F^{*}:\mathbb{R}^{n}\to\mathbb{R}\cup\{-\infty,\infty\}\) of a function \(F:\mathbb{R}^{n}\to\mathbb{R}\cup\{-\infty,+\infty\}\) is defined as

\[F^{*}(w)\colon=\sup_{x\in\mathbb{R}^{n}}\langle x,w\rangle-F(x)\,.\]
The Fenchel conjugate that is of particular interest to us is the conjugate of the function \(B_{\Psi}(y,z)\) with respect to the second argument, which we characterise with the following lemma.

**Lemma 1**.: _The Fenchel conjugate of \(F_{y}(z):=B_{\Psi}(y,z)\) with respect to the second argument \(z\) reads_

\[F_{y}^{*}(w)=\left(\frac{1}{2}\|\cdot\|^{2}+\Psi\right)(y+w)-\left(\frac{1}{2} \|\cdot\|^{2}+\Psi\right)(y)\,.\]

Proof.: From the definition of the Fenchel conjugate we observe

\[F_{y}^{*}(w) =\sup_{z\in\mathbb{R}^{m}}\left\langle z,w\right\rangle-F_{y}(z)\] \[=\sup_{z\in\mathbb{R}^{m}}\left\langle z,w\right\rangle-\left( \frac{1}{2}\|\cdot\|^{2}+\Psi\right)(y)-\left(\frac{1}{2}\|\cdot\|^{2}+\Psi \right)^{*}(z)+\left\langle y,z\right\rangle\] \[=-\left(\frac{1}{2}\|\cdot\|^{2}+\Psi\right)(y)+\sup_{z\in \mathbb{R}^{m}}\left\langle z,w+y\right\rangle-\left(\frac{1}{2}\|\cdot\|^{2 }+\Psi\right)^{*}(z)\] \[=-\left(\frac{1}{2}\|\cdot\|^{2}+\Psi\right)(y)+\left(\frac{1}{2 }\|\cdot\|^{2}+\Psi\right)(w+y)\,,\]

which concludes the proof.

Having defined the Burbea-Rao divergence and having established the Fenchel conjugate of \(B_{\Psi}(y,z)\) with respect to the second argument \(z\) for fixed \(y\), we can now present and verify our main result that is motivated by (Benning and Burger, 2011).

**Theorem 2**.: _Suppose \(R\) and \(\Psi\) satisfy the assumptions outlined in Section 3.1. Then, for data \(y^{\delta}\) and \(x^{\dagger}\) that satisfy \(B_{\Psi}(y^{\delta},Wx^{\dagger}+b)\leq\delta^{2}\) with \(\delta\geq 0\), a solution \(x^{\alpha}\in\mathcal{R}_{\alpha}(y^{\delta})\) of the variational regularisation problem (2), and a solution \(x^{\dagger}\) of the perceptron problem \(y=\sigma(Wx^{\dagger}+b)\) that satisfies \(x^{\dagger}\in\mathcal{S}(y)\) and (SC), we observe the error estimate_

\[\begin{split}(1-c)B_{\Psi}(y^{\delta},Wx_{\alpha}+b)+\alpha D_{R} ^{\text{symm}}(x_{\alpha},x^{\dagger})&\leq(1+c)\delta^{2}+ \frac{\alpha^{2}}{c}\|v^{\dagger}\|^{2}\\ &\quad+2cJ_{\Psi}\left(y^{\delta}+\frac{\alpha}{c}v^{\dagger},y^ {\delta}-\frac{\alpha}{c}v^{\dagger}\right)\end{split}, \tag{8}\]

_for a constant \(c\in(0,1]\)._

Proof.: Every solution \(x_{\alpha}\) that satisfies \(x_{\alpha}\in\mathcal{R}_{\alpha}(y^{\delta})\) can equivalently be characterised by the optimality condition

\[W^{*}\left(\sigma(Wx_{\alpha}+b)-y^{\delta}\right)+\alpha p_{\alpha}=0\,,\]

for any subgradient \(p_{\alpha}\in\partial R(x_{\alpha})\). Subtracting \(p^{\dagger}\in\partial R(x^{\dagger})\) from both sides of the equation and taking a dual product with \(x_{\alpha}-x^{\dagger}\) then yields

\[\left\langle\sigma(Wx_{\alpha}+b)-y^{\delta},Wx_{\alpha}-Wx^{\dagger}\right\rangle +\alpha D_{R}^{\text{symm}}(x_{\alpha},x^{\dagger})=-\alpha\langle p^{\dagger}, x_{\alpha}-x^{\dagger}\rangle\,. \tag{9}\]
We easily verify

\[D_{B_{\Psi}(y^{\delta},W\cdot+b)}(x^{\dagger},x_{\alpha})=B_{\Psi}(y^{\delta},Wx^ {\dagger}+b)-B_{\Psi}(y^{\delta},Wx_{\alpha}+b)-\left\langle\sigma(Wx_{\alpha}+b )-y^{\delta},Wx^{\dagger}-Wx_{\alpha}\right\rangle;\]

hence, we can replace \(\left\langle\sigma(Wx_{\alpha}+b)-y^{\delta},Wx_{\alpha}-Wx^{\dagger}\right\rangle\) with \(D_{B_{\Psi}(y^{\delta},W\cdot+b)}(x^{\dagger},x_{\alpha})+B_{\Psi}(y^{\delta}, Wx_{\alpha}+b)-B_{\Psi}(y^{\delta},Wx^{\dagger}+b)\) in (9) to obtain

\[D_{B_{\Psi}(y^{\delta},W\cdot+b)}(x^{\dagger},x_{\alpha})+B_{\Psi}(y^{\delta}, Wx_{\alpha}+b)+\alpha D_{R}^{\text{symm}}(x_{\alpha},x^{\dagger})=B_{\Psi}(y^{ \delta},Wx^{\dagger}+b)-\alpha\langle p^{\dagger},x_{\alpha}-x^{\dagger} \rangle\,.\]

We know \(0\leq D_{B_{\Psi}(y^{\delta},W\cdot+b)}(x^{\dagger},x_{\alpha})\) due to the convexity of \(B_{\Psi}(y^{\delta},W\cdot+b)\), and we also know that (SC) enables us to choose \(p^{\dagger}=W^{*}v^{\dagger}\). Hence, we can estimate

\[B_{\Psi}(y^{\delta},Wx_{\alpha}+b)+\alpha D_{R}^{\text{symm}}(x_{\alpha},x^{ \dagger})\leq B_{\Psi}(y^{\delta},Wx^{\dagger}+b)-\alpha\langle v^{\dagger},Wx _{\alpha}-Wx^{\dagger}\rangle\,.\]

Next, we introduce the constant \(c\in(0,1]\) to split the loss functions \(B_{\Psi}(y^{\delta},Wx_{\alpha}+b)\) and \(B_{\Psi}(y^{\delta},Wx^{\dagger}+b)\) into \((1-c)B_{\Psi}(y^{\delta},Wx_{\alpha}+b)+cB_{\Psi}(y^{\delta},Wx_{\alpha}+b)\) and \((1+c)B_{\Psi}(y^{\delta},Wx^{\dagger}+b)-cB_{\Psi}(y^{\delta},Wx^{\dagger}+b)\), respectively. This means we estimate

\[(1-c)B_{\Psi}(y^{\delta},Wx_{\alpha}+b)+\alpha D_{R}^{\text{symm }}(x_{\alpha},x^{\dagger}) \leq(1+c)B_{\Psi}(y^{\delta},Wx^{\dagger}+b)\] \[\quad+\left\langle\alpha v^{\dagger},Wx^{\dagger}+b\right\rangle -cB_{\Psi}(y^{\delta},Wx^{\dagger}+b)\] \[\quad-\left\langle\alpha v^{\dagger},Wx_{\alpha}+b\right\rangle -cB_{\Psi}(y^{\delta},Wx_{\alpha}+b)\,.\]

Next, we make use of Lemma 1 to estimate

\[\left\langle\alpha v^{\dagger},Wx^{\dagger}+b\right\rangle-cB_{\Psi}(y^{\delta },Wx^{\dagger}+b)\leq c\left(\left(\frac{1}{2}\|\cdot\|^{2}+\Psi\right)\left( y^{\delta}+\frac{\alpha}{c}v^{\dagger}\right)-\left(\frac{1}{2}\|\cdot\|^{2}+ \Psi\right)\left(y^{\delta}\right)\right)\,,\]

and

\[-\langle\alpha v^{\dagger},Wx_{\alpha}+b\rangle-cB_{\Psi}(y^{\delta},Wx_{ \alpha}+b)\leq c\left(\left(\frac{1}{2}\|\cdot\|^{2}+\Psi\right)\left(y^{ \delta}-\frac{\alpha}{c}v^{\dagger}\right)-\left(\frac{1}{2}\|\cdot\|^{2}+ \Psi\right)\left(y^{\delta}\right)\right)\,.\]

Adding both estimates together yields

\[\left\langle\alpha v^{\dagger},Wx^{\dagger}+b\right\rangle-cB_{ \Psi}(y^{\delta},Wx^{\dagger}+b)-\left\langle\alpha v^{\dagger},Wx_{\alpha}+b \right\rangle-cB_{\Psi}(y^{\delta},Wx_{\alpha}+b)\,,\] \[\leq\frac{\alpha^{2}}{c}\|v^{\dagger}\|^{2}+c\left(\Psi\left(y^{ \delta}+\frac{\alpha}{c}v^{\dagger}\right)+\Psi\left(y^{\delta}-\frac{\alpha} {c}v^{\dagger}\right)-2\Psi(y^{\delta})\right)\,,\] \[=\frac{\alpha^{2}}{c}\|v^{\dagger}\|^{2}+2cJ_{\Psi}\left(y^{\delta }+\frac{\alpha}{c}v^{\dagger},y^{\delta}-\frac{\alpha}{c}v^{\dagger}\right)\,,\]

which together with the error bound \(B_{\Psi}(y^{\delta},Wx^{\dagger}+b)\leq\delta^{2}\) concludes the proof.

**Remark 1**.: _We want to emphasise that for continuous \(\Psi\) and \(c>0\) we automatically observe_

\[\lim_{\alpha\to 0}J_{\Psi}\left(y^{\delta}+\frac{\alpha}{c}v^{\dagger},y^{\delta }-\frac{\alpha}{c}v^{\dagger}\right)=0\,,\]
in which case the important question from an error estimate point-of-view is if the term converges quicker to zero than \(\alpha\), as we would need to guarantee \(\lim_{\alpha\to 0}J_{\Psi}\left(y^{\delta}+\frac{\alpha}{c}v^{\dagger},y^{\delta}- \frac{\alpha}{c}v^{\dagger}\right)/\alpha=0\) in order to guarantee that the symmetric Bregman distance in (8) converges to zero for \(\alpha\to 0\)._

**Example 1** (ReLU perceptron).: _Let us consider a concrete example to demonstrate that (5) is a convergent regularisation with respect to the symmetric Bregman distance of \(R\). We know that for \(\sigma(z)=\text{prox}_{\Psi}(z)=\max(0,z)\) to hold true we have to choose \(\Psi(z)=\begin{cases}0&z\in[0,\infty)^{m}\\ \infty&\text{else}\end{cases}\). This means that for \(B_{\Psi}(y^{\delta},z)\) to be well-defined for any \(z\) we require \(y_{i}^{\delta}\geq 0\) for all \(i\in\{1,\ldots,m\}\). In order for the Burbea-Rao divergence to be well-defined, we further require_

\[-\frac{c}{\alpha}y_{i}^{\delta}\leq v_{i}^{\dagger}\leq\frac{c}{\alpha}y_{i}^ {\delta}\,,\]

_for all \(i\in\{1,\ldots,m\}\), or \(\|v^{\dagger}\|_{\infty}\leq(c\|y^{\delta}\|_{\infty}/\alpha)\) in more compact notation. If \(\|v^{\dagger}\|_{\infty}\leq(c\|y^{\delta}\|_{\infty}/\alpha)\) is guaranteed, we observe \(J_{\Psi}\left(y^{\delta}+\frac{\alpha}{c}v^{\dagger},y^{\delta}-\frac{\alpha} {c}v^{\dagger}\right)=0\). Hence, we can simplify the estimate (8) to_

\[\frac{1-c}{\alpha}B_{\Psi}(y^{\delta},Wx_{\alpha}+b)+D_{R}^{\text{ symm}}(x_{\alpha},x^{\dagger})\leq\frac{1+c}{\alpha}\delta^{2}+\frac{\alpha}{c}\|v^{ \dagger}\|^{2}\,,\]

_where we have also divided by \(\alpha\) on both sides of the inequality. If we choose \(\alpha(\delta)=\sqrt{c(1+c)}\delta/\|v^{\dagger}\|\), we obtain the estimate_

\[\frac{(1-c)\|v^{\dagger}\|}{\delta\sqrt{c(1+c)}}B_{\Psi}(y^{\delta},Wx_{\alpha (\delta)}+b)+D_{R}^{\text{symm}}(x_{\alpha(\delta)},x^{\dagger})\leq 2 \sqrt{\frac{1+c}{c}}\|v^{\dagger}\|\delta\,,\]

_as long as we can ensure_

\[\left\|\frac{v^{\dagger}}{\|v^{\dagger}\|}\right\|_{\infty}\leq\sqrt{\frac{c} {1+c}}\left\|\frac{y^{\delta}}{\delta}\right\|_{\infty}\,.\]

_Together with \(D_{R}^{\text{symm}}(x_{\alpha(\delta)},x^{\dagger})\geq D_{R}(x^{\dagger},x_{ \alpha(\delta)})\) we have established an estimate of the form (6), with constant \(C=2\sqrt{\frac{1+c}{c}}\|v^{\dagger}\|\). Hence, we have verified that the variational regularisation method (2) is not only a regularisation method but even a convergent regularisation method in this specific example._

We want to briefly comment on the extension of the convergence analysis to the general case \(L>1\) with the following remark.

**Remark 2**.: _The presented convergence analysis easily extends to a sequential, layer-wise inversion approach. Suppose we have \(L\) layers and begin with the final layer, then we can formulate the variational problem_

\[x_{L-1}^{\alpha}\in\operatorname*{arg\,min}_{x_{L-1}}\left\{B_{\Psi_{L}}(y^{ \delta},W_{L}x_{L-1}+b_{L})+\alpha_{L-1}\Psi_{L-1}(x_{L-1})\right\}\,,\]

_which is also of the form of (5), but where \(R\) has been replaced with \(\Psi_{L-1}\). Alternatively, one can also replace \(\Psi_{L-1}\) with another function \(R_{L-1}\) if good prior knowledge for the auxiliary variable \(x_{L-1}\) exists._
Once we have estimated \(x_{L-1}^{\alpha}\), we can recursively estimate

\[x_{l-1}^{\alpha}\in\operatorname*{arg\,min}_{x_{l-1}}\left\{B_{\Psi_{l}}(x_{l}^{ \alpha},W_{l}x_{l-1}+b_{l})+\alpha_{l-1}\Psi_{l-1}(x_{l-1})\right\}\,,\]

for \(l=L-1,\ldots,2\) and subsequently compute \(x_{\alpha}\) as a solution of (2) but with data \(x_{1}^{\alpha}\) instead of \(y^{\delta}\).

The advantage of such a sequential approach is that every individual regularisation problem is convex and the previously presented theorems and guarantees still apply. The disadvantage is that for this approach to work in theory, we require bounds for every auxiliary variable of the form \(B_{\Psi_{l}}(x_{l}^{\alpha},W_{l}x_{l-1}^{\alpha}+b_{l})\leq\delta_{l}^{2}\), which is a rather unrealistic assumption. Moreover, it is also not realistic to assume that good prior knowledge for the auxiliary variables exist.

Please note that showing that the simultaneous approach (2) is a (convergent) variational regularisation is beyond the scope of this work as it is harder and potentially requires additional assumptions for the following reason. The overall objective function in (2) is no longer guaranteed to be convex with respect to all variables simultaneously, which means that we cannot simply carry over the analysis of the single-layer to the multi-layer perceptron case.

This concludes the theoretical analysis of the perceptron inversion model. In the following section we focus on how to implement (5) and its more general counterpart (2).

## 4 Implementation

In this section, we describe how to computationally implement the proposed variational regularisation for both the single-layer and the multi-layer perceptron setting. More specifically, we show that the proposed variational regularisation can be efficiently solved via a generalised primal-dual hybrid gradient method and a coordinate descent approach.

### Inverting perceptrons

To begin with, we first consider the example of inverting a (single-layer) perceptron. For \(L=1\), Problem (2) reduces to (5), which for a composite regularisation function \(R\circ K\) reads

\[x^{\alpha}\in\operatorname*{arg\,min}_{x}\left\{B_{\Psi}\left(y^{\delta},f(x, \Theta)\right)+\alpha R(Kx)\right\}\,. \tag{10}\]

Here \(K\) is a matrix and \(\alpha R(Kx)\) denotes the regularisation function acting on the argument \(x\). The above Problem (10) can be reformulated to the saddle-point problem

\[\min_{x}\max_{z}B_{\Psi}\left(y^{\delta},f(x,\Theta)\right)+\left\langle z,Kx \right\rangle-\alpha R^{*}(z)\,, \tag{11}\]

where \(R^{*}\) denotes the convex conjugate of \(R\). Computationally, we can then solve the saddle-point problem with a generalised primal-dual hybrid gradient (PDHG) method (Zhu and Chan, 2008; Pock et al., 2009; Esser et al., 2010; Chambolle and Pock, 2011, 2016; Benning and Riis, 2021):

\[x^{k+1} =x^{k}-\tau_{x}\left(\left(\text{prox}_{\Psi}\left(f(x^{k},\Theta )\right)-y^{\delta}\right)\mathcal{J}_{f}^{x}(x^{k},\Theta)+\alpha K^{\top}z^ {k}\right) \tag{12a}\] \[z^{k+1} =\text{prox}_{\tau_{z}R^{*}}\left(z^{k}+\tau_{z}\alpha K\left(2x ^{k+1}-x^{k}\right)\right)\,. \tag{12b}\]
where we alternate between a descent step in the \(x\) variable and an ascent step in the dual variable \(z\). Since (10) is a convex minimisation problem, (12) is guaranteed to converge globally for arbitrary starting point, given that \(\tau_{x}\) and \(\tau_{z}\) are chosen such that \(\tau_{x}\tau_{z}<1/\|K\|\).

In this work, we will focus on the discrete total variation \(\|\nabla x\|_{p,1}\), (Rudin et al., 1992; Chambolle and Lions, 1997), as our regularisation function \(R(Kx)\), but other choices are certainly possible. If we consider a two-dimensional scalar-valued image \(x\in\mathbb{R}^{H\times W}\), we can define a finite forward difference discretisation of the gradient operator \(\nabla:\mathbb{R}^{H\times W}\to\mathbb{R}^{H\times W\times 2}\) as

\[(\nabla x)_{i,j,1}=\begin{cases}x_{i+1,j}-x_{i,j}\;\text{ if }1\leq i<H,\\ 0\;\text{ else,}\end{cases},\quad\quad(\nabla x)_{i,j,2}=\begin{cases}x_{i,j+ 1}-x_{i,j}\;\text{ if }1\leq j<W,\\ 0\;\text{ else.}\end{cases}.\]

The discrete total variation is defined as the \(\ell_{1}\) norm of the \(p\)-norm of the pixel-wise image gradients, i.e.

\[\|\nabla x\|_{p,1}=\sum_{i=1}^{H}\sum_{j=1}^{W}\left|\left(\nabla x\right)_{i,j}\right|_{p}=\sum_{i=1}^{H}\sum_{j=1}^{W}\left((\nabla x)_{i,j,1}^{p}+( \nabla x)_{i,j,2}^{p}\right)^{1/p}\,.\]

For our numerical results we consider the isotropic total variation and consequently choose \(p=2\). Hence for a perceptron with affine-linear transformation \(f(x,\Theta)=Wx+b\), and with \(\sigma=\text{prox}_{\Psi}\) denoting the activation function, the PDHG approach (12) of solving the perceptron inversion problem (5) can be summarised as

\[x^{k+1} =x^{k}-\tau_{x}\left(W^{\top}\left(\sigma(Wx^{k}+b)-y\right)- \alpha\text{div}z^{k}\right)\,, \tag{13a}\] \[z^{k+1} =\text{prox}_{\tau_{z}\|\cdot\|_{2,1}^{*}}\left(z^{k}+\tau_{z} \left(2\alpha\nabla x^{k+1}-\alpha\nabla x^{k}\right)\right)\,. \tag{13b}\]

Please note that we define the discrete approximation of the divergence div such that it satisfies \(\text{div}=-\nabla^{\top}\) in order to be the negative transpose of the discretised finite difference approximation of the gradient in analogy to the continuous case, which is why the sign in (13a) is flipped in comparison to (12a). The proximal map with regards to the convex conjugate of \(\|\cdot\|_{2,1}^{*}\) is simply the argument itself if the maximum of the Euclidean vector-norm per pixel is bounded by one or the projection onto this unit ball.

### Inverting multi-layer perceptrons

We now discuss the implementation of the inversion of multi-layer perceptrons with \(L\) layers as described in (2). Note that in this case in order to minimise for \(x\), we also need to optimise with respect to the auxiliary variables \(x_{1},\ldots,x_{L-1}\).

For the minimisation of (2) we consider an alternating minimisation approach, also known as _coordinate descent_(Beck and Tetruashvili, 2013; Wright, 2015; Wright and Recht, 2022). In this approach we minimise the objective with respect to one variable at a time. In particular, we focus on a semi-explicit coordinate descent algorithm, where we linearise with respect to the smooth functions of the overall objective function. This breaks down the overall minimisation problem into \(L\) sub-problems, where for \(x_{0}\) and each \(x_{l}\) variable
for \(l\in\{1,\ldots,L-1\}\), we have individual minimisation problems of the following form:

\[x_{0}^{k+1} = \tag{14a}\] \[x_{l}^{k+1} \tag{14b}\]

Note that one advantage for adopting this approach is that we exploit that the overall objective function is convex in each individual variable when all other variables are kept fixed. In the following, we will discuss different strategies to computationally solve each sub-problem.

When optimising with respect to the input variable \(x_{0}\), the structure of sub-problem (14a) is identical to the perceptron inversion problem that we have discussed in Section 4.1. Hence, we can approximate \(x_{0}^{k+1}\) with (11), but now with respect to \(x_{1}^{k}\) instead of \(y^{\delta}\), which yields the iteration

\[x_{0}^{k+1} =x_{0}^{k}-\tau_{x_{0}}\left(\left(\text{prox}_{\Psi}\left(f(x_{0} ^{k},\Theta_{1})\right)-x_{1}^{k}\right)\mathcal{J}_{f}^{x}(x_{0}^{k},\Theta_{ 1})+\alpha K^{\top}z^{k}\right)\,, \tag{15a}\] \[z^{k+1} =\text{prox}_{\tau_{z}R^{*}}\left(z^{k}+\tau_{z}\alpha K\left(2x ^{k+1}-x^{k}\right)\right)\,. \tag{15b}\]

For each auxiliary variable \(x_{l}\) with \(l\in\{1,\ldots,L-1\}\), the sub-problem associated with (14b) amounts to solving a proximal gradient step with suitable step-size \(\tau_{x_{l}}\), which we can rewrite to

\[x_{l}^{k+1}=\text{prox}_{\frac{\tau_{x_{l}}}{1+\tau_{x_{l}}}\Psi _{l}}\left(\frac{1}{1+\tau_{x_{l}}}\left(x_{l}^{k}-\tau_{x_{l}}\left(\left( \text{prox}_{\Psi_{l}}\left(f(x_{l}^{k},\Theta_{l+1})\right)-x_{l+1}^{k}\right) \mathcal{J}_{f}^{x}(x_{l}^{k},\Theta_{l+1})\right.\right.\right. \tag{16}\] \[\left.\left.\hskip 14.226378pt-f(x_{l-1}^{k+1},\Theta_{l}) \right)\right)\]

This concludes the discussion on the implementation of the regularised single-layer and multi-layer perceptron inversion. In the next section, we present some numerical results to demonstrate the effectiveness of the proposed approaches empirically.

## 5 Numerical Results

In this section, we present numerical results for the perceptron inversion problem implemented with the PDHG algorithm as outlined in (13), and for the multi-layer perceptron inversion problem implemented with the coordinate descent approach as described in (15) and (16). All results have been computed using PyTorch 3.7 on an Intel Xeon CPU E5-2630 v4.

### The Perceptron

We present results for two experiments: the first one is the perceptron inversion of the image of a circle from the noisy output of the perceptron, where we compare the Landweber regularisation and the total-variation-based variational regularisation (5). For the second experiment, we perform perceptron inversion for samples from the MNIST dataset and compare them with the performance of linear and nonlinear decoders.
**Circle** We begin with the toy example of recovering the image of a circle from noisy measurements of a ReLU perceptron. To prepare the experiment, we generate a circle image \(x^{\dagger}\in\mathbb{R}^{64\times 64}\), as shown in Figure 1. We construct a perceptron with ReLU activation function using random weights and biases where \(W\in\mathbb{R}^{512\times 4096},b\in\mathbb{R}^{512\times 1}\). The weights operates on the column-vector representation of x, where \(x\in\mathbb{R}^{4096\times 1}\). The noise-free data is generated via the forward operation of the model, i.e. \(y=\sigma(Wx^{\dagger}+b)\). We generate noisy data \(y^{\delta}\) by adding Gaussian noise with mean 0 and standard deviation \(0.005\). Note that we clip all the negative values of \(y^{\delta}\) to ensure \(y^{\delta}\in\text{dom}(\Psi)\).

A first attempt to solve this ill-posed perceptron inversion problem is via Landweber regularisation (Landweber, 1951). In Figure 2 we see the reconstructed image obtained with Landweber regularisation in combination with early stopping following Morozov's discrepancy principle Morozov (2012); Engl et al. (1996). Even though the Landweber regularised reconstruction matches the data up to the discrepancy value \(\|\sigma(Wx^{K}+b)-y^{\delta}\|\), the recovered image does not resemble the image \(x^{\dagger}\). We will discuss shortly the reason for this visually poor inversion. In comparison, we see a regularised inversion via the total variation regularisation approach following (13) in Figure 3. The regularisation parameter for this reconstruction is chosen as \(\alpha=1.5\times 10^{-2}\). Both \(x_{0}\) and \(z\) are initialised with zero vectors. The stepsize-parameters are chosen as \(\tau_{x}=1.99/\|W\|_{2}^{2}\) and \(\tau_{z}=1/(8\alpha)\), see (Chambolle, 2004). We stop the iterations when changes in \(x_{0}\) and \(z\) in norm are less than a threshold of \(10^{-5}\) or when we reach the maximum number of iterations, which we set to \(10000\). As shown in Figure 3, the TV-regularisation approach is capable of finding a (visually) more meaningful solution.

To explain why the Landweber iteration performs worse compared to the total variation regularisation for this specific example, we compare the \(\ell_{2}\) norms of each two solutions and the groundtruth image \(x^{\dagger}\). The \(\ell_{2}\) norm of the Landweber solution in Figure 2 measures \(6.58\) while the TV-regularised solution as in Figure 3 and the groundtruth image \(x^{\dagger}\) measure \(25.69\) and \(28.07\) respectively. This is not surprising, as the Landweber iteration is known to converge to a minimal Euclidean norm solution if the noise level converges to zero. On the other hand, when we compare the TV semi-norm of each solution, the groundtruth image in measures \(128.0\), while the Landweber solution in Figure 2 and TV-regularised solution in Figure 3 measure \(707.02\) and \(114.93\) respectively, suggesting that the TV-semi-norm is a more suitable regularisation function for the inversion of cartoon-like images such as \(x^{\dagger}\).

**MNIST** In this second example, we perform perceptron inversion on the MNIST dataset (LeCun et al., 1998). In particular, we consider the following experimental setup. We first train an autoencoder \(\mathcal{A}(x)=\mathcal{D}(\mathcal{E}(x,\mathbf{\Theta}_{\mathcal{E}}), \mathbf{\Theta}_{\mathcal{D}})\), where \(\mathcal{D}(\cdot,\mathbf{\Theta}_{\mathcal{D}})\) and \(\mathcal{E}(\cdot,\mathbf{\Theta}_{\mathcal{E}})\) denotes the decoder and the encoder, parametrised by
parameters \(\mathbf{\Theta}_{\mathcal{D}}\) and \(\mathbf{\Theta}_{\mathcal{E}}\) respectively. We pre-train the autoencoder \(\mathcal{A}\), compute the code \(\mathcal{E}(x,\mathbf{\Theta}_{\mathcal{E}})\) and assign it to the noise-free data variable \(y\), and solve the inverse problem for the input \(x\) from the perturbed code \(y^{\delta}\)

\[\mathcal{E}(x,\mathbf{\Theta}_{\mathcal{E}})=y^{\delta}\.\]

To be more precise, we first train a two-layer fully connected autoencoder \(y=W_{2}(\sigma(W_{1}x+b_{1}))+b_{2}\) using the vanilla stochastic gradient method (SGM) by minimising the mean squared error (MSE) on the MNIST training dataset. We set the code dimension to 100 and use ReLU as the activation function. Hence \(\mathbf{\Theta}_{\mathcal{E}}=(W_{1},b_{1})\) where \(W_{1}\in\mathbb{R}^{784\times 100}\) and \(b_{1}\in\mathbb{R}^{100\times 1}\).

All MNIST images are centred as a means of pre-processing. Algorithmically, we follow (13) to computationally solve (10). The stepsize-parameters are chosen at \(\tau_{x}=1.99/\|W_{1}\|_{2}^{2}\) and \(\tau_{z}=1/(8\alpha)\). We choose the regularisation parameter \(\alpha\) in the range \([10^{-4},10^{-2}]\) and set to \(5\times 10^{-3}\) for all sample images from the training set, and set to \(\alpha=5\times 10^{-2}\) for all sample images from the validation set. These choices work well with regards to the visual quality of the inverted images.

In Figure 4 and Figure 5, we show visualisations of five sample images from the training set, and from the validation set respectively. In comparison, we have also visualised the decoder output. As can be seen, using the code that contains the same compressed information, the inverted images show more clearly

Figure 4: Groundtruth input images from the MNIST training dataset, together with the MNIST validation dataset, together with the corresponding autoencoder output images and inverted input images of the perceptron
defined edges and better visual quality than the decoded outputs. This is to be expected as we compare a nonlinear regularised inversion method with a linear decoder.

### Multi-layer perceptrons

In this section, we present numerical results for inverting multi-layer perceptrons. In particular, we consider feedforward neural networks with convolutional layers (CNN), where in the network architecture two-dimensional convolution operations are used to represent the linear operations in the affine-linear functions \(f(x,\Theta)\). Similar to the experimental design described in Section 5.1, we consider a multi-layer neural network inversion problem where we infer input image \(x\) from a noise perturbed code \(y^{\delta}\).

More specifically, we first train a six-layer convolutional autoencoder on the MNIST training dataset via stochastic gradient method to minimise the MSE. The encoder \(\mathcal{E}(x,\mathbf{\Theta}_{\mathcal{E}})\) consists of two convolutional layers, both with \(4\times 4\) convolutions with stride \(2\), each followed by the application of a ReLU activation function. As image spatial dimension reduce by half, we double the number of feature channels from 8 to 16. We use a fully-connected layer with weights \(W_{3}\in\mathbb{R}^{300\times 784}\) and bias \(b_{3}\in\mathbb{R}^{300\times 1}\) to generate the code. The decoder network first expands the code with an affine-linear transformation with weights \(W_{4}\in\mathbb{R}^{784\times 300}\) and bias \(b_{4}\in\mathbb{R}^{784\times 1}\). This is followed by two layers of transpose convolutionals with kernel size \(4\times 4\), where each is followed by a ReLU activation function. The number of feature channels halves each time as we double the spatial dimension.

Figure 6: Inverted input images of the CNN from the MNIST training dataset, together with from the MNIST validation dataset, together well-trained autoencoder output images and with well-trained autoencoder output images and groudtruth input images
Following the implementation details outlined in Section 4.2, we iteratively compute the update steps (15) and (16). For the PDHG method, we choose the stepsize-parameters as \(\tau_{x}=1.99/\|W_{1}\|_{2}^{2}\) and \(\tau_{z}=1/(8\alpha)\). The initial values \(x_{0}\) and \(z\) are both zero. The update steps stop either after reaching the maximum iterations of \(1500\) or when the improvements on \(x_{0}\) and \(z\) are less than \(10^{-5}\) in norm. For the coordinate descent algorithm, the stepsize-parameters are set to \(\tau_{x_{l}}=1.99/\|W_{l+1}\|_{2}^{2}\) for each layer.

In Figure 6 and Figure 7, we visualise the inverted images, the decoder output images, along with the groundtruth images, from the training dataset and validation dataset respectively. For each image, \(\alpha\) is chosen in the range \([10^{-4},10^{-2}]\) and set at \(9\times 10^{-3}\) for both training sample images and validation sample images for best visual inversion quality.

In Figure 8 we further compare how total variation regularisation and decoder respond to different levels of data noise. The noisy data is produced by adding Gaussian noise to perturb the code of each image. We start with zero mean Gaussian noise with standard deviation 0.33 and gradually reduce the noise level, this translates to decreasing \(\delta^{2}\) from 6.80 down to 0.00.

Please note that for each noise level the regularisation factor \(\alpha\) is manually selected in the range \([10^{-4},10^{-2}]\) for the best PSNR value. As we can see, for the noise level with standard deviation \(0.33\) where \(\delta^{2}\) is at 6.80, the decoder is only capable of producing a blurry distorted output, while the inverted image shows the structure of the digit more clearly. When we decrease the noise level down to \(0.00\), the inverted image becomes more clean-cut while the decoded image is still less sharply defined.

Figure 9 plots the PSNR value of the decoded and inverted image against decreasing noise level. We want to emphasise that it would be more rigorous to compute and compare \(D_{R}^{\text{symm}}(x_{\alpha(\delta)},x^{\dagger})\) as suggested in the error estimate bound in (8), but empirically the PSNR value does also support the notion of a convergent regularisation.

## 6 Conclusions & Outlook

We have introduced a novel variational regularisation framework based on a lifted Bregman formulation for the stable inversion of feed-forward neural networks (also known as multi-layer perceptrons). We have proven that the proposed framework is a convergent regularisation for the single-layer perceptron case under the mild assumption that the inverse problem solution has to be in the range of the regularisation operator. We have derived a general error estimate as well as a specific error estimate for the case that the activation function is the ReLU activation function. We have also addressed the extension of the theory to

Figure 8: Visualisation of the comparison between inverted image and decoded image against various levels of noise. **Top:** Decoded output image from the trained convolutional autoencoder. **Bottom:** Inverted input image from the CNN with total variation regularisation.
the multi-layer perceptron case, which can be carried out sequentially, albeit under unrealistic assumptions. We have discussed implementation strategies to solve the proposed scheme computationally, and presented numerical results for the regularised inversion of the image of a circle and piecewise constant images of hand-written digits from single- and multi-layer perceptron outputs with total variation regularisation.

Despite all the positive achievements presented in this work, the proposed framework also has some limitations. The framework is currently restricted to feed-forward architectures with affine-linear transformations and proximal activation functions. While it is straight-forward to extend the framework to other architectures such as ResNets (He et al., 2016) or U-Nets (Ronneberger et al., 2015), it is not straight-forward to include nonlinear operations that cannot be expressed as proximal maps of convex functions, such as max-pooling. However, for many examples there exist remedies, such as using average pooling instead of max-pooling in the previous example.

An open question is how a convergence theory without restrictive, unrealistic assumptions can be established for the multi-layer case. One issue is the non-convexity of the proposed formulation. A remedy could be the use of different architectures that lead to lifted Bregman formulations that are jointly convex in all auxiliary variables.

And last but not least, one would also like to consider other forms of regularisation, such as iterative regularisation, data-driven regularisations (Kabri et al., 2022), or even combinations of both (Aspri et al., 2020). However, a convergence analysis for such approaches is currently an open problem.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Figure 9: Comparison of PSNR values of total variation-based reconstruction and decoder output per noise level. Each curve reports the change of PSNR value over gradually decreasing levels of Gaussian noise, with \(\delta^{2}\) ranging from 0.00 to 6.80.
## Author Contributions

XW has programmed and contributed all numerical results as well as Section 4 and Section 5. MB has contributed the introduction (Section 1) as well as the theoretical results (Section 3). Both authors have contributed equally to Section 2 and Section 6.

## Acknowledgments

The authors acknowledge support from the Cantab Capital Institute for the Mathematics of Information, the Cambridge Centre for Analysis (CCA) and the Alan Turing Institute (ATI).

## Data Availability Statement

The programming code for this study can be found in the University of Cambridge data repository at 10.17863/CAM.94404.
