We propose a non-stationary iterated network Tikhonov (iNETT) method for the solution of ill-posed inverse problems. The iNETT employs deep neural networks to build a data-driven regularizer, and it avoids the difficult task of estimating the optimal regularization parameter. To achieve the theoretical convergence of iNETT, we introduce uniformly convex neural networks to build the data-driven regularizer. Rigorous theories and detailed algorithms are proposed for the construction of convex and uniformly convex neural networks. In particular, given a general neural network architecture, we prescribe sufficient conditions to achieve a trained neural network which is component-wise convex or uniformly convex; moreover, we provide concrete examples of  realizing convexity and uniform convexity in the modern U-net architecture. With the tools of convex and uniformly convex neural networks, the iNETT algorithm is developed and a rigorous convergence analysis is provided. Lastly, we  show applications of the iNETT algorithm in 2D computerized tomography, where numerical examples illustrate the efficacy of the proposed algorithm.

Keywords: iterated network Tikhonov; uniformly convex neural networks;  data-driven regularizer; U-net;  regularization of inverse problem.

MSC2020: 47A52; 65F22; 68T07.

§ INTRODUCTION

Consider the discretized form of an ill-posed linear problem,

\begin{equation}\label{model_eq1}
F \bx  = \by,
\end{equation}

where $X$ and $Y$ denote finite dimensional normed spaces, i.e. $X= \left(\R^N,\|\cdot\|_X\right)$, $Y= \left(\R^M,\|\cdot\|_Y\right)$, and $F : X \to  Y$ is the  discretization of an ill-posed linear operator $\F$. The inverse problem aims to recover $\bx$ from the observed data $\by^\delta$ contaminated by unknown error with bounded norm,

\[
\by^\delta = \by+\boldsymbol{\eta}, \quad \mathrm{where}\ \|\boldsymbol{\eta}\|_Y\leq \delta\,.
\]

As $\by^\delta$ is not necessarily in the range of $F$, i.e. $\by^\delta \notin \operatorname{Rg}(F)$, we consider a variational approach to solve the inverse problem,

\begin{equation}\label{model_eq2}
\bx^\delta\coloneqq \underset{\bx \in X}{\argmin} \|F\bx - \by^\delta  \|^2_Y\,.
\end{equation}

For the ill-posed inverse problem, regularization techniques should be introduced when solving equation (<ref>). The regularization aims to provide prior knowledge and improve stability of the solution. For example, a typical choice of regularization in imaging is the $\ell^p$-norm of $\bx$, with $p\geq 1$, which can be weighted by the Laplacian operator with Dirichlet or Neumann boundary condition <cit.>. Recently, deep learning approaches are introduced to develop data-driven regularization terms in the solutions of inverse problems. In <cit.>, the authors propose a network Tikhonov (NETT) approach, which combines deep neural networks with a Tikhonov regularization strategy. The general form of NETT can be summarized as follows,

\begin{equation}\tag*{NETT}\label{NETT}
x^\delta_\alpha \coloneqq \underset{x \in\dom(\F)\cap \dom\left(\Phi_\Theta\right)}{\argmin} \A(\F x,y^\delta) + \alpha\psi\left(\Phi_{\Theta} (x)\right),%\\
%& \T_{\alpha}(y^\delta) = \A(\F x,y^\delta) + \alpha\Reg(x).
\end{equation}

where $\A(\F x,y^\delta)\geq 0$ is the data-fidelity term which measures misfits between the approximated and measurement data, $\alpha>0$ is a regularization parameter, and $\psi\left(\Phi_\Theta (x)\right)$ is the regularization (or penalty) term including a neural network architecture $\Phi_{\Theta}$. In particular, $\psi$ is a nonnegative functional, and the neural network $\Phi_{\Theta}$ is trained to penalize artifacts in the recovered solution. By training $\Phi_{\Theta}$ in an appropriate way, the neural-network based regularization term is able to capture the feature of solution errors due to data noises and the inexact iterative scheme, so that it can provide penalization on the artifacts of solutions in an adaptive manner. This data-driven regularization strategy shows many advantages in solving inverse problems, and related studies can be found in <cit.> and <cit.> as well.

Motivated by <ref>, we propose an iterated network Tikhonov (iNETT) method which combines the data-driven regularization strategy with an iterated Tikhonov method. In a Tikhonov-like method as <ref>, the regularization parameter $\alpha$ plays an important role since it controls the trade-off between the data-fidelity term and the regularization term. The value of $\alpha$ relies on the noise level $\delta$, and it will affect the proximity of the recovered solution to the minimizer of the data-fidelity term. Poor choices of $\alpha$ can lead to very poor solutions, and it is well known that an accurate estimate of the optimal $\alpha$ is difficult to achieve and it typically relies on heuristic assumptions (e.g., <cit.>). As a result, a natural strategy is to consider an iterated Tikhonov method with non-stationary values of $\alpha$ in the iteration. The non-stationary iterated Tikhonov method is able to avoid exhaustive tuning of the regularization parameter, and it achieves better convergence rates in many applications. For example, we refer the readers to <cit.> for the applications of iterated Tikhonov in Hilbert spaces, and <cit.> in Banach spaces.

Combining the strategy of neural-network based regularizer with the non-stationary iterated Tikhonov method, the  iNETT method has the following general form,

\begin{equation}\tag*{iNETT}\label{iterNETT}
\begin{cases}
\bx^\delta_n:= \underset{\bx \in X}{\argmin} \frac{1}{r}\|F \bx - \by^\delta\|_Y^r + \alpha_n\Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1}),\\
\boldsymbol{\xi}^\delta_{n}:= \boldsymbol{\xi}^\delta_{n-1} - \frac{1}{\alpha_{n}} F^TJ_r\left(F\bx_n^\delta - \by^\delta\right),\\
\bx_0 \in X, \, \boldsymbol{\xi}_0 \in \partial \Reg(\bx_0),
\end{cases}
\end{equation}

where $\Reg\coloneqq \Phi_{\Theta}^{uc} \colon X\to (\R, |\cdot|)$ is a uniformly convex neural network, $\Breg_{\bxi^\delta_{n-1}}(\cdot,\cdot)$ is the Bregman distance induced by $\Reg$ in the direction $\bxi^\delta_{n-1}\in \partial \Reg(\bx^\delta_{n-1})$, $J_r$ denotes the duality map for $r \in (1,\infty)$, and $\{\alpha_n\}_n$ is a sequence of positive real numbers. The value of $\alpha_n$ controls the amount of regularization, and it plays the role of regularization parameter. By taking a decreasing sequence of $\{\alpha_n\}_n$ and considering the standard discrepancy principle as the stopping rule, the <ref> algorithm can automatically determines the amount of regularization. We will provide the details of <ref> in Section <ref>, including a rigorous convergence analysis and many implementation details.

In the formula of <ref>, the neural network $ \Phi_{\Theta}^{uc}$ is employed to build the regularization term, and it is required to be uniformly convex. The property of uniform convexity is demanded in the convergence analysis of the iterated Tikhonov method <cit.>. As a result, another important aspect of the paper is the modeling of convex and uniformly convex neural networks. In Section <ref>, we provide an exact mathematical modeling for the general architecture of neural networks. Our modeling can express the modern convolutional neural networks, where the operations like skip connection and concatenation are included. In Section <ref>, we propose rigorous theories for the convex and uniformly convex neural networks. Given a general neural network $\Phi_\Theta : X \to Z$, we prescribe sufficient conditions to obtain a related neural network which is component-wise convex or uniformly convex. The main idea comes from some recent works on convex neural networks, e.g. <cit.>, but we largely extend them to build modern architectures which can embrace state-of-the-art neural networks. In Section <ref>, we provide particular examples of convex and uniformly convex U-net architectures. The U-net is a convolutional neural network widely used in image processing and related imaging science <cit.>. We give rigorous formulas for the U-net architecture, and explain the approaches to obtain convex and uniformly convex U-net architectures according to the general theories proposed in section <ref>. In Section <ref> and Section <ref>, we provide implementation details as we employ the convex U-net to build a uniformly convex  regularizer for the <ref> algorithm. The proposed method is successfully applied to computerized tomography in Section <ref>. The tool of convex and uniformly convex neural networks is actually a by-product when designing the <ref> algorithm, but it seems more interesting than the algorithm itself. The tool of convex neural networks shall have many interesting applications in the future study.

§ NOTATION AND SETTING

We collect here most of the notations and definitions we will use through this work. As main references, the reader can look at <cit.>. First of all, let us fix $X:= \left(\R^N, \|\cdot \|_X\right)$ and $Y:= \left(\R^M, \|\cdot \|_Y\right)$,  where $N,M \in \mathbb{N}$, and $\|\cdot\|_X$ and $\|\cdot\|_Y$ are some norms on $\R^N$ and $\R^M$, respectively. In the case of standard $\ell^p$ spaces, with $p\geq 1$, then we will indicate the corresponding norm with the usual notation $\|\cdot\|_p$.

We will indicate in bold any finite dimensional (column) vector, e.g. $\bx := (x_1,\ldots, x_N)^T \in \R^N$, where $T$ denotes the transpose operation, and we will use the notation $\bx\leq \hat{\bx}$ meaning that $x_i\leq \hat{x}_i$ for every $i=1,\ldots, N$.  With abuse of language, given a real-valued function $\sigma \colon \R^N \to \R$, we will say that $\sigma$ is monotone nondecreasing if $\sigma(\bx)\leq \sigma(\hat{\bx})$ for every $\bx\leq \hat{\bx}$. In case of a function with multivariate output, $\sigma \colon \R^N \to \R^D$, we will indicate with $\sigma_d \colon \R^N \to \R$ its components, for $d=1,\ldots, D$.

For a fixed $\bz \in Z\coloneqq(\R^D,\|\cdot\|_Z)$, we indicate with  $C(\bz,\cdot) : X \to Z\times X$ the “concatenation” operator, that is,

\begin{equation}\label{concatenation}
C(\bz,\bx) \coloneqq (z_1,\ldots, z_D,x_1,\ldots,x_N)^T.
\end{equation}

Fix now a matrix $F : X \to Y$, which is the discretization of an ill-posed linear operator between Banach spaces. We will assume that, given the unperturbed and observed data $\by \in Y$ and  $\by^\delta\in Y$, respectively, then

\begin{equation}\tag{H0}\label{hypothesis0}
\by \in \operatorname{Rg}(F), \quad \mbox{that is}, \quad 	F\bx = \by \quad \mbox{is solvable},
\end{equation}

$$\by^\delta=\by+\boldsymbol{\eta}  \quad \mbox{where} \quad \|\boldsymbol{\eta}\|_Y\leq \delta.$$

We recall that  a Banach space $Y$ is uniformly smooth if its modulus of smoothness

\begin{equation*}
\rho(\tau) \coloneqq \sup\left\{\frac{\|\by + \tau\hat{\by}\|+\|\by - \tau\hat{\by}\|}{2}-1  \mid \|\by\|=\|\hat{\by}\|=1\right\}, \quad \tau>0
\end{equation*}

satisfies $\lim_{\tau\to 0^+}\rho(\tau)/\tau=0$.  Examples of uniformly smooth spaces are all the $\ell^p$-spaces for $p\in(1,\infty)$.

Given an extended real-valued function, $\Reg  \colon \dom(\Reg)\subseteq X \to (-\infty,+\infty]$, then $\Reg$ is uniformly convex if there exists a nonnegative map $h\colon [0,+\infty)\to [0,+\infty]$ such that $h(s)=0$ if and only if $s=0$ and

\begin{equation*}\label{eq:uni_convex}
\Reg(t\bx + (1-t)\hat{\bx}) + t(1-t)h(\|\bx-\hat{\bx}\|_X) \leq t\Reg(\bx)+ (1-t)\Reg(\hat{\bx}), \qquad \forall t\in[0,1]\mbox{ and } \forall\bx,\hat{\bx} \in \dom(\Reg).
\end{equation*}

Finally, we recall that $\Reg$ is coercive if it is bounded below on bounded sets and

\begin{equation*}
\liminf_{\|\bx\|_X\to \infty}\frac{\Reg(\bx)}{\|\bx\|_X}=\infty.
\end{equation*}

§.§ Bregman distance
Given a convex function
$$\Reg : \dom\left(\Reg\right)\subseteq X \to (-\infty,+\infty],$$
$\Reg$ is called proper if $\dom\left(\Reg\right):= \left\{\bx \in X \, : \, \Reg(\bx)<+\infty\right\}\neq \emptyset$. For every $\hat{\bx} \in \dom\left(\Reg\right)$, a subgradient of $\Reg$ at $\hat{\bx}$ is an element $\bxi$ of the dual space $X^*$ such that

\begin{equation*}
\Reg(\bx) - \Reg(\hat{\bx}) - \langle \bxi, \bx-\hat{\bx}\rangle \geq 0 \quad \forall \bx \in X,
\end{equation*}

where the bracket is the evaluation of $\bxi$ at $\bx-\hat{\bx}$. Clearly, since $X$ is finite dimensional, then $X$ is reflexive and $\langle \cdot, \cdot \rangle$ is the standard inner product, that is,

\begin{equation*}
\langle \bxi, \bx-\hat{\bx}\rangle = \sum_{i=1}^N \xi_i(x_i-\hat{x}_i).
\end{equation*}

The collection of all subgradients of $\Reg$ at $\hat{\bx}$ is denoted by $\partial \Reg (\hat{\bx})$. The subdifferential of $\Reg$ is the multi-valued map $\partial \Reg : \dom\left(\partial \Reg\right)\subseteq X \to 2^{X^*}$ such that

\begin{align*}
& \dom\left(\partial \Reg\right):= \left\{\hat{\bx} \in \dom(\Reg) \, : \, \partial \Reg (\hat{\bx})\neq \emptyset \right\},\\
& \hat{\bx} \mapsto \partial \Reg(\hat{\bx}).
\end{align*}

Let us recall that if $\dom(\Reg)=X$, then $\dom\left(\partial \Reg\right)=X$ (e.g. <cit.>).

Finally, for every $\hat{\bx} \in\dom\left(\partial \Reg\right)$ and $\bxi \in \partial \Reg (\hat{\bx})$, the Bregman distance $\Breg_\bxi(\cdot,\hat{\bx}): X \to [0,+\infty)$ induced by $\Reg$ at $\hat{\bx}$ in the direction $\bxi$ is defined by

\Breg_\bxi(\bx,\hat{\bx}):= \Reg(\bx) - \Reg(\hat{\bx}) - \langle \bxi, \bx-\hat{\bx}\rangle.

It is straightforward to check that if $\Reg$ is uniformly convex then $\Breg_\bxi(\cdot,\hat{\bx})$ is uniformly convex too, for any fixed $\bxi$ and $\hat{\bx}$. Moreover, since $X$ is reflexive, then $\Breg_\bxi(\cdot,\hat{\bx})$ is coercive. See for example <cit.>

We can now introduce the definition of solution of the model problem (<ref>), with respect to the Bregman distance from a reference initial guess.

Fix $\bx_0 \in \dom\left(\partial \Reg\right)$, $\bxi_0 \in\partial \Reg (\bx_0)$. 	An element $\bx^\dagger \in  \dom(\Reg)$ is called a $\Breg_{\bxi_0}$-minimizing solution of (<ref>) if $F\bx^\dagger=\by$ and

\begin{equation*}
\Breg_{\bxi_0}(\bx^\dagger,\bx_0) = \min \left\{  \Breg_{\bxi_0}(\bx,\bx_0) \, : \, \bx \in \dom(\Reg),\, F\bx=\by \right\}.
\end{equation*}

As a last piece of notation, we introduce the duality map. For every fixed $r\in (1,\infty)$, the duality map $J_r \colon X \to 2^{X^*}$ is given by

\begin{equation*}
J_r(\bx) \coloneqq \{ \bxi \in X^* \mid \|\bxi\|=\|\bx\|_X^{r-1} \mbox{ and } \langle \bxi, \bx\rangle =\|\bx\|_X^r \}.
\end{equation*}

In particular, $J_r$ is the subdifferential of the map $\bx\mapsto \frac{\|\bx\|_X^r}{r}$.

If $X$ is an $\ell^p$ space, then $J_r$ is single-valued and for $r=2$ it holds

\begin{equation*}
J_2(\bx) = \operatorname{sgn}(\bx)\frac{|\bx|^{p-1}}{\|\bx\|_2^{p-2}}.
\end{equation*}

More generally, if a Banach space is uniformly smooth, then $J_r$ is single-valued.

In view of the above remark and for the well-posedness of the <ref> method (see Section <ref>), we will assume that

\begin{equation}\tag{H1}\label{hypothesis1}
Y \mbox{ is uniformly smooth}.
\end{equation}

§.§ Neural networks

A neural network is a chain of compositions of  affine operators and  nonlinear operators. For an introduction to neural networks from an applied mathematical point of view, we refer to <cit.>, whereas we refer to <cit.> for a focus on deep learning techniques for inverse problems.  We present here the basic architecture upon which we will devise the neural networks to be implemented in <ref>. There are several many choices for the linear and the nonlinear operators, and each of them generate a different neural network. We do not focus now on those choices, which will be made only later (see  Sections <ref>, <ref> and <ref>), to keep here a more general setting.

First, fix a set of parameters $\Theta:=\{\bb^k;A^{k,j_k};W^k\}_{k=1}^{L},$ where $\bb^k$ are vectors commonly called bias terms, and $A^{k,j_k}$ and $W^k$ are matrices. Second, fix a collection $\{\sigma^k\}_{k=1}^{L}$ of possibly nonlinear operators.

Finally, define  $\Phi_\Theta: X \to Z$ such that

\begin{equation}\tag{NN}\label{DNN}
\begin{cases}
\Phi_\Theta(\bx)\coloneqq \bz \in Z &\mbox{where } \bz= \bz^{L+1},\\
\bz^{k+1}= \sigma^{k}\left(\hat{\bb}^{k} +  W^{k}\hat{\bz}^{k}\right) \in Z^{k+1} &\mbox{for }k=1,\ldots,L,\\
\hat{\bz}^{k}= C^k(\bz^k)\coloneqq
\begin{cases}
\bz^k \mbox{ or} \\
\end{cases} & \mbox{for } i_k\in\{1,\ldots,k\},\\
\hat{\bb}^{k}	= \bb^{k} + A^{k,j_{k}}\bz^{j_{k}} & \mbox{for } j_{k}\in \{1,\ldots, k\},\\
\bz^1\coloneqq \bx \in X,
\end{cases}
\end{equation}

where $Z^{L+1}=Z\coloneqq(\R^D,\|\cdot\|_Z)$ and $Z^k\coloneqq (\R^{D_k},\|\cdot\|_{Z_k})$ are  finite dimensional normed vector spaces, and $C(\cdot,\cdot)$ is the concatenation operator (<ref>).

The operators $C(\bz^{i_{k}},\cdot)$ and $A^{k,j_k}$  represent the skip connections, that is, some of the data in the previous iterations are used in future iterations, skipping intermediate steps.
See Figure <ref> for a visual representation. The integer $L$ is referred to as the depth of the neural network, and $\sigma^{k}\left(\hat{\bb}^{k} +  W^{k}\hat{\bz}^{k}\right)$ as the $k$-th layer. When $L>2$, then the neural network is commonly called deep neural network.

The set $\Theta$ is made by the disjoint union of two subsets, the set of free parameters $\Theta_{free}$ and the set of frozen parameters $\Theta_{frozen}$, that is

\Theta= \Theta_{free}\sqcup \Theta_{frozen}.

The set of free parameters $\Theta_{free}$ is typically initialized to a starting set of values and then it is trained  by minimizing a loss function over a training sample.  Vice-versa, the set of frozen parameters $\Theta_{frozen}$ is fixed and unaffected by the training process. For example, some of the matrices $W^k$ can be fixed to be the identity matrix $I$ or the bias terms $\bb^k$ to be the zero vector $\boldsymbol{0}$. $\Theta_{frozen}$ can be empty, that is, all the parameters are trainable. About the specific training strategy we will employ, see Subsection <ref>.

In the case that $A^{k,j_k}$ is fixed to be the zero matrix and $\hat{\bz}^{k}= \bz^{k}$, for every $k$, then we have a feedforward neural network. For examples of simple architectures of feedforward neural networks of convolutional type, see <cit.>.  For examples of more involved neural networks described by (<ref>), see ResNet <cit.>, DenseNet <cit.> and U-Net <cit.>.

max totalsize=.9.9,center
[auto,node distance =3 cm ,on grid ,
thick ,
inner sep=2pt,
top color = amber,
bottom color = amber
minimum width=2cm,
minimum height=1cm,
top color = azure,
bottom color = azure
input/.style =draw, circle, top color = darkpastelgreen, bottom color =darkpastelgreen, black, text=black,inner sep=0.5pt,
ghost/.style = circle ,top color =white , bottom color = white ,
draw,white , text=white , minimum width =1 cm,
bias/.style = circle ,top color =white , bottom color = white ,
draw, white , text=black , inner sep=0.5pt,
output/.style =draw, circle, top color = amethyst, bottom color =amethyst, black, text=black,inner sep=0.5pt,
inner sep=2pt,

top color = azure,
bottom color = azure,

[square2, top color=white, bottom color=white] (0) [] $\bx$;
[square1] (1/2) [below =1.5cm of 0] $\sigma^1$;
[square2] (1) [below =1.5cm of 1/2] $\bz^2$;
[square1] (3/2) [below =1.5cm of 1] $\sigma^2$;
[square2] (2) [below =1.5cm of 3/2] $\bz^3$;
[square1] (5/2) [below =1.5cm of 2] $\sigma^3$;
[square2] (3) [below =1.5cm of 5/2] $\bz^4$;
[square1] (7/2) [right =2.5cm of 3] $\sigma^4$;
[square2] (4) [right =2.5cm of 7/2] $\bz^5$;
[square1] (9/2) [right =2.5cm of 4] $\sigma^5$;
[square2] (5) [right =2.5cm of 9/2] $\bz^6$;
[square1] (11/2) [above =1.5cm of 5] $\sigma^6$;
[square2] (6) [above =1.5cm of 11/2] $\bz^7$;
[square1] (13/2) [above =1.5cm of 6] $\sigma^7$;
[square2] (7) [above =1.5cm of 13/2] $\bz^8$;
[square1] (15/2) [above =1.5cm of 7] $\sigma^8$;
[square2, top color=darkpastelgreen, bottom color=darkpastelgreen] (8) [above =1.5cm of 15/2] $\bz$;

(0) edge [black,-] node[right] $W^1, \bb^1$ (1/2);
(1/2) edge [black,->] node[right]  (1);
(1) edge [black,-] node[right] $W^2, \bb^2$ (3/2);
(3/2) edge [black,->] node[right]  (2);
(2) edge [black,-] node[right] $W^3, \bb^3$ (5/2);
(5/2) edge [black,->] node[right]  (3);

(0) edge [black,-,bend right=25] node[left] $A^{2,1}$ (1);
(0) edge [black,-,bend right=45] node[left] $A^{4,1}$ (3);

(3) edge [black,-] node[above] $W^4, \bb^4$ (7/2);
(7/2) edge [black,->] node[below]  (4);
(4) edge [black,-] node[above] $W^5,\bb^5$ (9/2);
(9/2) edge [black,->] node[below]  (5);

(5) edge [black,-] node[left] $W^6, \bb^6$ (11/2);
(11/2) edge [black,->] node[left]  (6);
(6) edge [black,-] node[left] $W^7, \bb^7$ (13/2);
(13/2) edge [black,->] node[left]  (7);
(7) edge [black,-] node[left] $W^8, \bb^8$ (15/2);
(15/2) edge [black,->] node[left]  (8);

(3) edge [black,-, bend left= 30] node[below] $A^{6,4}$ (5);

(2) edge [black,double,-] node[below] $C(\bz^3,\bz^7)$ (6);

at (5,-2) [] (layer)
max totalsize=.35.35,center
[auto,node distance =3 cm ,on grid ,
thick ,
inner sep=2pt,
top color = amber,
bottom color = amber
input/.style =draw, circle, top color = azure, bottom color =azure, black, text=black,inner sep=0.5pt,
ghost/.style = circle ,top color =white , bottom color = white ,
draw,white , text=white , minimum width =1 cm,
bias/.style =circle,top color =white , bottom color = white ,
draw, white , text=black , inner sep=0.5pt,
output/.style =draw, circle, top color = azure, bottom color =azure, black, text=black,inner sep=0.5pt,
inner sep=2pt,
top color = amethyst,
bottom color = amethyst,

[state] (E) [] $W^2_1$;
[ghost] (F) [below =of E] ;
[state] (G) [below =of F] $W^2_2$;

[input] (H) [left =of F] $z^2_2$;
[input] (I) [above =of H] $z^2_1$;
[input] (L) [below =of H] $z^2_3$;

[ghost] (A) [right =1.5cm of E] ;
[ghost] (B) [right =1.5cm of G] ;

[ghost] (C) [right =0.3cm of A] ;
[square] (D) [below = of C] $\sigma^2$;

[square] (S1) [right =1.5cm of C] $\sigma^2_1$;
[square] (S2) [below =2cm of S1] $\sigma^2_2$;
[square] (S3) [below =2cm of S2] $\sigma^2_3$;
[square] (S4) [below =2cm of S3] $\sigma^2_4$;

[output] (Y1) [right =1.5cm of S1] $z^3_1$;
[output] (Y2) [below =2cm of Y1] $z^3_2$;
[output] (Y3) [below =2cm of Y2] $z^3_3$;
[output] (Y4) [below =2cm of Y3] $z^3_4$;

[bias] (Eb) [below  =2cm of E] $\hat{b}^2_{\textcolor{dartmouthgreen}{1}}$;
[bias] (Gb) [above  =2cm of G] $\hat{b}^2_{\textcolor{dartmouthgreen}{2}}$;

(H) edge [black,-] node[fill=white, anchor=center, pos=0.5,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{1}\textcolor{azure}{2}}$ (E);
(H) edge [black,-] node[fill=white, anchor=center, pos=0.5,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{2}\textcolor{azure}{2}}$ (G);

(I) edge [black,-] node[fill=white, anchor=center, pos=0.5,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{1}\textcolor{azure}{1}}$ (E);
(I) edge [black,-] node[fill=white, anchor=center, pos=0.70,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{2}\textcolor{azure}{1}}$ (G);

(L) edge [black,-] node[fill=white, anchor=center, pos=0.70,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{1}\textcolor{azure}{3}}$ (E);
(L) edge [black,-] node[fill=white, anchor=center, pos=0.5,font=] $w^{\textcolor{dartmouthgreen}{2}}_{\textcolor{dartmouthgreen}{2}\textcolor{azure}{3}}$ (G);

[black,-] (E) -| (D);
[black,-] (G) -| (D);

[black,-] (D) – ($(D)+(0.5,0)$) |-  (S1);
[black,-] (D) – ($(D)+(0.5,0)$) |-  (S2);
[black,-] (D) – ($(D)+(0.5,0)$) |-  (S3);
[black,-] (D) – ($(D)+(0.5,0)$) |-  (S4);

(S1) edge [black,->] node[]  (Y1);
(S2) edge [black,->] node[]  (Y2);
(S3) edge [black,->] node[]  (Y3);
(S4) edge [black,->] node[]  (Y4);

(Eb) edge [black,o->] node[]  (E);
(Gb) edge [black,o->] node[]  (G);

[thick,dotted]     ($(1.north west)+(-0.28,0.10)$) rectangle ($(3/2.south east)+(0.93,-0.55)$);

[thick,dotted]     ($(5,-2)+(-3,2.5)$) rectangle ($(5,-2)+(3,-2.5)$);
[below=0.5pt] (S) at (layer.south) layer;
[thick,dotted,-Triangle[open]] ($(1)+(1.3,-1)$) – ($(1)+(+2,-1)$);

[below=0.7cm] (S) at (layer.south) Neural Network;

Example of a simple neural network architecture of depth $L=8$, with input $\bx$ and output $\bz$ described by the white and green rectangles, respectively. As it can be visually checked, the operators $A^{2,1}$, $A^{4,1}$ $A^{6,4}$, and $C(\bz^3,\bz^7)$ bring over data of previous iterations skipping some connections. To differentiate the concatenation $C(\cdot,\cdot)$ from the other matrix-vector operators, we use a double edge to denote it. Inside the dotted square, it is depicted a blow-up of the second layer: The rows of the matrix $W^2=(w^2_{ij})\in \R^{2\times 3}$ are represented by the violet circles. Note that in that layer the bias term $\hat{\bb}^2$ is given by $\hat{\bb}^2=\bb^2+ A^{2,1}\bx$.

§ CONVEX AND UNIFORMLY CONVEX NEURAL NETWORKS

We provide some sufficient conditions to guarantee that a neural network, with the architecture described by (<ref>) in Subsection <ref>, is component-wise convex or component-wise uniformly convex. In particular, given any neural network of depth $L$, it is possible to generate a component-wise uniformly convex neural network of depth $L+1$ which embeds the original neural network.

Fix a neural network of depth $L$, $\Phi_\Theta: X \to Z$,  and consider the following constraints:

\begin{align}\tag{C1}\label{hypothesis2}
&	\sigma_d^k \mbox{ is continuous, convex and monotone nondecreasing} \quad  \forall d=1,\ldots,D_k,\, \forall k=1,\dots,L;\\
\tag{C2}\label{hypothesis3}
W^k \mbox{ has nonnegative entries } \forall k=2,\dots,L,\\
A^{k,j_k} \mbox{ has nonnegative entries } \forall j_k\geq 2.
\end{cases}
\end{align}

In the next proposition we show that $\Phi_{\Theta}$ is component-wise convex. This result is a generalization of the input convex neural networks presented in <cit.>, but we largely extend them to build modern architectures which can embrace state-of-the-art neural networks.

Under the constraints (<ref>) and	(<ref>), any neural network described by the set of equations (<ref>) is continuous and component-wise convex.

$\Phi_{\Theta}$ is trivially continuous because composition of continuous  functions. We want to show that given $\Phi_\Theta=(\Phi_{\Theta,1},\ldots, \Phi_{\Theta,D})^T$, then

\begin{equation*}
\Phi_{\Theta,d}(t\bx +(1-t)\bw)\leq t\Phi_{\Theta,d}(\bx) + (1-t)\Phi_{\Theta,d}(\bw) \quad \mbox{for every }d=1,\ldots,D \mbox{ and }\bx,\bw\in X.
\end{equation*}

We proceed by induction. Fix $L=1$ and $\bz^1\coloneqq  t\bx +(1-t)\bw \in X$: Then, by (<ref>),  it holds that

\begin{equation*}
\Phi_{\Theta,d}(\bz^1)=\sigma^1_d\left(\hat{\bb}^{1} +  W^{1}\hat{\bz}^1\right)= \sigma^1_d\left(\bb^{1} + A^{1,1}\bz^1 +  W^{1}C^1(\bz^1)\right).
\end{equation*}

Let us remember that $C^1$ is such that

\begin{cases}
C^1 : X \to X,\\
\end{cases}\quad \mbox{or} \quad \begin{cases}
C^1 : X \to X\times X,\\
\end{cases}

In the latter case,

\begin{align*}
C^1(\bz^0)&= (t\bx +(1-t)\bw, t\bx +(1-t)\bw)^T\\
&= t(\bx,\bx)^T+ (1-t)(\bw,\bw)^T\\
&=tC^1(\bx) + (1-t)C^1(\bw).
\end{align*}

Then, by the linearity of $A^{1,1}$ and $W^1$, and independently of the choice of $C^1$, it is clear that

\begin{align*}
\bb^{1} + A^{1,1}\bz^1 +  W^{1}C^1(\bz^1)&= t(\bb^{1} + A^{1,1}\bx +  W^1C^1(\bx)) + (1-t)(\bb^{1} + A^{1,1}\bw + W^1C^1(\bw))\\
&= t\bu + (1-t)\bv,
\end{align*}

\bu\coloneqq\bb^{1} + A^{1,1}\bx +  W^1C^1(\bx), \qquad \bv\coloneqq\bb^{1} + A^{1,1}\bw + W^1C^1(\bw).

Using now the convexity of $\sigma^1_d$, we obtain

\begin{align*}
\Phi_{\Theta,d}(t\bx +(1-t)\bw) &=\sigma^1_d\left(\bb^{1} + A^{1,1}\bz^1 +  W^{1}C^1(\bz^1)\right)\\ &= \sigma^1_d\left( t\bu + (1-t)\bv\right)\\
&\leq t\sigma^1_d(\bu) + (1-t)\sigma^0_d(\bv)\\
&= t\Phi_{\Theta,d}(\bx) + (1-t)\Phi_{\Theta,d}(\bw),
\end{align*}

namely, $\Phi_\Theta$ is component-wise convex. Assume now that every neural network $\Phi_{\Theta}$ described by (<ref>) is component-wise convex when the depth $L$ is $L=1,\ldots, n-1$, for any fixed $n\geq2$. We need to prove that $\Phi_{\Theta}$ is component-wise convex when $L=n$. Fix again $\bz^1\coloneqq t\bx +(1-t)\bw$. Observe that, at each step $k=1,\ldots,L-1$, every $\bz^{k+1}$ in (<ref>) is a neural network $\Phi_{\Theta}^{k} : X \to Z^{k+1}$ of depth $L_{k}\leq L-1$ and with the same initial input $\bz^1$, that is $\bz^{k+1}=\Phi_{\Theta}^{k}(\bz^1)$. Then, by the induction hypothesis,

\begin{equation}\label{eq:induction}
\Phi_{\Theta}^{k}(\bz^1)=\Phi_{\Theta}^{k}(t\bx +(1-t)\bw)\leq t\Phi_{\Theta}^{k}(\bx) +(1-t)\Phi_{\Theta}^{k}(\bw) \qquad \forall\, k=1,\ldots,L-1.
\end{equation}

\begin{align*}
\Phi_{\Theta,d}(\bz^1)=\bz^{L+1}&= \sigma^{L}_d\left(\hat{\bb}^{L} +  W^{L}\hat{\bz}^{L}\right)\\
&= \sigma^{L}_d\left(\bb^{L} + A^{L,j_{L}}\bz^{j_{L}} +  W^{L}C^{L}(\bz^{L})\right),\\
&= \sigma^{L}_d\left(\bb^{L} + A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bz^1) +  W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bz^1))\right).
\end{align*}

We need to take care of the concatenation operator. Again, there are two cases:
\begin{cases}
C^{L}\colon  Z^{L} \to Z^{L},\\
C^{L}(\Phi_{\Theta}^{L-1}(\bz^1)) = \Phi_{\Theta}^{L-1}(\bz^1),
\end{cases}\quad \mbox{or} \quad
\begin{cases}
C^{L}\colon  Z^{L} \to Z^{i_{L}}\times Z^{L},\\
\end{cases}

where $i_{L} \in \{1,\ldots,L\}$, and with the convention that $\Phi_{\Theta}^0(\bz^1)=\bz^1$. Clearly,

C(\bu_1,\bu_2)\leq C(\bv_1,\bv_2) \quad \mbox{for every}\quad \bu_1\leq \bv_1,\;  \bu_2\leq \bv_2.

Therefore, if we are in the latter case, by  (<ref>)

\begin{align*}
C(\Phi_{\Theta}^{i_{L-1}}(\bz^1),\Phi_{\Theta}^{L-1}(\bz^1)) &= C(\Phi_{\Theta}^{i_{L-1}}(t\bx +(1-t)\bw),\Phi_{\Theta}^{L-1}(t\bx +(1-t)\bw))\nonumber\\
&\leq C(t\Phi_{\Theta}^{i_{L-1}}(\bx) + (1-t)\Phi_{\Theta}^{i_{L-1}}(\bw),t\Phi_{\Theta}^{L-1}(\bx) + (1-t)\Phi_{\Theta}^{L-1}(\bw))\nonumber\\
&=tC(\Phi_{\Theta}^{i_{L-1}}(\bx),\Phi_{\Theta}^{L-1}(\bx)) + (1-t)C(\Phi_{\Theta}^{i_{L-1}}(\bw),\Phi_{\Theta}^{L-1}(\bw))\\
&= tC^{L}(\Phi_{\Theta}^{L-1}(\bx)) +(1-t)C^{L}(\Phi_{\Theta}^{L-1}(\bw)).
\end{align*}

As a consequence, independently from the case we are in,

\begin{equation*}
C^{L}(\Phi_{\Theta}^{L-1}(\bz^1)) \leq tC^{L}(\Phi_{\Theta}^{L-1}(\bx)) +(1-t)C^{L}(\Phi_{\Theta}^{L-1}(\bw)),
\end{equation*}

and then, by (<ref>),

\begin{equation}\label{prop1:eq1}
W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bz^1))\leq tW^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bx)) + (1-t)W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bw)).
\end{equation}

Now, instead, if $j_{L}=1$, then $\Phi_{\Theta}^{j_{L}-1}(\bz^1)=\bz_1$ and by linearity

\begin{equation}\label{prop1:eq2}
A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bz^1)=tA^{L,j_{L}}\bx + (1-t)A^{L,j_{L}}\bw.
\end{equation}

If $j_{L}\geq 2$, by (<ref>) and (<ref>),

\begin{equation}\label{prop1:eq3}
A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bz^1)\leq tA^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bx) + (1-t)A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bw).
\end{equation}

Let us write

\begin{align*}
&\bu\coloneqq t(\bb^{L} + A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bx) +  W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bx))), \\
&\bv\coloneqq (1-t)(\bb^{L} + A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bw) +  W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bw))).
\end{align*}

Combining (<ref>), (<ref>) and (<ref>) with (<ref>), we get

\begin{align*}\Phi_{\Theta,d}(\bz^1)=\sigma^{L}_d\left(\bb^{L} + A^{L,j_{L}}\Phi_{\Theta}^{j_{L}-1}(\bz^1) +  W^{L}C^{L}(\Phi_{\Theta}^{L-1}(\bz^1))\right)&\leq \sigma^{L}_d(t\bu+(1-t)\bv)\\
&\leq t\sigma^{L}_d(\bu) + (1-t)\sigma^{L}_d(\bv)\\
&=t\Phi_{\Theta,d}(\bx) + (1-t)\Phi_{\Theta,d}(\bw).\end{align*}

This concludes the proof.

max totalsize=.8.8,center
[auto,node distance =3 cm ,on grid ,
thick ,
isosceles triangle,
isosceles triangle apex angle=60,
inner sep=1pt,
top color = amber,
bottom color = amber
minimum width=2cm,
minimum height=1cm,
top color = azure,
bottom color = azure
input/.style =draw, circle, top color = darkpastelgreen, bottom color =darkpastelgreen, black, text=black,inner sep=0.5pt,
ghost/.style = circle ,top color =white , bottom color = white ,
draw,white , text=white , minimum width =1 cm,
bias/.style = circle ,top color =white , bottom color = white ,
draw, white , text=black , inner sep=0.5pt,
output/.style =draw, circle, top color = amethyst, bottom color =amethyst, black, text=black,inner sep=0.5pt,
inner sep=2pt,

top color = azure,
bottom color = azure,

[square2, top color=white, bottom color=white] (0) [] $\boldsymbol{x}$;
[square1, shape border rotate=270] (1/2) [below =1.5cm of 0] $\sigma^1$;
[square2] (1) [below =1.5cm of 1/2] $\boldsymbol{z}^2$;
[square1, shape border rotate=270] (3/2) [below =1.5cm of 1] $\sigma^2$;
[square2] (2) [below =1.5cm of 3/2] $\boldsymbol{z}^3$;
[square1, shape border rotate=270] (5/2) [below =1.5cm of 2] $\sigma^3$;
[square2] (3) [below =1.5cm of 5/2] $\boldsymbol{z}^4$;
[square1] (7/2) [right =2.5cm of 3] $\sigma^4$;
[square2] (4) [right =2.5cm of 7/2] $\boldsymbol{z}^5$;
[square1] (9/2) [right =2.5cm of 4] $\sigma^5$;
[square2] (5) [right =2.5cm of 9/2] $\boldsymbol{z}^6$;
[square1, shape border rotate=90] (11/2) [above =1.5cm of 5] $\sigma^6$;
[square2] (6) [above =1.5cm of 11/2] $\boldsymbol{z}^7$;
[square1, shape border rotate=90] (13/2) [above =1.5cm of 6] $\sigma^7$;
[square2] (7) [above =1.5cm of 13/2] $\boldsymbol{z}^8$;
[square1, shape border rotate=90] (15/2) [above =1.5cm of 7] $\sigma^8$;
[square2, top color=darkpastelgreen, bottom color=darkpastelgreen] (8) [above =1.5cm of 15/2] $\boldsymbol{z}$;

(0) edge [black,-,bend right=25] node[left] $A^{2,1}$ (1);
(0) edge [black,-,bend right=45] node[left] $A^{4,1}$ (3);

(0) edge [black,-] node[right] $W^1, \bb^1$ (1/2);
(1/2) edge [black,->] node[right]  (1);
(1) edge [black,-] node[right] $W^2_{\geq0}, \bb^2$ (3/2);
(3/2) edge [black,->] node[right]  (2);
(2) edge [black,-] node[right] $W^3_{\geq0}, \bb^3$ (5/2);
(5/2) edge [black,->] node[right]  (3);

(3) edge [black,-] node[above] $W^4_{\geq0}, \bb^4$ (7/2);
(7/2) edge [black,->] node[below]  (4);
(4) edge [black,-] node[above] $W^5_{\geq0}, \bb^5$ (9/2);
(9/2) edge [black,->] node[below]  (5);

(5) edge [black,-] node[left] $W^6_{\geq0}, \bb^6$ (11/2);
(11/2) edge [black,->] node[left]  (6);
(6) edge [black,-] node[left] $W^7_{\geq0}, \bb^7$ (13/2);
(13/2) edge [black,->] node[left]  (7);
(7) edge [black,-] node[left] $W^8_{\geq0}, \bb^8$ (15/2);
(15/2) edge [black,->] node[left]  (8);

(3) edge [black,-, bend left= 30] node[below] $A^{6,4}_{\geq0}$ (5);

(2) edge [black,double,-] node[below] $C(\boldsymbol{z}^3,\bz^7)$ (6);

at (5,-1) [] (layer) ;

[below=4cm] (S) at (layer.south) Convex Neural Network;

Continuation of the example in Figure <ref>. All the functions $\sigma^k$ are now assumed to be component-wise convex and monotone nondecreasing, and this is represented by the triangular-shaped nodes. Observe that the matrices $W^k$, for $k=2,\ldots,8$, and $A^{6,4}$ have nonnegative entries, while $W^1$, $A^{2,1}$, $A^{4,1}$ and the bias terms $\bb^k$ do not have any constraints.

See Figure <ref> for a visual representation of a component-wise  convex neural network. Now, given a component-wise convex neural network $\Phi^c_{\Theta}\colon X \to Z$ of depth $L$,  define

\begin{equation}\tag{ucNN}\label{eq:ucNN}
\begin{cases}
\Phi^{uc}_\Theta: X \to Z',\\
\Phi^{uc}_\Theta (\bx)\coloneqq \Sigma_{\boldsymbol{f},\boldsymbol{g} }\left(\bx,\Phi^c_{\Theta}(\bx)\right),
\end{cases}
\end{equation}

where $Z'\coloneqq (\R^{D'},\|\cdot\|_{Z'})$ and $\Sigma_{\boldsymbol{f},\boldsymbol{g}} \colon X\times Z \to Z'$ is given by

\begin{equation*}
\Sigma_{\boldsymbol{f},\boldsymbol{g} }(\bx,\bz) \coloneqq (f_1(\bx)+g_1(\bz),\ldots, f_{D'}(\bx)+g_{D'}(\bz))^T,
\end{equation*}

and $f_d: X\to \R$, $g_d:Z \to \R$ are generic real-valued functions for $d=1,\ldots,D'$. In particular, $\Phi_{\Theta}^{uc}$ is  a  neural network $\Phi^{uc}_\Theta$ of depth $L+1$, whose last layer $\sigma^{L+1}(\hat{\bb}^{L+1}+W^{L+1}\hat{\bz}^{L+1})$ is given by

\begin{cases}
\hat{\bz}^{L+1}=C(\bz^1,\bz^{L+1}),\quad	W^{L+1}=I, \quad \bb^{L+1}=A^{L+1,j_{L+1}}=\boldsymbol{0},\\
\sigma^{L+1}=\Sigma_{\boldsymbol{f},\boldsymbol{g}},
\end{cases}

W^{L+1}, \bb^{L+1}, A^{L+1,j_{L+1}} \in \Theta_{frozen}.

max totalsize=.8.8,center
[auto,node distance =3 cm ,on grid ,
thick ,
isosceles triangle,
isosceles triangle apex angle=60,
inner sep=1pt,
top color = amber,
bottom color = amber
isosceles triangle,
isosceles triangle apex angle=60,
inner sep=1pt,
top color = Red2,
bottom color = Red2
minimum width=2cm,
minimum height=1cm,
top color = azure,
bottom color = azure
minimum width=2cm,
minimum height=1cm,
top color = Red2,
bottom color = Red2
input/.style =draw, circle, top color = darkpastelgreen, bottom color =darkpastelgreen, black, text=black,inner sep=0.5pt,
ghost/.style = circle ,top color =white , bottom color = white ,
draw,white , text=white , minimum width =1 cm,
bias/.style = circle ,top color =white , bottom color = white ,
draw, white , text=black , inner sep=0.5pt,
output/.style =draw, circle, top color = amethyst, bottom color =amethyst, black, text=black,inner sep=0.5pt,
inner sep=2pt,

top color = azure,
bottom color = azure,

[square2, top color=white, bottom color=white] (0) [] $\boldsymbol{x}$;
[square1, shape border rotate=270] (1/2) [below =1.5cm of 0] $\sigma^1$;
[square2] (1) [below =1.5cm of 1/2] $\boldsymbol{z}^2$;
[square1, shape border rotate=270] (3/2) [below =1.5cm of 1] $\sigma^2$;
[square2] (2) [below =1.5cm of 3/2] $\boldsymbol{z}^3$;
[square1, shape border rotate=270] (5/2) [below =1.5cm of 2] $\sigma^3$;
[square2] (3) [below =1.5cm of 5/2] $\boldsymbol{z}^4$;
[square1] (7/2) [right =2.5cm of 3] $\sigma^4$;
[square2] (4) [right =2.5cm of 7/2] $\boldsymbol{z}^5$;
[square1] (9/2) [right =2.5cm of 4] $\sigma^5$;
[square2] (5) [right =2.5cm of 9/2] $\boldsymbol{z}^6$;
[square1, shape border rotate=90] (11/2) [above =1.5cm of 5] $\sigma^6$;
[square2] (6) [above =1.5cm of 11/2] $\boldsymbol{z}^7$;
[square1, shape border rotate=90] (13/2) [above =1.5cm of 6] $\sigma^7$;
[square2] (7) [above =1.5cm of 13/2] $\boldsymbol{z}^8$;
[square1, shape border rotate=90] (15/2) [above =1.5cm of 7] $\sigma^8$;
[square2, top color=darkpastelgreen, bottom color=darkpastelgreen] (8) [above =1.5cm of 15/2] $\boldsymbol{z}$;
[square1uc, shape border rotate=90] (17/2) [above =1.4cm of 8] $\Sigma_{\boldsymbol{f},\boldsymbol{g}}$;
[square2uc] (9) [above =1.5cm of 17/2] $\boldsymbol{z}'$;

(0) edge [black,-,bend right=25] node[left] $A^{2,1}$ (1);
(0) edge [black,-,bend right=45] node[left] $A^{4,1}$ (3);
(0) edge [Red2,double,-] node[above] $C(\bx,\bz)$ (8);

(0) edge [black,-] node[right] $W^1, \boldsymbol{b}^1$ (1/2);
(1/2) edge [black,->] node[right]  (1);
(1) edge [black,-] node[right] $W^2_{\geq0}, \boldsymbol{b}^2$ (3/2);
(3/2) edge [black,->] node[right]  (2);
(2) edge [black,-] node[right] $W^3_{\geq0}, \boldsymbol{b}^3$ (5/2);
(5/2) edge [black,->] node[right]  (3);

(3) edge [black,-] node[above] $W^4_{\geq0}, \boldsymbol{b}^4$ (7/2);
(7/2) edge [black,->] node[below]  (4);
(4) edge [black,-] node[above] $W^5_{\geq0}, \boldsymbol{b}^5$ (9/2);
(9/2) edge [black,->] node[below]  (5);

(5) edge [black,-] node[left] $W^6_{\geq0}, \boldsymbol{b}^6$ (11/2);
(11/2) edge [black,->] node[left]  (6);
(6) edge [black,-] node[left] $W^7_{\geq0}, \boldsymbol{b}^7$ (13/2);
(13/2) edge [black,->] node[left]  (7);
(7) edge [black,-] node[left] $W^8_{\geq0}, \boldsymbol{b}^8$ (15/2);
(15/2) edge [black,->] node[left]  (8);
(8) edge [black,-] node[left] $I$ (17/2);
(17/2) edge [black,->] node[left]  (9);

(3) edge [black,-, bend left= 30] node[below] $A^{6,4}_{\geq0}$ (5);

(2) edge [black,double,-] node[below] $C(\boldsymbol{z}^3,\bz^7)$ (6);

at (5,-1) [] (layer) ;

[below=4cm] (S) at (layer.south) Uniformly Convex Neural Network;

Continuation of the examples in Figures <ref> and <ref>. This architecture depicts a component-wise uniformly convex neural network $\Phi^{uc}_\Theta$ obtained from the component-wise convex neural network $\Phi^{c}_\Theta$ in Figure <ref>. The differences are depicted in red. Following equation (<ref>), it has been added a  concatenation operator $C(\bx,\bz)$ and a new layer after $\bz=\Phi^{c}_\Theta(\bx)$.

Consider the following assumptions, for every $d=1,\ldots,D'$:

\begin{align}
&f_d \quad  \mbox{continuous and uniformly convex};\tag{C3}\label{hypothesis4}\\
&g_d\quad  \mbox{continuous, convex and monotone nondecreasing on the image } \Phi^c_{\Theta}(X)\subseteq Z.\tag{C4}\label{hypothesis5}
\end{align}

Under the constraints (<ref>) and (<ref>), the neural network $\Phi^{uc}_\Theta$ defined in (<ref>) is component-wise uniformly convex.

Given a component-wise convex map $\Phi^c_{\Theta} : X \to Z$, by direct computation it is immediate to check that $f_d + g_d\circ\Phi^c_{\Theta}$ is uniformly convex for every $d=1,\ldots, D'$, under the constraints (<ref>) and (<ref>).

Combining Theorem <ref> and Theorem <ref>, given a neural network $\Phi_{\Theta} : X \to Z$, it is now easy to build a continuous and component-wise uniformly convex neural network $\Phi^{uc}_{\Theta} : X \to Z'$ by employing the following steps:

*  Fix a neural network $\Phi_{\Theta} : X \to Z$ whose architecture is given by (<ref>).

*  Impose the constraints (<ref>)-(<ref>) to get  $\Phi^c_{\Theta} : X \to Z$.

*  Define $\Phi^{uc}_{\Theta} : X \to Z'$ by equation (<ref>) and impose the constraints (<ref>)-(<ref>).

See Figure <ref> for a visual representation of a component-wise uniformly convex neural network.

§ EXAMPLE: CONVEX AND UNIFORMLY CONVEX U-NET WITH BATCH NORMALIZATION

Here we provide particular examples of convex and uniformly convex U-net, which we will use in Subsection <ref> and Section <ref>. U-net is a convolutional neural network originally designed for biomedical image segmentation <cit.>. The network consists of a contracting path and a symmetric expanding path so that it has a U-shaped architecture. Several variants of the U-net architecture have been efficiently implemented in the solutions of inverse problems, e.g. <cit.>.

We propose a modified U-net architecture as shown in Figure <ref>, and we will explain the approaches to obtain convex and uniformly convex U-net architectures according to the general theories proposed in section <ref>.
In Subsection <ref>, we give a brief description of the batch normalization which is one of the extra features we add to the modified U-net. In Subsection <ref>, we provide a mathematical modeling of the modified U-net architecture using the notations of general neural networks introduced in  Subsection <ref>. Then in Subsections <ref> and <ref>, we explain the approaches to achieve convex and uniformly convex U-net architectures.

A modified U-net architecture.

§.§ Batch Normalization (BN)

Batch normalization (BN) is a state-of-the-art technique for accelerating the training speed and improving the performance of optimization <cit.>. The batch normalization (BN) works differently in the training process and in the inference process after training. Hereby we give a brief description of the BN operation according to our notations for general neural networks in equation (<ref>).

With a given batch size $n_0$ and a mini-batch of the initial input,
consider the output after convolution at the $k$-th layer of the neural network $\Phi_\Theta$,

\mathcal{B}^k= \{ \bu^k_1,\ldots, \bu^k_{n_0}\}, \quad \mathrm{where}\quad \bu^k_j\coloneqq \hat{\bb}^k+ W^k\hat{\bz}^k_j,\quad j=1,\cdots,n_0\,.

The operation of $\operatorname{BN}$ reads as follows,

\begin{equation*}
\operatorname{BN}(\bu)\coloneqq \boldsymbol{\beta}^k + \operatorname{diag}(\boldsymbol{\gamma}^k)  \frac{\bu - \mathds{E}(\mathcal{B}^k)}{\sqrt{\operatorname{Var}(\mathcal{B}^k)+\epsilon}},\quad \forall\, \bu\in\mathcal{B}^k
\end{equation*}

where $\boldsymbol{\gamma}^k$ and $\boldsymbol{\beta}^k$ are column vectors and the free parameters to be learned, $\mathds{E}(\mathcal{B}^k)$ and $\operatorname{Var}(\mathcal{B}^k)$ denote mean value and variance of the elements in $\mathcal{B}^k$, respectively, and  $\epsilon>0$ is a small constant introduced for numerical stability.

In the training process, since $\mathds{E}(\mathcal{B}^k)$ and $\operatorname{Var}(\mathcal{B}^k)$ will involve the variable $\bu\in\mathcal{B}^k$,
BN is not necessarily component-wise convex nor monotone nondecreasing. After the training process is completed, $\mathds{E}(\mathcal{B}^k)$ and $\operatorname{Var}(\mathcal{B}^k)$ will be fixed using the full population rather than the mini-batch, and so BN reduces to an affine transformation.

§.§ The modified U-net architecture

We provide a mathematical modeling of the modified U-net architecture $\Phi_{\Theta}$ as shown in Figure <ref>. Here we will make use of the notations of general neural networks introduced in Subsection <ref>. Before proceeding further, it is necessary to briefly introduce the following operators.

$\blacksquare$ Convolution: It is a linear operator which can be represented by a tensor $K$ of dimension $f_1\times f_2\times c_1\times c_2$, where $f_1\times f_2$ is the filter size of the convolution, $c_1$ and $c_2$ denote the number of channels of the input and output images, respectively. More details about the convolution operator can be found in <cit.> and <cit.>. In the U-net architecture, a bias term is added after every convolution.

$\blacksquare$ Zero padding: It is a linear operator which can be achieved by embedding a given matrix of dimension $n\times m$ into a zero matrix of dimension $n'\times m'$, with $n'\geq n$ and $m'\geq m$. By reshaping the matrices before and after zero padding to be column vectors, zero padding can be expressed as a matrix product, $\mathbf{M}=P_0\mathbf{N}$, where $P_0$ is a binary matrix with entries from $\{0,1\}$.

$\blacksquare$ Rectified linear unit (ReLU): $\bx\mapsto \max\{0,\bx\}$, where the maximum is applied component-wise.

$\blacksquare$$\operatorname{max-pool}$: Given a matrix of real numbers, partition it into $i\times j$ sub-matrices; the $\operatorname{max-pool}_{i\times j}(\cdot)$ returns a matrix that calculates the maximum value for each sub-matrix in the $i\times j$ partition.

$\blacksquare$ Batch normalization (BN): Details are given in Subsection <ref>. We emphasis that BN reduces to an affine transformation after training.

In the modified U-net architecture $\Phi_{\Theta}$, every convolution operator is combined with a zero padding, such that the input and output images have the same size. In addition, we add a Batch normalization (BN) between each convolution and ReLU function. As shown in Figure <ref>, the modified U-net architecture can be constructed by several unit blocks: right arrow (involving concatenation), down arrow, and up arrow. Accordingly, we provide mathematical modeling for these unit blocks.

* Right arrow [baseline=(a.base)](a) at (0pt,0pt) [] ;
[-Triangle[width=8pt,length=4pt], color=arrow, line width=5pt]($(a)+(0pt,4pt)$) – ($(a)+(12pt,4pt)$); : A $k$-th right arrow reads as

\operatorname{ReLU}(\operatorname{BN}(\bb^k+K^kP_0\hat{\bz}^k)).

According to the notations in equation (<ref>), the block of right arrow can be modeled as the composition of two layers,

\begin{equation} \label{eqn_RA}
\begin{cases}
\bz^{k+1}=\sigma^{k}(\hat{\bb}^{k} + W^{k}\hat{\bz}^{k}),\ \hat{\bz}^{k+1}=\bz^{k+1} \\
\bz^{k+2}=\sigma^{k+1}(\hat{\bb}^{k+1} + W^{k+1}\hat{\bz}^{k+1})
\end{cases}\,,
\end{equation}

\[
\begin{cases}W^k=K^kP_0,\quad \hat{\bb}^k=\bb^k \\
\sigma^k(\bu)= \frac{\bu - \mathds{E}(\mathcal{B}^k)}{\sqrt{\operatorname{Var}(\mathcal{B}^k)+\epsilon}} \\
K^k,\,\bb^k\in \Theta_{free}, \quad P_0\in \Theta_{frozen}
\end{cases}\,, \quad \mathrm{and} \quad
\begin{cases}W^{k+1}=\operatorname{diag}(\boldsymbol{\gamma}^{k}),\quad \hat{\bb}^{k+1}=\boldsymbol{\beta}^{k} \\
\sigma^{k+1}(\bu)= \operatorname{ReLU}(\bu) \\
\boldsymbol{\gamma}^{k},\,\boldsymbol{\beta}^{k}\in \Theta_{free}
\end{cases}\,.
\]

$P_0$ denotes the matrix of zero padding, which is a binary matrix with entries from $\{0, 1\}$. The filter size of every convolution $K^k$ is $3\times 3$ and the stride is $1$. In the ascending part of the U-net architecture, the concatenation operation is involved before the first convolution of each layer, and it can be modeled by the formulation of general neural networks as shown in Subsection <ref>, i.e.  $\hat{\bz}^{k}=C(\bz^{i_k},\bz^k)$.

* Down arrow [baseline=(a.base)](a) at (0pt,0pt) [] ;
[-Triangle[width=8pt,length=4pt], color=arrow, line width=5pt]($(a)+(0pt,12pt)$) – ($(a)+(0pt,0pt)$); : An $i$-th down arrow reads as

It can be simply modeled by the formula of general neural networks as shown in equation (<ref>),

\begin{equation}\label{eqn_DA}
\sigma^i(\hat{\bb}^i+W^i\hat{\bz}^i)\,,
\end{equation}

where $ \hat{\bb}^i=\boldsymbol{0}$, $W^i=I$, and $\sigma^i(\bu)= \operatorname{max-pool}_{2\times2}(\bu)$; $\sigma^i$ is continuous, component-wise convex and monotone nondecreasing.

* Up arrow [baseline=(a.base)](a) at (0pt,0pt) [] ;
[-Triangle[width=8pt,length=4pt], color=arrow, line width=5pt]($(a)+(0pt,0pt)$) – ($(a)+(0pt,12pt)$); : The up arrow has the same structure of the right arrow, and it can be modeled in the same way by equation (<ref>). There are some minor differences in the convolution $K^k$: (i) the filter size of $K^k$ is $2\times 2$; (ii) the way of zero padding is different from that in the right arrow, but since it is still a zero padding, the matrix $P_0$ is still a binary matrix with entries from $\{0,1\}$.

* Last right arrow [baseline=(a.base)](a) at (0pt,0pt) [] ;
[-Triangle[width=8pt,length=4pt], color=last-arrow, line width=5pt]($(a)+(0pt,4pt)$) – ($(a)+(12pt,4pt)$); :  The last right arrow is just a convolution with filter size $1\times 1$. It can be modeled as part of equation (<ref>):

\begin{equation} \label{eqn_last_RA}
\bz^{L+1}=\sigma^{L}(\hat{\bb}^{L} + W^{L}\hat{\bz}^{L})\,,
\end{equation}

where $W^L=K^L P_0$, $\hat{\bb}^L=\bb^L$, and  $\sigma^L(\bu)=\bu$. The output $\bz^{L+1}$ is of size $256\times256\times1$, where the channel number is $1$. As a result, the bias term $\hat{\bb}^L=\bb^L$ is essentially a scalar, i.e. the components of $\bb^L$ are all equal to $b^L_0\in\R$.

§.§ Convex U-net

In Subsection <ref>, we showed that the modified U-net architecture $\Phi_{\Theta}$ can be modeled by the set of equations (<ref>). Then according to Theorem <ref>, we can achieve a convex U-net by imposing the constraints (<ref>) and (<ref>). Since the skip connection described by $A^{k,j_k}$ is not involved in $\Phi_{\Theta}$, we only need to consider the properties of $\sigma^k$ and $W^k$.

The architecture modeled by equation (<ref>) automatically satisfies the constraints (<ref>) and (<ref>). For the architecture modeled by equation (<ref>) (which includes the formulation in equation (<ref>) as well), $\sigma^k$ and $\sigma^{k+1}$ are continuous, convex, and monotone nondecreasing after training, so that the constraint (<ref>) is satisfied. To impose the constraint (<ref>) on $W^k$ and $W^{k+1}$, we have to impose nonnegativity constraint on the entries of $K^k$ and $\boldsymbol{\gamma}^{k}$. Above all, we have the following conclusion.

By imposing nonnegativity constraint on $\{K^k\mid\forall k\ge2\}$ and $\{\boldsymbol{\gamma}^{k}\mid\forall k\ge1\}$, the modified U-net $\Phi_{\Theta}$ turns out to be a convex U-net in the inference process after training.

In next subsection, we will build a uniformly convex U-net by employing the architecture of convex U-net. We will ask the output of the convex U-net to be nonnegative, which is helpful to construct the uniformly convex U-net in a very simple way. Considering the convex U-net obtained in Proposition <ref>, thanks to the ReLU function applied just before the last convolution operator, and the nonnegativity of the entries of the last convolution operator itself, the output $\bz^{L+1}$ must be nonnegative if we further impose nonnegativity constraint on the last bias term $\bb^L$.

By imposing nonnegativity constraint on $\{K^k\mid\forall k\ge2\}$, $\{\boldsymbol{\gamma}^{k}\mid\forall k\ge1\}$ and $\bb^L$, the modified U-net $\Phi_{\Theta}$ turns out to be a convex U-net $\Phi_{\Theta}^c$ having nonnegative output, i.e. the components of $\bz^{L+1}=\Phi_{\Theta}^c(\bx)$ are all nonnegative.

§.§ Uniformly convex U-net

With the trained convex U-net $\Phi_{\Theta}^c$ having nonnegative output, one can construct a uniformly convex U-net according to the following formula,

\begin{equation}\label{eq:Reg}
\Phi_{\Theta}^{uc}(\bx)\coloneqq a\|\bx\|_p^p+\| \Phi^{c}_\Theta(\bx) \|_q^q\,,\qquad p\ge2,\ q\ge1\,,
\end{equation}

where $\|\cdot\|_p$ and $\|\cdot\|_q$ denote the standard $l^p$-norm and $l^q$-norm, respectively, and $a>0$ is a small positive constant.

Indeed, formula (<ref>) is in the form of equation (<ref>), with

\[
D'=1,\quad f_1(\bx)=a\|\bx\|_p^p,\quad\mathrm{and}\quad g_1(\bz)=\|\bz\|_q^q\,.
\]

The function $f_1(\bx)=a\|\bx\|_p^p$ is uniformly convex for $p\geq 2$; e.g., see <cit.> and <cit.>. The function $g_1(\bz)=\|\bz\|_q^q\ (q\ge1)$ is continuous, convex and monotone nondecreasing in the domain $\R^N_{\geq0}:=\{\bz\in\R^N\mid z_i\ge0,\ i=1,\cdots,N\}$. Moreover, as shown in Proposition <ref>, the output of the convex U-net $\Phi_{\Theta}^c$ satisfies $\Phi_{\Theta}^c\in\R^N_{\geq0}$. As a result, $\Phi_{\Theta}^{uc}$ in equation (<ref>) is uniformly convex according to Theorem <ref>. Equation (<ref>) gives a very simple way to construct the uniformly convex U-net.

§ ITERATED NETWORK TIKHONOV (INETT) METHOD

In this section we propose the iterated network Tikhonov (iNETT) method. First, fix a continuous, uniformly  convex and well trained neural  network

\begin{equation*}
\Phi^{uc}_\Theta : X \to (\R, |\cdot|),
\end{equation*}

for example by applying the steps <ref>, <ref> and <ref> at the end of Section <ref>, and define

$$\Reg\coloneqq \Phi_{\Theta}^{uc}.

For the numerical experiments in Section <ref>, we will adopt the modified U-net architecture presented in Section <ref> as equation (<ref>).

Fix $r\in(1,\infty)$, $\bx_0\in X$, $\xi_0 \in \partial \Reg(\bx_0)$  and a sequence of positive real numbers $\{\alpha_n\}_n$ such that

\sum_{n=1}^\infty \alpha_n^{-1}=\infty \quad \mbox{and} \quad \alpha_n\leq c\alpha_{n+1}, \quad \mbox{where } c>0 \mbox{ is a constant}.

The iNETT method is then given by

\begin{equation}\tag*{iNETT}\label{iterNETT2}
\begin{cases}
\bx^\delta_n:= \underset{\bx \in X}{\argmin} \frac{1}{r}\|F \bx - \by^\delta\|_Y^r + \alpha_n\Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1}),\\
\boldsymbol{\xi}^\delta_{n}:= \boldsymbol{\xi}^\delta_{n-1} - \frac{1}{\alpha_{n}} F^TJ_r\left(F\bx_n^\delta - \by^\delta\right),\\
\bx_0 \in X, \,	\boldsymbol{\xi}_0 \in \partial \Reg(\bx_0).
\end{cases}
\end{equation}

The stopping rule associated to <ref> is the standard discrepancy principle. Specifically, we stop the iterations at the first step $n_\delta=n$ such that

\begin{equation}\label{dicrepancy_principle}
\left\|F\bx_{n}^\delta  - \by^\delta\right\|_Y\leq \tau\delta< \left\|F\bx_{n-1}^\delta  - \by^\delta\right\|_Y,
\end{equation}

where $\tau>1$ is a fixed constant, and $\delta$ is the parameter which estimates the level of noise $\boldsymbol{\eta}$, viz. $\|\boldsymbol{\eta}\|_Y=\|\by^\delta-\by\|_Y\leq \delta$. Following we will provide convergence analysis and implementation details of the iNETT method.

§.§ Convergence analysis

We present the proof of well-posedness, convergence and stability  of the <ref> method. The proof is mainly a straightforward adaptation, in this finite dimensional setting, of  standard results of general convex regularization theory (e.g. <cit.>) and <cit.>. For the convenience of the reader and to make the theoretical treatment self-contained, we develop here the main points.

Let $F : X \to Y$ be a linear operator between finite dimensional normed spaces and suppose that  (<ref>) and (<ref>) are valid.  Let $\Reg\coloneqq \Phi_\Theta^{uc}$ be a continuous and uniformly convex neural network. Then the method <ref> is well-posed and there exists a unique $\Breg_{\bxi_0}$-minimizing solution $\bx^\dagger$. The method <ref>, coupled with the stopping rule (<ref>), stops in finite steps $n_\delta < \infty$ and it converges to  $\bx^\dagger$ as $\delta\to 0$. In particular,

\begin{equation*}
\bx_{n_\delta}^\delta \to \bx^\dagger, \quad \Reg(\bx_{n_\delta}^\delta)\to \Reg(\bx^\dagger), \quad \Breg_{\boldsymbol{\xi}^\delta_{n_\delta}}(\bx^\dagger,\bx_{n_\delta}^\delta)\to 0
\end{equation*}

as $\delta \to 0$.

About the existence of $\bx^\dagger$ in Definition <ref>, observe that

\begin{equation*}
0\leq c\coloneqq \inf \left\{ \Breg_{\bxi_0}(\bx,\bx_0) \mid x\in X,\; Fx = y  \right\}
\end{equation*}

is well-defined because of (<ref>).  Since $\Breg_{\boldsymbol{\xi}_0}(\bx,\bx_0)$ is continuous and coercive (see Remark <ref>), the level sets $\{\bx \mid \Breg_{\bxi_0}(\bx,\bx_0)\leq M\}$ are compact for any $M>0$, and therefore there exists a minimizer $\bx^\dagger$, by standard topological arguments, which is unique because of the linearity of $F$. In the same fashion, and because of the convexity of $\Breg_{\boldsymbol{\xi}^\delta_{n-1}}$, $\bx_n^\delta$ is well-defined and unique at each step.

Observe now that $J_r$ is a single-valued map thanks to (<ref>), and that $\bxi_n^\delta \in \partial \Reg(\bx^\delta_n)$ in virtue of the minimality of $\bx_n^\delta$. Indeed, if we define

\mathcal{G}(\bx)\coloneqq  \frac{1}{r}\|F \bx - \by^\delta\|_Y^r + \alpha_n\Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1}),

by subgradient calculus (see for example <cit.>), and because $J_r$ is the subgradient of the map $\by \mapsto \frac{1}{r}\|\by\|_Y^r$, it holds that

\begin{equation*}
0\in\partial \mathcal{G}(\boldsymbol{\bx}^\delta_{n}) = F^TJ_r(F\bx^\delta_{n}-\by^\delta) + \alpha_n\partial\Reg(\bx^\delta_{n}) -\alpha_n\bxi^\delta_{n-1}.
\end{equation*}

The rest of the proof is a straightforward application of <cit.>.

§.§ Implementation

Since our applications in Section <ref> will be based on computerized tomography, from now on we will consider $X$ and $Y$ as vector spaces of digital images. We will use the standard $\ell^2$-norm for $X$, that is, $X=(\R^N,\|\cdot\|_2)$ and a normalized $\ell^2$-norm for $Y= \left(\R^M,\|\cdot\|_Y\right)$ as $\|\by\|_Y\coloneqq\frac{1}{\sqrt{M}}\|\by\|_2$, and we will fix $r=2$. Therefore, the  <ref> iteration  will assume the following form (see also Remark <ref>),

\begin{equation}\tag*{iNETT}\label{iterNETT3}
\begin{cases}
\bx^\delta_n= \underset{\bx \in X}{\argmin}\frac{1}{2M}\|F \bx - \by^\delta\|_2^2+ \alpha_n\Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1})\,,\\
\boldsymbol{\xi}^\delta_{n}= \boldsymbol{\xi}^\delta_{n-1} - \frac{1}{\alpha_{n}}\frac{1}{M}  F^T\left(F\bx_n^\delta - \by^\delta\right)\,, \\
\bx_0 \in X, \,	\boldsymbol{\xi}_0 \in \partial \Reg(\bx_0)\,.
\end{cases}
\end{equation}

The minimum problem in <ref> can be solved by a gradient descent approach. Recalling that the Bregman distance of the regularization term is defined as

\[
\Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1})=\Reg(\bx)-\Reg(\bx^\delta_{n-1})-\langle \boldsymbol{\xi}^\delta_{n-1}, \bx-\bx^\delta_{n-1}\rangle\,,
\]

the gradient descent algorithm for $\bx^\delta_n$ has the following form,

\begin{equation*} \label{eq:gradient_descent}
\bx^\delta_{n,k+1}=\bx^\delta_{n,k}-s\cdot\left[ \frac{1}{M}F^T\left(F\bx_{n,k}^\delta - \by^\delta\right)+\alpha_{n}\left(\frac{\partial\Reg(\bx)}{\partial\bx}{\Big|}_{\bx=\bx^\delta_{n,k}}-\boldsymbol{\xi}^\delta_{n-1} \right) \right]\,, \quad k\in\mathbb{N}
\end{equation*}

where $s>0$ denotes the step size of the gradient descent iteration. The initial guess for the iteration can be taken as $\bx^\delta_{n,0}=\bx^\delta_{n-1}$.

§.§ Building a uniformly convex neural network regularizer

We will adopt the uniformly convex neural network (<ref>) presented in Subsection <ref>, with $p=q=2$. That is,

\begin{equation*}
\Reg(\bx)\coloneqq \Phi_{\Theta}^{uc}(\bx)=a\|\bx\|_2^2+\| \Phi^{c}_\Theta(\bx) \|_2^2,
\end{equation*}

where $\Phi^{c}_\Theta(\bx)$ is a well trained convex U-net having nonnegative output as shown in Proposition <ref>, and $a>0$ is a small positive constant.

Since the regularizer is designed to penalize artifacts in the solution, the trained U-net $\Phi^{c}_\Theta(\bx)$ should have small output for artifact-free images and have large output for images with artifacts. We follow the strategy proposed in  <cit.> for the preparation of training data. Let $\{\bx_s^*\mid s=1,\cdots,N_1+N_2\}$ denote a set of images in $X$. The training images are constructed in the following way,

\begin{equation}\label{train_input}
\boldsymbol{z}_s=\left\{
\begin{array}{lcr}
F^\dagger (F\bx_s^*+\boldsymbol{\eta}) &,& s=1,\cdots,N_1 \\
\bx_s^*&,& s=N_1+1,\cdots,N_1+N_2
\end{array}
\right.
\end{equation}

where $\boldsymbol{\eta}$ denotes noise perturbations on the data, and $F^{\dagger}$ denotes the pseudo inverse of the forward operator $F$. In the set of training images $\{\boldsymbol{z}_s\mid s=1,\cdots,N_1+N_2\}$, the first $N_1$ $\boldsymbol{z}_s$'s simulate the images with artifacts, while the last $N_2$ $\boldsymbol{z}_s$'s simulate the artifact-free images. The output label is the error between $\boldsymbol{z}_s$ and $\bx_s^*$,

\begin{equation}\label{train_output}
\boldsymbol{r}_s=|\bx_s^*-\boldsymbol{z}_s|,\quad s=1,\cdots,N_1+N_2\,,
\end{equation}

where the absolute value is applied component-wise, and the training set is composed of the pairs of input images and output labels:
$$\{(\boldsymbol{z}_s,\boldsymbol{r}_s)\mid s=1,\cdots,N_1+N_2\}.$$

The loss function of one training sample is defined as follows,

\[
L_s(\Theta)=\|\Phi^c_{\Theta}(\boldsymbol{z}_s)- \boldsymbol{r}_s \|_2^2+\lambda\|\Theta_{free}\|_2^2\,,
\]

\[
\Theta_{free}= \{\boldsymbol{\beta}^{k};\boldsymbol{\gamma}^{k};\bb^k;K^k\}_{k}
\]

considering the modified U-net architecture as shown in Section <ref>. We propose to use a mini-batch optimization approach in the training process, so that the cost function for updating $\Theta_{free}$ has the following form,

\begin{equation*}\label{eq:cost-function}
\mathcal{L}_t(\Theta)=\frac{1}{|\mathcal{I}_t|}\sum_{s\in\mathcal{I}_t} L_s(\Theta) = \frac{1}{|\mathcal{I}_t|}\sum_{s\in\mathcal{I}_t} \|\Phi^c_{\Theta}(\boldsymbol{z}_s)- \boldsymbol{r}_s \|_2^2+\lambda\|\Theta_{free}\|_2^2\,,
\end{equation*}

where $\mathcal{I}_t\subset\{1,2,\cdots,N_1+N_2\}$ denotes a mini-batch of the index set, and $|\mathcal{I}_t|$ denotes the batch size. The parameter $\lambda$ controls the amount of regularization applied. We utilize a random shuffling strategy <cit.> for the sampling of the mini-batch set $\mathcal{I}_t$. Then the set of training parameters $\Theta_{free}$ is updated according to the Adam algorithm <cit.>. After each updating, we impose non-negativity constraint on $\boldsymbol{\gamma}^{k}$ and $K^k$ for the convexity of the U-net $\Phi^c_{\Theta}$, with the only exception of the convolution matrix $K^1$ in the first layer. In addition, to make the convex U-net $\Phi^c_{\Theta}$ have nonnegative output, we impose non-negativity constraint on the bias term of the last layer $\bb^L$. After the training process is completed, we fix the elements $\mathds{E}(\mathcal{B}^k)$ and $\operatorname{Var}(\mathcal{B}^k)$ which define the batch normalization operators. Algorithm <ref> summarizes the details of the training process.

[!h]Training the convex U-net $\Phi^c_{\Theta}$ having nonnegative output
[1]
Initialize the  free parameters set $\Theta_{free}=\{\boldsymbol{\beta}^{k};\boldsymbol{\gamma}^{k};\bb^k;K^k\}_{k}$.
epoch $e=1,2,\cdots,N_e$
Sample a random permutation $\{\pi(1), \pi(2),\ldots, \pi(N_1+N_2)\}$ of the index set $\{1,2,\ldots,N_1+N_2\}$.
With a given batch size $n_0$, partition the shuffled index set into $n=\frac{N_1+N_2}{n_0}$ subsets: $\mathcal{I}_t=\{\pi((t-1)n_0+1),\ldots,\pi(t\,n_0)\}$, $t=1,2,\cdots,n$.
$t=1$ to $n$
Update $\Theta_{free}$ according to the Adam algorithm and mini-batch cost function  $\mathcal{L}_t(\Theta)$.
Impose nonnegativity constraint: $K^k=\max\{0,K^k\}$, $\forall k\ge2$;  $\boldsymbol{\gamma}^{k}=\max\{0,\boldsymbol{\gamma}^{k}\}$,

$\forall k\ge1$;  $\bb^L=\max\{0,\bb^L\}$, where $\bb^L$ denotes the bias term of the last layer.
Fix $\mathds{E}(\mathcal{B}^k)$ and $\operatorname{Var}(\mathcal{B}^k)$ in batch normalization (BN) using the full population rather than the mini-batch.

§.§ Some remarks

As final comments for this section, we provide discussions on two subjects: (1) similarities and differences between <ref> and <ref>; (2) alternative training strategies for the <ref> method.

§.§.§ NETT and iNETT

The <ref> method in <cit.> is developed in the setting of infinite dimensional Banach spaces. The <ref> method can be generalized to the infinite dimensional case as well, where one may require the Banach space X to be reflexive, and re-define the neural network architecture (<ref>) in a proper way. In this paper, we propose <ref> in a finite dimensional case mainly for personal taste. Since we focus on the inverse problems in imaging science, e.g. computerized tomography, the setup is naturally of finite dimension. The images are defined in a finite dimensional vector space, and the neural networks are trained and deployed in the set of images with fixed finite pixels. When considering the accuracy of discretization, it will be necessary to study <ref> in the infinite dimensional setup. A main difficulty will be the convergence analysis: It is unclear under what hypotheses a convex neural network, coupled with training and optimization rules, converges to an infinite dimensional operator.

A major difference between <ref> and <ref> lies in the assumption of coercivity or convexity on the data-driven regularizer. For the convergence analysis of <ref> method, the regularization term $\Reg(\bx)\coloneqq \psi(\Phi_\Theta(\bx))$ is assumed to be coercive <cit.>. On the other hand, in our <ref> method, the neural network regularizer $\Reg(\bx)\coloneqq\Phi_{\Theta}^{uc}(\bx)$ is required to be uniformly convex. In fact, at each iteration of  <ref>, the penalty term $\mathcal{P}_n(\bx)\coloneqq \Breg_{\boldsymbol{\xi}^\delta_{n-1}}(\bx,\bx^\delta_{n-1})$ is coercive, as a consequence of the uniform convexity of $\Reg$ (see Remark <ref>). The assumption of convexity on $\Reg$ is from the iterated Tikhonov method. For example, considering the proofs of <cit.>, a crucial role to prove the (strong) convergence is played by the following inequality,

\Breg_{\boldsymbol{\xi}}(\bx,\hat{\bx})\geq f(\|\bx -\hat{\bx}\|_X),

where $f \colon [0,+\infty) \to [0,+\infty)$ is a strictly increasing continuous function such that $f(0)=0$. The above inequality is true if and only if $\Reg$ is uniformly convex (see <cit.>). It would be interesting to understand if the assumption of uniform convexity on $\Reg$ can be relaxed, say for example just asking coercivity.  To the best of our knowledge, it is yet an open problem.

§.§.§ Alternative training strategies for <ref>

As described in section <ref>, the training strategy of <ref> is adapted from that of <ref> <cit.>. This is a natural choice when we consider the evolution of <ref> into an iterated Tikhonov method. In this way, we can have a simple comparison between the reconstruction results of <ref> and <ref>, which will be provided in the next section.

Nevertheless, it would be useful and interesting to study alternative training strategies for <ref>. In <cit.>, it proposes an adversarial training strategy for the neural network regularizer. The adversarial regularizer learns to discriminate between the distribution of ground truth images and the distribution of unregularized reconstructions. The training can be performed even if only unsupervised training data are available <cit.>. In <cit.>, the adversarial strategy is adopted to successfully build a learned convex regularizer.
This strategy deserves elaborate study in the future work, and it provides a direction to further improve the <ref> algorithm.

Moreover, in <cit.>, it develops a deep ADMM-Net for compressive sensing MRI, which is one of the earliest learned reconstruction methods in imaging science. The ADMM-net maps the ADMM iterative procedure to a data flow graph, and then generalizes the operations of ADMM to have learnable parameters as network layers. This strategy inspires us to generalize <ref> to a fully learned architecture in the future study, e.g. learn the neural network regularizer $\Reg(\bx)$ in the optimization process, which would provide more flexibility of building the data-driven regularizer, even though the convergence analysis can become much more difficult.

§ APPLICATION TO COMPUTERIZED TOMOGRAPHY

We provide applications of the <ref> algorithm, as presented in Subsection <ref>, in 2D computerized tomography (CT). The imaging problem can be modeled as an inverse problem of the Radon transform:

\begin{equation}\label{app1}
y(s,\phi)=\int_{\mathbb R} x(s\cos\phi-t\sin\phi, s\sin\phi+t\cos\phi)\,\mathrm{d}t\,,\quad  (s,\phi)\in\mathbb{R}\times[0,2\pi]\,,
\end{equation}

where $x(\cdot,\cdot)$ denotes the attenuation image to be recovered, and $y$ denotes the projection data. After discretization, equation (<ref>) reduces to a linear system,

\begin{equation}\label{eqn8.2}
\by=F\bx,
\end{equation}

where $\bx=(x_1,\ldots,x_N)^T$ is the vector of the attenuation image, $\by=(y_1,\ldots,y_M)^T$ simulates the measurement projection data, and $F=(a_{ij})_{M\times N}$ denotes the projection matrix; the component $a_{ij}$ represents the contribution of the $j$-th pixel of the attenuation image to the $i$-th projection datum. The inverse problem is to recover $\bx$ from the measurement data $\by^{\delta}$ perturbed by unknown noises, viz. $\by^{\delta}=\by+\boldsymbol{\eta}$. In our CT application, we consider the attenuation image with $256\times256$ pixels, and the sinogram of projection data is of size $256$ pixels by $60$ views. As a result, we have $N=256^2,\, M=256\times60$.

§.§ Training convex U-net

The convex U-net for <ref> is trained according to the approach described in Subsection <ref>. The attenuation images are taken as synthetic phantoms consisting of randomly generated piecewise constant ellipses. We prepared 2600 synthetic phantoms, denoted as $\{\bx_s^*\mid s=1,\cdots,2600\}$, and Figure <ref> shows six of them.

six of the 2600 synthetic phantoms.

The pairs of input images and output labels $\{(\boldsymbol{z}_s,\boldsymbol{r}_s)\mid s=1,\cdots,2600\}$ are constructed according to equation (<ref>) and equation (<ref>), where we take $N_1=N_2=1300$, so that the first 1300 samples simulate the images with artifacts and the second 1300 samples simulate the artifact-free images. The input images of the first 1300 samples are constructed according to

\[
\boldsymbol{z}_s=F^{\dagger}(F\bx_s^*+\boldsymbol{\eta}),
\]

where we add $0-10\%$ Gaussian noise to simulate the real situation of practical measurements. To be specific,

\boldsymbol{\eta}=\eta\cdot(F\bx_s^*)\cdot \mathcal{N}(0,1),

where $\eta$ denotes a random number between 0 and 0.1, and $\mathcal{N}(0,1)$ denotes the Gaussian noise with zero mean and unit variance. Figure <ref> (a) shows the simulated measurement data with Gaussian noises for the six synthetic phantoms displayed in Figure <ref>. The pseudo inverse $F^\dagger$ for the radon transform is achieved by the algebraic reconstruction technique (ART) <cit.>, which is a Kaczmarz iterative scheme for the linear system (<ref>). The ART algorithm has the following formula:

\begin{equation*}\label{eqn3.2}
\bx^{(n)}=\bx^{(n-1)}+\frac{y^\delta_i-\mathbf{a}_i^T\bx^{(n-1)}}{\|\mathbf{a}_i\|_2^2}\,\mathbf{a}_i,\qquad n\in\mathbb{N}^+,\ i=1,2,\cdots,M\,,
\end{equation*}

where $\mathbf{a}_i$ denotes the $i$-th row of the projection matrix $F$, $\mathbf{a}_i=(a_{i1},\cdots,a_{iN})^T$, and $y^\delta_i$ denotes the $i$-th component of the measurement data. One round of ART iterations implies one complete sweep of the measurement data $y^\delta_i$, with $i$ going from $1$ to $M$. In our application, the pseudo inverse $F^\dagger$ is constructed by 5 rounds of ART iterations, and the initial guess is taken as $\bx^{(0)}=\left(\frac{1}{N}\right)_{N\times1}$. Figure <ref> (b) shows the input images $\boldsymbol{z}_s=F^{\dagger}(F\bx_s^*+\boldsymbol{\eta})$ corresponding to the six synthetic phantoms displayed in Figure <ref>. Since ART and <ref> are both iterative approaches, the artifacts in the images generated by ART is consistent with the artifacts to be regularized in the <ref> algorithm.

Measurement data and input images with artifacts corresponding to the six synthetic phantoms displayed in Figure <ref>. (a) Measurement data with Gaussian noises: $F\bx_s^*+\boldsymbol{\eta}$; (b) input images with artifacts constructed by 5 round of ART iterations on the measurement data.

We separate the 2600 prepared samples into three sets: the training set of 2000 samples, $\{(\boldsymbol{z}_s,\boldsymbol{r}_s)\mid s\in \{1,\cdots,1000\}\cup \{1301,\cdots2300 \} \}$, the validation set of 400 samples, $\{(\boldsymbol{z}_s,\boldsymbol{r}_s)\mid s\in \{1001,\cdots,1200\}\cup \{2301,\cdots2500 \} \}$, and the test set of 200 samples, $\{(\boldsymbol{z}_s,\boldsymbol{r}_s)\mid s\in \{1201,\cdots,1300\}\cup \{2501,\cdots2600 \} \}$. Each set consists of two types of input images: half of images with artifacts, and half of artifact-free images. The convex U-net $\Phi^c_{\Theta}$ is trained according to Algorithm <ref>, where the regularization parameter $\lambda$ is taken as $5\times10^{-4}$, the batch size is set as $n_0=|\mathcal{I}_t|=10$, and the learning rate is $5\times10^{-4}$. Figure <ref> shows the convergence plot in the training process, where the blue curve illustrates the convergence of $\mathcal{L}_t(\Theta)$ on training set, and the red curve shows its performance on validation set. Since we employ a mini-batch optimization approach with random shuffling, and the non-negativity constraint is imposed after every iteration, the cost function $\mathcal{L}_t(\Theta)$ has an oscillating behavior. Finally, we evaluate the trained convex U-net on the test set of 200 samples. The mean squared error (MSE) is used to measure the performance of reconstructions, and it is defined in the following way,

\[
\mathrm{MSE}=\frac{1}{N}\|\Phi^c_{\Theta}(\boldsymbol{z}_s)- \boldsymbol{r}_s \|_2^2\,,\qquad \boldsymbol{z}_s,\, \boldsymbol{r}_s \in\mathbb{R}^N\,.
\]

Convergence plot of $\mathcal{L}_t(\Theta)$ in the training process. The blue curve illustrates the convergence on training set, and the red curve shows the performance on validation set.

Mean squared error (MSE) on the test set of 200 samples. The first 100 samples correspond to input images with artifacts, and the second 100 samples correspond to artifact-free images with $\boldsymbol{r}_s=\boldsymbol{0}$. (a) MSE for the convex U-net $\Phi_{\Theta}^c$; (b) MSE for the general U-net $\Phi_{\Theta}$.

Figure <ref> (a) shows the values of MSE for the trained convex U-net on the 200 test samples. In next subsection, we will illustrate the performance of <ref> and compare it with the standard NETT algorithm <cit.>. Since the NETT  algorithm does not require convexity of the neural network, we train a general U-net without the constraint of convexity, denoted as $\Phi_{\Theta}$. The training set, validation set, and test set are taken the same as those of convex U-net $\Phi_{\Theta}^c$. Figure <ref> (b) shows the values of MSE for the general U-net on the 200 test samples. It concludes that the performance of the convex U-net $\Phi_{\Theta}^c$ is almost comparable with that of the general U-net $\Phi_{\Theta}$ on the test set.

§.§ Reconstruction results

The <ref>  algorithm is implemented according to Subsection <ref>.  With the convex U-net $\Phi_{\Theta}^c$, we build the uniformly convex regularizer $\Reg(\bx)$ for <ref> as shown in Subsection <ref>, where we set $a=10^{-3}$. The sequence of regularization parameters $\alpha_n$ is taken as $2^{-n}$, $n\in\mathbb{N}^+$; the initial guess is $\bx_0=\frac{1}{N}$ with $N=256^2$, and $\bxi_0 \in\partial \Reg(\bx_0)$. The stopping rule is the discrepancy principle as shown in equation (<ref>), where we take the constant $\tau=1.01$.

§.§.§ Example 1: synthetic phantom with non-overlapping ellipses

Figure <ref> (a) shows the true model of attenuation image. It is a synthetic phantom of the same type as the phantom images in the training set and validation set, but it is not contained in those phantom images.

Example 1: synthetic phantom. (a) True model of attenuation image; (b) measurement data with $5\%$ Gaussian noises; (c) reconstruction result by <ref>; (d) reconstruction result by ART.

Figure <ref> (b) shows the simulated measurement data, where we add $5\%$ Gaussian noises with zero mean and unit variance. The <ref> algorithm is performed to reconstruct the attenuation image from the noisy measurement data. Figure <ref> (c) shows the reconstruction result. Figure <ref> (d) provides the reconstruction result by ART, where we have performed 5 rounds of ART iterations to get the solution. To quantitatively evaluate the performance of reconstruction, we compute the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM), where for both a higher value means a better reconstruction (note that the range of SSIM is between $0$ and $1$); the results are shown in the second and third columns of Table <ref>. The <ref> algorithm provides much better reconstruction for the synthetic phantom. It is capable of removing artifacts due to under-sampling and data noises while preserving the resolution of ellipses in the reconstructed image.

As a comparison, we provide reconstruction results by a standard iterated Tikhonov (SIT) method and by the NETT algorithm <cit.>, respectively. In the standard iterated Tikhonov (SIT) method, the regularizer $\Reg$ is taken as $\Reg(\bx)=\|\bx\|_2^2$ without the convex neural network in <ref>. The iteration is the same as that of <ref>, while it has a closed form in this standard situation:

\begin{equation}\tag*{SIT}\label{SIT}
\bx_n^\delta=\bx_{n-1}^\delta-\left(F^TF+\hat{\alpha}_n I\right)^{-1}F^T\left(F\bx_{n-1}^\delta-\by^\delta\right),
\end{equation}

where $F^T$ denotes the transpose of $F$, $I$ is the identity matrix, and $\hat{\alpha}_n=2M \alpha_n$; recall that $\alpha_n=2^{-n}$, and $M=256\times60$. The initial guess and the stopping rule are the same as those of <ref>. Figure <ref> shows the reconstruction result by <ref>, and Table <ref> lists the values of PSNR and SSIM for this reconstruction. It shows that the <ref> algorithm has a better performance than ART, but the improvement is inadequate comparing to the performance of <ref>.

In the NETT algorithm, the attenuation image is recovered by solving the following minimization problem,

\begin{equation}\tag*{NETT}\label{NETT2}
\bx_\alpha^\delta=\underset{\bx\in X}{\argmin}\frac{1}{2M}\|F \bx - \by^\delta\|_2^2+\alpha\|\Phi_{\Theta}(\bx) \|_2^2,
\end{equation}

where $\Phi_\Theta$ denotes the general U-net trained in section <ref>, and $\alpha$ is a fixed parameter controlling the amount of regularization. As stated above, $\Phi_\Theta$ is trained without the constraint of convexity, since the <ref> algorithm does not require convexity of neural network in general. We perform the <ref> algorithm with a series of different values of $\alpha$, although we are not trying to tune the parameter $\alpha$ exhaustively. Figure <ref> shows the reconstruction results with $\alpha=0.001, 0.01, 0.05$ and $0.1$, respectively. Table <ref> lists the values of PSNR and SSIM for those reconstruction results. It shows that the <ref> algorithm achieves the best reconstruction among the four solutions as $\alpha=0.05$. Comparing the results in Figure <ref> (c) and Figure <ref>, and considering the values of PSNR and SSIM listed in Table <ref>, we conclude that the <ref> algorithm can achieve a comparable reconstruction as <ref> without the procedure of tuning the regularization parameter manually.

Example 1: synthetic phantom. Reconstruction result by the standard iterated Tikhonov (SIT) method.

Example 1: synthetic phantom. Reconstruction results by <ref> with $\alpha=0.001, 0.01, 0.05$ and $0.1$, respectively. (a) $\alpha=0.001$; (b) $\alpha=0.01$; (c) $\alpha=0.05$; (d) $\alpha=0.1$.

2*       2*iNETT       2*ART      2*SIT      4cc]@c@NETT

$\alpha=0.001$     $\alpha=0.01$      $\alpha=0.05$      $\alpha=0.1$

PSNR                  21.30       14.48      17.97        16.85     18.23     19.34     18.67

SSIM                  0.84         0.54        0.66         0.62      0.73      0.79      0.79

(Example 1) PSNR and SSIM for the reconstruction results by <ref>, ART, <ref> and <ref>, respectively.

§.§.§ Example 2: synthetic phantom with overlapping ellipses

Figure <ref> (a) shows the true model of attenuation image. Again, it is a synthetic phantom of the same type as the phantom images in the training and validation sets, but it is not contained in them.

Example 2: synthetic phantom. (a) True model of attenuation image; (b) measurement data with $5\%$ Gaussian noises; (c) reconstruction result by <ref>; (d) reconstruction result by ART.

Figure <ref> (b) shows the simulated measurement data, where we add $5\%$ Gaussian noises with zero mean and unit variance. Figure <ref> (c) shows the reconstruction result by the <ref> algorithm. Figure <ref> (d) provides the reconstruction result by ART, where we have performed 5 rounds of ART iterations to get the solution. To quantitatively evaluate the performance of reconstruction, we compute the values of PSNR and SSIM as well; the results are shown in the second and third columns of Table <ref>. It shows that the <ref> algorithm provides much better reconstruction for the synthetic phantom. It is capable of removing artifacts in the reconstructed image and achieving higher values of PSNR and SSIM.

As a comparison, we provide reconstruction results by <ref> and <ref>, respectively. The solution by <ref> is shown in Figure <ref>, and the values of PSNR and SSIM are listed in the fourth column of Table <ref>. The <ref> algorithm has a better performance than ART, but the improvement is inadequate comparing to the performance of <ref>. The reconstruction results by <ref> are shown in Figure <ref>, where we perform the <ref> algorithm with $\alpha=0.001, 0.01, 0.05$ and $0.1$, respectively. The PSNR and SSIM of these reconstruction results are listed in Table <ref>. Comparing the results in Figure <ref> (c) and Figure <ref>, and considering the values of PSNR and SSIM listed in Table <ref>, we conclude that the <ref> algorithm can achieve a comparable (if not better) reconstruction as <ref> without tuning the regularization parameter exhaustively.

Example 2: synthetic phantom. Reconstruction result by the standard iterated Tikhonov (SIT) method.

Example 2: synthetic phantom. Reconstruction results by <ref> with $\alpha=0.001, 0.01, 0.05$ and $0.1$, respectively. (a) $\alpha=0.001$; (b) $\alpha=0.01$; (c) $\alpha=0.05$; (d) $\alpha=0.1$.

2*       2*iNETT       2*ART      2*SIT      4cc]@c@NETT

$\alpha=0.001$     $\alpha=0.01$      $\alpha=0.05$      $\alpha=0.1$

PSNR                  25.35       17.26      21.85        19.51     22.46     22.36     21.66

SSIM                  0.91         0.63        0.77         0.70      0.81      0.84      0.83

(Example 2) PSNR and SSIM for the reconstruction results by <ref>, ART, <ref> and <ref>, respectively.

§.§.§ Example 3: attenuation images of lung and myocardium

In this example, we consider reconstructions of attenuation images which are quite different from the synthetic phantoms in the training process.

Example 3: lung imaging. (a) Attenuation image of lung; (b) measurement data with $5\%$ Gaussian noises; (c) reconstruction result by <ref>; (d) reconstruction result by ART.

The true models of attenuation images are shown in Figure <ref> (a) and Figure <ref> (a), respectively, where Figure <ref> (a) shows the attenuation image of lung, and Figure <ref> (a) shows the attenuation image of myocardium. In the <ref> algorithm, we employ the same convex U-net $\Phi^c_{\Theta}$ trained on synthetic phantoms as described in section <ref>. The reconstruction results of lung imaging are provided in Figure <ref>, where Figure <ref> (b) shows the measurement data with $5\%$ Gaussian noises, Figure <ref> (c) shows the recovered solution of <ref>, and Figure <ref> (d) shows the recovered solution of ART. Correspondingly, the reconstruction results of myocardium imaging are provided in Figure <ref>.

Example 3: myocardium imaging. (a) Attenuation image of myocardium; (b) measurement data with $5\%$ Gaussian noises; (c) reconstruction result by <ref>; (d) reconstruction result by ART.

In addition, the values of PSNR and SSIM for these reconstruction results are listed in Table <ref>. It shows that the <ref> algorithm provides much better solutions in both lung imaging and myocardium imaging, removing artifacts in the reconstructed images and achieving higher values of PSNR and SSIM. In practical applications, one can further improve the performance of <ref> by re-training the convex neural network $\Phi^c_{\Theta}$ according to the types of attenuation images to be reconstructed, e.g. using images of lungs to train $\Phi^c_{\Theta}$ in lung imaging and using images of myocardia to train $\Phi^c_{\Theta}$ in myocardium imaging. In this example, we use the convex U-net trained on synthetic phantoms for the reconstruction of lung and myocardium images, which is partly due to the lack of appropriate training data, and partly to illustrate the generalization ability of the <ref> algorithm.

2*     2cc]@c@ lung imaging        2cc]@c@ myocardium imaging

iNETT       ART                         iNETT       ART

PSNR                             20.02        16.46                           19.95     16.08

SSIM                              0.70         0.51                               0.69      0.46

(Example 3) List of PSNR and SSIM for the reconstruction results.

§ CONCLUSION

We have proposed a non-stationary iterated network Tikhonov (iNETT) method for the solution of ill-posed inverse problems. Motivated by the network Tikhonov (NETT) method <cit.>, we employ a data-driven regularization term including a neural network architecture, where the neural network is trained to penalize artifacts in the recovered solution. The iNETT method is then developed by combining the strategy of neural-network based regularization with the non-stationary iterated Tikhonov method. The main advantage of iNETT is that it avoids the difficult task of estimating the optimal regularization parameter, while keeps the virtue of the data-driven regularizer in NETT. We provide numerical examples to illustrate applications of iNETT in 2D computerized tomography. It shows that the iNETT method can achieve comparable (if not better) results as NETT without tuning the regularization parameter exhaustively.

To achieve the theoretical convergence of iNETT, we require that the neural-network based regularizer in iNETT is uniformly convex. As a result, we develop theories and algorithms for the construction of convex and uniformly convex neural networks. Given a general neural network, we prescribe sufficient conditions to achieve a related neural network which is component-wise convex or uniformly convex. Our formulation can embrace state-of-the-art convolutional neural networks. As a concrete example, we provide rigorous formulas for the U-net architecture, and explain the approaches to obtain convex and uniformly convex U-net architectures, which are successfully used in iNETT for the computations of 2D computerized tomography. The tool of convex and uniformly convex neural networks shall have many interesting applications in the future study.

For the iNETT method itself, there are several directions to further improve the algorithm. (1) Study more efficient neural network architectures with different data-fidelity norms, which can be more stable under random noises with non-Gaussian distributions, e.g. the $\ell^p$-$\ell^q$ norm regularization <cit.>. (2) Test alternative training strategies for the uniformly convex neural network regularizer, such as the adversarial strategy in <cit.>. (3) Learn the sequence $\{\alpha_n\}_n$ in iNETT, mimicking  the procedure developed in <cit.>. (4) Devise nonlocal differential-like operators  in the regularization term to enforce some a-priori knowledge on the recovered solution <cit.>. (5) Generalize iNETT to a fully learned network architecture, in the same spirit of <cit.>, e.g. learn the neural network regularizer $\Reg(\bx)$ in the optimization process.

§ ACKNOWLEDGMENTS

Davide Bianchi is supported by NSFC (grant no. 12250410253). Wenbin Li is supported by Natural Science Foundation of Shenzhen (grant no. JCYJ20190806144005645) and NSFC (grant no. 41804096).

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

titleMoDL: Model-based deep learning architecture for inverse
journaltitleIEEE Trans. Med. Imag.

booktitleProceedings of the 34th International Conference on
Machine Learning
seriesProceedings of Machine Learning Research
titleInput Convex Neural Networks

titleDigital image restoration

Multidisciplinary Digital Publishing Institute

titleDiscretization of learned NETT regularization for solving
inverse problems
journaltitleJ. Imaging

Taylor & Francis

titleDeep learning for photoacoustic tomography from sparse data
journaltitleInverse Probl. Sci. Eng.

booktitlePhotons Plus Ultrasound: Imaging and Sensing 2019
titleNETT regularization for compressed sensing photoacoustic

De Gruyter

titleA data and knowledge driven approach for SPECT using
convolutional neural networks and iterative algorithms
journaltitleJ. Inverse Ill-Posed Probl.

Cambridge University Press

titleSolving inverse problems using data-driven models
journaltitleActa Numer.

IOP Publishing

titleIterative total variation schemes for nonlinear inverse
journaltitleInverse Problems

IOP Publishing

titleDeep unfolding of a proximal interior point method for image
journaltitleInverse Problems

titleGraph Laplacian for image deblurring

IOP Publishing

titleIterated fractional Tikhonov regularization
journaltitleInverse Problems

booktitle2021 21st International Conference on Computational
Science and Its Applications (ICCSA)
titleGraph approximation and generalized Tikhonov regularization
for signal deblurring

titleOn generalized iterated Tikhonov regularization with
operator-dependent seminorms

titleUniformly convex functions on Banach spaces
journaltitleProc. Amer. Math. Soc.

Wiley Online Library

titleIterated Tikhonov regularization with a general penalty term
journaltitleNumer. Linear Algebra Appl.

IOP Publishing

titleModulus-based iterative methods for constrained
$\ell^p$–$\ell^q$ minimization
journaltitleInverse Problems

titleRegularization preconditioners for frame-based image
deblurring with reduced boundary artifacts
journaltitleSIAM J. Sci. Comput.

Optical Society of America

titleLow-dose CT via convolutional neural network
journaltitleBiomed. Opt. Express

titleOn nondecreasing sequences of regularization parameters for
nonstationary iterated Tikhonov
journaltitleNumer. Algor.

titleA guide to convolution arithmetic for deep learning

titleOn the choice of the regularization parameter for iterated
Tikhonov regularization of III-posed problems
journaltitleJ. Approx. Theory

Springer Science & Business Media

titleRegularization of inverse problems

titleNonlocal operators with applications to image processing
journaltitleMultiscale Model. Simul.

Academic Press

titleAlgebraic reconstruction techniques (ART) for
three-dimensional electron microscopy and X-ray photography
journaltitleJ. Theoret. Biol.

titleLimitations of the L-curve method in ill-posed problems

titleNonstationary iterated Tikhonov regularization
journaltitleJ. Optim. Theory Appl.

titleAnalysis of discrete ill-posed problems by means of the
journaltitleSIAM Rev.

booktitleProc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern
Recognit. (CVPR)
titleDeep residual learning for image recognition

De Gruyter

titleResidual U-net convolutional neural network architecture for
low-dose CT denoising
journaltitleCurr. Dir. Biomed. Eng.

titleDeep Learning: An Introduction for Applied Mathematicians
journaltitleSIAM Rev.

family=Van DerMaaten,

booktitleProc. IEEE Comput. Soc. Conf. Comput. Vis. Pattern
Recognit. (CVPR)
titleDensely connected convolutional networks

booktitleInternational conference on machine learning
titleBatch normalization: Accelerating deep network training by
reducing internal covariate shift

IOP Publishing

titleNonstationary iterated Tikhonov regularization for ill-posed
problems in Banach spaces
journaltitleInverse Problems

titleNonstationary iterated Tikhonov regularization in Banach
spaces with uniformly convex penalty terms
journaltitleNumer. Math.

titleAdam: A method for stochastic optimization

IOP Publishing

titleNETT: Solving inverse problems with deep neural networks
journaltitleInverse Problems

IOP Publishing

titleA stochastic gradient descent approach with
partitioned-truncated singular value decomposition for large-scale inverse
problems of magnetic modulus data
journaltitleInverse Problems

Springer Science & Business Media

titleClassical Banach spaces II: function spaces

titleAdversarial regularizers in inverse problems
journaltitleAdvances in neural information processing systems

Multidisciplinary Digital Publishing Institute

titleA Green Prospective for Learned Post-Processing in
Sparse-View Tomographic Reconstruction
journaltitleJ. Imaging

titleLearned convex regularizers for inverse problems

IOP Publishing

titleAugmented NETT regularization of inverse problems
journaltitleJ. Phys. Commun.

IOP Publishing

titleDeep synthesis network for regularizing inverse problems
journaltitleInverse Problems

titleAn iterative regularization method for total variation-based
image restoration
journaltitleMultiscale Model. Simul.

titleOld and new parameter choice rules for discrete ill-posed
journaltitleNumer. Algorithms

given=Roger J-B,

Springer, Berlin, Heidelberg

editionCorrected 3rd printing 2009
seriesGrundlehren der mathematischen Wissenschaften
titleVariational analysis

booktitleProc. Med. Image Comput. Comput. Assist. Interv.
titleU-net: Convolutional networks for biomedical image

booktitleProceedings of Thirty Third Conference on Learning
seriesProceedings of Machine Learning Research
titleHow Good is SGD with Random Shuffling?

titleVariational methods in imaging

booktitleJoint European Conference on Machine Learning and
Knowledge Discovery in Databases
titleThe curious case of convex neural networks

titleDeep ADMM-Net for compressive sensing MRI
journaltitleAdvances in neural information processing systems

titleData-Driven Mirror Descent with Input-Convex Neural Networks

titleDeep-learning schemes for full-wave nonlinear inverse
scattering problems
journaltitleIEEE Trans. Geosci. Remote Sensing

booktitleProceedings of the IEEE conference on computer vision and
pattern recognition
titleISTA-Net: Interpretable optimization-inspired deep network
for image compressive sensing

World scientific

titleConvex analysis in general vector spaces

titleOn uniformly convex functions
journaltitleJ. Math. Anal. Appl.