# Stability Bounds for the Unfolded Forward-Backward Algorithm

Emilie Chouzenoux, Cecile Della Valle, Jean-Christophe Pesquet

► To cite this version:

Emilie Chouzenoux, Cecile Della Valle, Jean-Christophe Pesquet. Stability Bounds for the Unfolded Forward-Backward Algorithm. Journal of Mathematical Imaging and Vision, 2025, 67 (48), 10.1007/s10851-025-01258-6. hal-04854255v2

HAL Id: hal-04854255

https://inria.hal.science/hal-04854255v2

Submitted on 4 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stability Bounds for the Unfolded Forward-Backward Algorithm

Emilie Chouzenoux¹, Cécile Della Valle²,
Jean-Christophe Pesquet¹*

¹Center for Visual Computing, Inria, CentraleSupélec, University Paris Saclay, Gif-sur-Yvette, 91190, France.

²Université de Paris, Paris Sorbonne Université, Paris, 75005, France.

*Corresponding author(s). E-mail(s):

jean-christophe.pesquet@centralesupelec.fr;

Contributing authors: emilie.chouzenoux@inria.fr;

cecile.della-valle@parisdescartes.fr;

## Abstract

We consider a neural network architecture designed to solve inverse problems where the degradation operator is linear and known. This architecture is constructed by unrolling a forward-backward algorithm derived from the minimization of an objective function that combines a data-fidelity term, a Tikhonov-type regularization term, and a potentially nonsmooth convex penalty. The robustness of this inversion method to input perturbations is analyzed theoretically. Ensuring robustness complies with the principles of inverse problem theory, as it ensures both the continuity of the inversion method and the resilience to small noise - a critical property given the known vulnerability of deep neural networks to adversarial perturbations. A key novelty of our work lies in examining the robustness of the proposed network to perturbations in its bias, which represents the observed data in the inverse problem. Additionally, we provide numerical illustrations of the analytical Lipschitz bounds derived in our analysis.

**Keywords:** Stability analysis, neural networks, deep unrolling, forward-backward, proximal-gradient, inverse problems
# 1 Introduction

## Inverse problems.

A large variety of inverse problems consists of inverting convolution operators such as signal/image restoration [1, 2], tomography [3], Fredholm equation of the first kind [4], or inverse Laplace transform [5].

Let $T$ be a bounded linear operator from a separable Hilbert space $\mathcal{X}$ to a Hilbert space $\mathcal{Y}$. The problem consists in estimating $\bar{x} \in \mathcal{X}$ from observed data

$$y = T\bar{x} + w \quad (1)$$

where $w$ corresponds to an additive measurement noise. The above problem is often ill-posed i.e., a solution might not exist, might not be unique, or might not depend continuously on the data.

## Variational problem.

The well-posedness of the inverse problem defined by (1) is retrieved by regularization. Here we consider Tikhonov-type regularization. Let $\tau \in ]0, +\infty[$ be the regularization parameter. Solving the inverse problem (1) with such regularization, leads to the resolution of the following optimization problem:

$$\min_{x \in C} J_{\tau}(x), \quad (2)$$

where

$$(\forall x \in \mathcal{X}) \quad J_{\tau}(x) = \frac{1}{2} \|Tx - y\|^2 + \frac{\tau}{2} \|Dx\|^2, \quad (3)$$

$C$ is a nonempty closed convex subset of $\mathcal{X}$ encoding some prior knowledge, e.g. some range constraint or sparsity pattern, while $D$ acts as a derivative operator of some order $r \ge 0$. Often, we have an a priori of smoothness on the solution which justifies the use of such a derivative-based regularization.<fnref n="1" /> Problem (2) is actually an instance of the more general class of problems stated below, encountered in many signal/image processing tasks:

$$\min_{x \in \mathcal{X}} J_{\tau}(x) + \mu g(x), \quad (4)$$

where $\mu \in [0, +\infty[$ is an additional regularization constant and $g$ is a proper lower-continuous convex function from Hilbert space $\mathcal{X}$ to $]-\infty, +\infty[$]. Indeed, Problem (2) corresponds to the case when $g$ is the indicator function $\iota_C$ of set $C$.

## Neural network.

We focus our attention on seeking for a solution to the addressed inverse problem through nonlinear approximation techniques making use of neural networks. Thus, instead of considering the solution to the regularized problem (4), we define the solution to the inverse problem (1) as the output of a neural network, whose structure is similar to a recurrent network [6, 7].

<footnote n="1">¹ $\| \cdot \|$ denotes the norm of the considered Hilbert space.</footnote>
Namely, by setting an initial value $x_0$, we are interested in the following m-layers neural network where $m \in \mathbb{N} \setminus \{0\}$:

$$ 
\left\{
\begin{array}{ll}
\textbf{Initialization:} & \\
\qquad b_0 = T^*y, & (5) \\
\textbf{Layer } n \in \{1, \dots, m\}: & \\
\qquad x_n = \text{prox}_{\lambda_n \mu_n g}(x_{n-1} - \lambda_n (T^*T + D_n^* D_n)x_{n-1} + \lambda_n b_0). &
\end{array}
\right.
$$

Hereabove, $\text{prox}_\varphi$ stands for the proximity operator of a lower-semicontinuous proper convex function $\varphi$ (see Section 2) and, for every $n \in \{1, \dots, m\}$, $D_n$ is a bounded linear operators from $\mathcal{X}$ to some Hilbert space $\mathcal{Z}$, and $\lambda_n$ is positive real. Throughout this paper, $L^*$ denotes the adjoint of a bounded linear operator $L$ defined on suitable Hilbert spaces. The overall structure of the network is shown in Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Global architecture of the m-layers neural network (5).</figcaption></figure>

A simple choice consists in setting, for every $n \in \{1, \dots, m\}$, $D_n = \sqrt{\tau_n} D$ where $\tau_n$ is a positive constant. The constants $(\mu_n)_{1 \le n \le m}$, $(\tau_n)_{1 \le n \le m}$, and possibly $(\lambda_n)_{1 \le n \le m}$ can be learned during training. Then, Model (5) can be viewed as unrolling m iterations of an optimization algorithm, so leading to Algorithm 1. If additionally, we set $\mu_n \equiv \mu$ and $\tau_n \equiv \tau$, Algorithm 1 identifies with the forward-backward algorithm [8, 9] (also called proximal gradient algorithm) applied to the variational problem (4).

<caption>**Algorithm 1** Proximal forward-backward splitting method</caption>

1: Set $x_0$,
2: **for** $n = 1, 2, \dots, m$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;Set $\lambda_n, \tau_n, \mu_n$,
4: &nbsp;&nbsp;&nbsp;&nbsp;$x_n = \text{prox}_{\lambda_n \mu_n g}(x_{n-1} - \lambda_n \nabla J_{\tau_n}(x_{n-1}))$,
5: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
6: **return** $x_m$

Let $\bar{\chi} \ge 0$ denote the modulus of strong convexity of the regularization term $g$ (with $\bar{\chi} > 0$ if $g$ is strongly convex). We have thus $g = g_0 + \bar{\chi} \|\cdot\|^2/2$ where $g_0$ is a proper lower-semicontinuous convex function from $\mathcal{X}$ to $]-\infty, +\infty]$, and Model (5)
can be equivalently expressed as

$$ 
\left\{
\begin{array}{l}
\text{Initialization:} \\
b_0 = T^*y, \\
\text{Layer } n \in \{1, \dots, m\}: \\
x_n = R_n(W_n x_{n-1} + V_n b_0),
\end{array}
\right.
\qquad (6)
$$

where, for every $n \in \{1, \dots, m\}$,

$$ R_n = \underset{\lambda_n \in \{\lambda_n\}}{\operatorname{prox}} \frac{\lambda_n \mu_n}{1 + \lambda_n \chi_n} g_0 \qquad (7) $$

$$ W_n = \frac{1}{1 + \lambda_n \chi_n} (\mathbf{1} - \lambda_n (T^*T + D_n^*D_n)) \qquad (8) $$

$$ V_n = \frac{\lambda_n}{1 + \lambda_n \chi_n} \mathbf{1} \qquad (9) $$

$$ \chi_n = \mu_n \bar{\chi}, \qquad (10) $$

and $\mathbf{1}$ denotes the identity operator.

From a theoretical standpoint, little is known about the theoretical properties of Model (6) in relation with the minimization of the original regularized objective function (4). The main challenges are that (i) the number of iterations $m$ is predefined, and (ii) the parameters can vary along the iterations. However, insight can be gained by quantifying the robustness of Model (6) to perturbations on its initialization $x_0$ and on its bias $b_0$, through an accurate estimation of its Lipschitz properties.

## Related works and contributions.

There has been a plethora of techniques developed to invert models of the form (1). Among these methods, Tikhonov-type methods are attractive from a theoretical viewpoint, especially because they provide good convergence rate as the noise level decreases, as shown in [10] or [11]. Optimization techniques [9] are then classically used to solve Problem (4). However, limitations of such methods may be encountered in their implementation. Indeed, certain parameters such as gradient descent stepsizes or the regularization coefficient need to be set, as discussed in [12] or [13]. The latter parameter depends on the noise level, as shown in [14], which is not always easy to estimate. In practical scenarios, methods like the L-curve method (see [15]) can be used to set the regularization parameter. However, when dealing with more than one parameter, their application becomes computationally expensive. Moreover, incorporating constraints on the solution may be difficult in such approaches, and often reduces to projecting the resulting solution onto the desired set. These reasons justify the use of a neural network structure to avoid laborious calibration of the parameters and to easily incorporate constraints on the solution.

The use of neural networks for solving inverse problems has become increasingly popular, especially in the image processing community. A rich panel of approaches have been proposed, either adapted to the sparsity of the data [16, 17], or mimicking variational models [18, 19], or iterating learned operators [20–24], or adapting
Tikhonov method [25]. The successful numerical results of the aforementioned works raise two theoretical questions: in cases when these methods are based on the iterations of an optimization algorithm, do they converge (in the algorithmic sense)? Are these inversion methods stable or robust?

In iterative approaches, a regularization operator is learned, either in the form of a proximity (or denoiser) operator as [20, 21, 24], of a regularization term [25], of a pseudodifferential operator [26], or of its gradient [3, 27]. Strong connections also exist with Plug-and-Play methods [22, 28, 29], where the regularization operator is a pretrained neural network. Such strategies have in particular enabled high-quality imaging restoration or tomography inversion [3]. Here, the nonexpansiveness of the neural network is a core property to establish convergence of the algorithm [3, 22], although recent works have provided convergence guarantees using alternative tools [30].

Other recent works solve linear inverse problems by unrolling the optimization iterative process in the form of a network architecture as in [31–33]. Here the number of iterations is fixed, instead of iterating until convergence, and the network is often trained in an end-to-end fashion. Neural network frameworks offer powerful differential programming capabilities for learning hyper-parameters in an unrolled optimization algorithm as in [34–38]. With respect to a traditional optimization method applied to (4)-(3), the use of these data-driven techniques provides an increased flexibility leading to a better solution of the inverse problem.

All of the above strategies have shown very good numerical results. However, few studies have been conducted on their theoretical properties, especially on their stability. The study of the robustness of such a structure is often empirical, based on a series of numerical tests, as performed in [38, 39]. In [25], stability properties are established at low noise regime when using neural networks regularization. A fine characterization of the convergence conditions of recurrent neural network and of their stability via the estimation of a Lipschitz constant is done in [40, 41]. In particular, the Lipschitz constant estimated in [41] is more accurate than in basic approaches which often rely in computing the product of the norms of the linear weight operators of each layer as in [42, 43]. Thanks to the aforementioned works, proofs of stability/convergence have been demonstrated on specific neural networks applied to classification (e.g., [44, 45]), or inverse problems (e.g., [22, 23, 34]). The analysis carried out in this article is in the line of these references.

Our contributions in this paper are the following.

1. We propose a neural network architecture to solve the inverse problem (1), incorporating the proximity operator of a potentially nonsmooth convex function as its activation function. This design enables the imposition of constraints on the sought solution. A key advantage of this architecture lies in its foundation on the unrolling of a forward-backward algorithm, which ensures the network structure is both interpretable and frugal in terms of the number of parameters to be learned.

2. We study theoretically and numerically the stability of the so-built neural network. The sensitivity analysis is conducted with respect to the observed data $y$, which corresponds to a bias term in each layer of Model (6). This analysis is more general than the one performed in [34], in which only the impact of the initialization was considered.
## Outline.

The outline of the paper is as follows. In Section 2, we recall the theoretical background of our work and specify the notation, which follows the framework in [46]. Section 3 introduces a class of dynamical systems with a leakage factor, providing a more general framework than the neural network defined by (6). In Section 4, we establish the stability of the corresponding neural network defined on the product space $\mathcal{X} \times \mathcal{X}$, building upon the results from [40] and [41]. By stability, we mean that the network output remains controlled with respect to both its initial input $x_0$ and the bias term $T^*y$, ensuring that small differences or errors in these vectors are not amplified through the network. We also provide sufficient conditions for the averagedness of the network. Finally, we show how stability properties for network (6) can be deduced from these results. Section 5 presents numerical illustrations of the derived Lipschitz bounds, followed by concluding remarks in Section 6.

## 2 Notation

We introduce the elements of convex analysis we will be dealing with, in Hilbert spaces. We also cover the bits of monotone operator theory that will be needed throughout.

Let us consider the Hilbert space $\mathcal{X}$ endowed with the norm $\| \cdot \|$ and the scalar product $\langle \cdot, \cdot \rangle$. In the following, $\mathcal{X}$ shall refer to the underlying signal space. The notation $\| \cdot \|$ will also refer to the operator norm of bounded operators from $\mathcal{X}$ onto $\mathcal{X}$.

An operator $S: \mathcal{X} \to \mathcal{X}$ is nonexpansive if it is 1-Lipschitz, that is

$$ (\forall (x, y) \in \mathcal{X} \times \mathcal{X}) \quad \|Sx - Sy\| \le \|x - y\|. $$

Moreover, $S$ is said to be

1. firmly nonexpansive if

$$ (\forall (x, y) \in \mathcal{X} \times \mathcal{X}) \quad \|Sx - Sy\|^2 + \|(1-S)x - (1-S)y\|^2 \le \|x - y\|^2; $$

2. a Banach contraction if there exists $\kappa \in ]0, 1[$ such that

$$ (\forall (x, y) \in \mathcal{X} \times \mathcal{X}) \quad \|Sx - Sy\| \le \kappa \|x - y\|. \quad (11) $$

If $S$ is a Banach contraction, then the iterates $(S^n x)_{n \in \mathbb{N}}$ converge linearly to a fixed point of $S$ according to Picard's theorem. On the other hand, when $S$ is nonexpansive, the convergence is no longer guaranteed. A way of recovering the convergence of the iterates is to assume that $S$ is averaged, i.e., there exists $\alpha \in ]0, 1[$ and a nonexpansive operator $R: \mathcal{X} \to \mathcal{X}$ such that $S = (1 - \alpha)1 + \alpha R$. In particular, $S$ is $\alpha$-averaged if and only if

$$ (\forall (x, y) \in \mathcal{X} \times \mathcal{X}) \quad \|Sx - Sy\|^2 + \frac{1 - \alpha}{\alpha} \|(1-S)x - (1-S)y\|^2 \le \|x - y\|^2. $$
If $S$ has a fixed point and it is averaged, then the iterates $(S^n x)_{n \in \mathbb{N}}$ converge weakly to a fixed point. Note that $S$ is firmly nonexpansive if and only if it is $1/2$-averaged and that, if $S$ satisfies (11) with $\kappa \in ]0, 1[$, then it is $(\kappa + 1)/2$-averaged.

Let $\Gamma_0(\mathcal{X})$ be the set of proper lower semicontinuous convex functions from $\mathcal{X}$ to $]-\infty, +\infty]$. Then we define the proximity operator [46, Chapter 24] as follows.

**Definition 1.** Let $f \in \Gamma_0(\mathcal{X})$, $x \in \mathcal{X}$, and $\gamma > 0$. Then $\text{prox}_{\gamma f}(x)$ is the unique point that satisfies

$$ \text{prox}_{\gamma f}(x) = \underset{y \in \mathcal{X}}{\operatorname{argmin}} \left( f(y) + \frac{1}{2\gamma} \|x - y\|^2 \right). $$

The function $\text{prox}_{\gamma f}: \mathcal{X} \to \mathcal{X}$ is the proximity operator of $\gamma f$.

Finally, the proximity operator has the following property.

**Proposition 2.1.** [46, Proposition 12.28] The operators $\text{prox}_{\gamma f}$ and $1 - \text{prox}_{\gamma f}$ are firmly nonexpansive.

In the proposed neural network (6), the activation operator is a proximity operator. In practice, this is the case for most activation operators, as shown in [40]. The neural network (6) is thus a cascade of firmly nonexpansive operators and linear operators. If the linear part is also nonexpansive, bounds on the effect of a perturbation of the neural network or its iterates can be established.

# 3 Considered neural network architecture

In this section, Model (6) is reformulated as a virtual network, which takes as inputs the classical ones on top of a new one, which is a bias parameter.

## 3.1 Virtual neural network with leakage factor

To facilitate our theoretical analysis, we will introduce a virtual network making use of new variables $(z_n)_{n \in \mathbb{N}}$ in the product space $\mathcal{X} \times \mathcal{X}$. For every $n \in \mathbb{N} \setminus \{0\}$, we define the n-th layer of our virtual network as follows:

$$ z_n = \begin{pmatrix} x_n \\ b_n \end{pmatrix}, \quad z_n = Q_n(U_n z_{n-1}), \quad \text{with} \quad \left\{ \begin{array}{l} Q_n = \begin{pmatrix} R_n \\ 1 \end{pmatrix}, \\ U_n = \begin{pmatrix} W_n & V_n \\ 0 & \eta_n 1 \end{pmatrix}, \end{array} \right. \qquad (12) $$

and $W_n$ (resp. $V_n$) defined by (8) (resp. (9)). In our context, $R_n$ is defined by (7), but all the results in our stability analysis remain valid for any firmly nonexpansive operator such as the resolvent of a maximally monotone operator. Note that, in order to gain more flexibility, we have included positive multiplicative factors $(\eta_n)_{n \ge 1}$ on
the bias. Since, for every $n \ge 1$, $b_n = \eta_n \dots \eta_1 b_0$, cascading $m$ such layers yields

$$ 
\left\{
\begin{array}{l}
\text{Initialization:} \\
b_0 = T^*y, \\
\text{Layer } n \in \{1, \dots, m\}: \\
x_n = R_n(W_n x_{n-1} + \tilde{V}_n b_0) \\
\quad = \text{prox}_{\lambda_n \mu_n g}(x_{n-1} - \lambda_n (T^*T + D_n^* D_n)x_{n-1} + \lambda_n \eta_{n-1} \cdots \eta_0 b_0),
\end{array}
\right.
\qquad (13)
$$

where

$$ \tilde{V}_n = \frac{\lambda_n}{1 + \lambda_n \chi_n} \eta_{n-1} \cdots \eta_0 \mathbf{1} \quad \text{and} \quad \eta_0 = 1. \qquad (14) $$

We thus see that the network defined by Model (6) is equivalent to the virtual one when all the factors $\eta_n$ are equal to one. When $n \ge 1$ and $\eta_n < 1$, the parameters $\eta_n$ can be interpreted as a leakage factor.

**Remark 3.1.** In the original forward-backward algorithm, the introduction of $(\eta_n)_{n \ge 1}$ amounts to introducing an error $e_n$ in the gradient step, at iteration $n$, which is equal to

$$ e_n = \frac{\lambda_n}{1 + \lambda_n \chi_n} (\eta_{n-1} \cdots \eta_0 - 1)b_0. \qquad (15) $$

From known properties concerning the forward-backward algorithm [9], the standard convergence results for the algorithm are still valid provided that

$$ \sum_{n=2}^{+\infty} \frac{\lambda_n}{1 + \lambda_n \chi_n} |\eta_{n-1} \cdots \eta_0 - 1| < +\infty. \qquad (16) $$

## 3.2 Multilayer architecture

When we cascade the layers of the virtual neural network (12), the following triangular linear operator plays a prominent role:

$$ U = U_m \circ \dots \circ U_1 = \begin{pmatrix} W_{1,m} & \widetilde{W}_{1,m} \\ 0 & \eta_{1,m} \mathbf{1} \end{pmatrix}, \qquad (17) $$

where, for every $n \in \{1, \dots, m\}$ and $i \in \{1, \dots, n\}$

$$ \widetilde{W}_{i,n} = \sum_{j=i}^{n} \frac{\lambda_j}{1 + \lambda_j \chi_j} \eta_{i,j-1} W_{j+1,n} \qquad (18) $$

and, for every $i \in \{1, \dots, m+1\}$ and $j \in \{0, \dots, m\}$,

$$ W_{i,j} = \begin{cases} W_j \circ \dots \circ W_i & \text{if } j \ge i \\ \mathbf{1} & \text{otherwise,} \end{cases} \qquad (19) $$

$$ \eta_{i,j} = \begin{cases} \eta_j \cdots \eta_i & \text{if } j \ge i \\ 1 & \text{otherwise.} \end{cases} \qquad (20) $$
In particular, $W_{i,i} = W_i$ and $\eta_{i,i} = \eta_i$.

# 4 Stability and $\alpha$-averagedness

This section presents our stability results. Section 4.1 introduces the main assumptions and defines key quantities relevant to our analysis. Section 4.2 focuses on the stability analysis of the Virtual Neural Network defined in (12), considered as an operator $S$ from $\mathcal{X} \times \mathcal{X}$ to itself, with input $(x_0, b_0)$ and output $(x_m, b_m)$. The core result in Proposition 4.2 provides an estimate $\theta_m$ of its Lipschitz constant, which can be computed with quadratic complexity. Additionally, we derive simpler upper and lower bounds for the Lipschitz constant. Proposition 4.3 offers sufficient conditions for the VNN to be averaged. Notably, in the absence of a leakage factor, the VNN cannot be nonexpansive, as shown in Proposition 4.4. First estimates of the Lipschitz constant for network (6) are derived in Section 4.3. Section 4.4.1 aims to refine the previous Lipschitz bounds by focusing on the first output of the VNN, so coming back to the original network (6). This enables, in Section 4.4.2, an investigation of the stability of the network when its initialization $x_0$ is linked to the observation $y$ through a filtering relation—thus defining an operator from $\mathcal{X}$ to $\mathcal{X}$ (see Proposition 4.6).

## 4.1 Assumptions

We will make the following assumption on the degradation operator $T$ and the regularization operators $(D_n)_{1 \le n \le m}$:

Operators $T^*T$ and $(D_n^*D_n)_{1 \le n \le m}$ can be diagonalized in an orthonormal basis $(v_p)_p$ of $\mathcal{X}$.

The existence of such an eigenbasis is satisfied in the following cases:

* $T$ and $(D_n)_{1 \le n \le m}$ are compact operators (which guarantees the existence of their singular value decompositions), and $(D_n^*D_n)_{0 \le n \le m}$ commute pairwise if we set $D_0 = T$.

* $\mathcal{X} = \mathcal{Y}$ is the space $L^2([0, T])$ of square summable complex-valued functions defined on $[0, T]$, $T$ is a filter (i.e., circular convolutive operator) and, for every $n \in \{1, \dots, m\}$, $D_n$ is a single-input $c$-output filter that maps every signal $x \in \mathcal{X}$ to a vector $(D_{n,1}x, \dots, D_{n,c}x) \in \mathcal{X}^c$. In this case, the diagonalization is performed in the Fourier basis defined as

$$ (\forall p \in \mathbb{Z})(\forall t \in [0, T]) \quad v_p(t) = \frac{1}{\sqrt{T}} \exp \left( 2\pi i \frac{pt}{T} \right). $$

* $\mathcal{X} = \mathcal{Y}$ is the space $l^2(N)$ of complex-valued sequences $x = (x_n)_{0 \le n \le N-1}$ equipped with the standard Hermitian norm, $T$ is a discrete periodic filter and, for every $n \in \{1, \dots, m\}$, $D_n$ is a single-input $c$-output periodic filter that maps every signal $x \in \mathcal{X}$ to a vector $(D_{n,1}x, \dots, D_{n,c}x) \in \mathcal{X}^c$. The eigenbasis is associated with the discrete Fourier transform:

$$ (\forall (p, n) \in \{0, \dots, N-1\}^2) \quad v_p(n) = \frac{1}{\sqrt{N}} \exp \left( 2\pi i \frac{pn}{N} \right). $$
• The two previous examples extend to d-dimensional signals (d = 2 for images), defined on $[0, T_1] \times \dots \times [0, T_d]$ or $\{0, \dots, N_1 - 1\} \times \dots \times \{0, \dots, N_d - 1\}$.

Based on the above assumption, we define the respective nonnegative real-valued eigenvalues $(\beta_{T,p})_p$ and $(\beta_{D_n,p})_p$ of $T^*T$ and $D_n^*D_n$ with $n \in \{1, \dots, m\}$, as well as the following quantities, for every eigenspace index $p$ and $i \in \{1, \dots, n\}$,

$$ \beta_p^{(n)} = \frac{1}{1 + \lambda_n \chi_n} (1 - \lambda_n (\beta_{T,p} + \beta_{D_n,p})), \quad (21) $$

$$ \beta_{i,n,p} = \prod_{j=i}^{n} \beta_{p}^{(j)}, \quad (22) $$

$$ \tilde{\beta}_{i,n,p} = \sum_{j=i}^{n-1} \beta_{p}^{(n)} \dots \beta_{p}^{(j+1)} \frac{\lambda_j \eta_{i,j-1}}{1 + \lambda_j \chi_j} + \frac{\lambda_n \eta_{i,n-1}}{1 + \lambda_n \chi_n}, \quad (23) $$

with the convention $\sum_{i=n}^{n-1} = 0$. As limit cases, we have

$$ \beta_{n,n,p} = \beta_{p}^{(n)}, \quad \tilde{\beta}_{n,n,p} = \frac{\lambda_n}{1 + \lambda_n \chi_n}. \quad (24) $$

Note that $(\beta_p^{(n)}, v_p)_p$, $(\beta_{i,n,p}, v_p)_p$, and $(\tilde{\beta}_{i,n,p}, v_p)_p$ are the eigensystems of $W_n$, $W_{i,n}$ and $\tilde{W}_{i,n}$, defined by (8), (19), and (18), respectively.

## 4.2 Stability results for the virtual network

We first recall some existing results on the stability of neural networks [40, Proposition 3.6(iii)] [41, Theorem 4.2].

**Proposition 4.1.** Let $m > 1$ be an integer, let $(\mathcal{H}_i)_{0 \le i \le m}$ be non null Hilbert spaces. For every $n \in \{1, \dots, m\}$, let $U_n$ be a bounded linear operator from $\mathcal{H}_{n-1}$ to $\mathcal{H}_n$ and let $Q_n: \mathcal{H}_n \to \mathcal{H}_n$ be a firmly nonexpansive operator. Set $U = U_m \circ \dots \circ U_1$ and

$$ \theta_m = \|U\| + \sum_{k=1}^{m-1} \sum_{1 \le j_1 < \dots < j_k \le m-1} \|U_m \circ \dots \circ U_{j_k+1}\| \|U_{j_k} \circ \dots \circ U_{j_{k-1}+1}\| \dots \|U_{j_1} \circ \dots \circ U_1\|. \quad (25) $$

Let $S = Q_m \circ U_m \circ \dots \circ Q_1 \circ U_1$. Then the following hold:

1. $\theta_m/2^{m-1}$ is a Lipschitz constant of $S$.

2. Let $\alpha \in [1/2, 1]$. If $\mathcal{H}_m = \mathcal{H}_0$ and

$$ \|U - 2^m(1 - \alpha)\mathbf{1}\| - \|U\| + 2\theta_m \le 2^m \alpha, \quad (26) $$

then $S$ is $\alpha$-averaged.
In light of these results, we will now analyze the properties of the virtual network (12) based on the spectral quantities defined at the end of Section 4.1, and the parameters $(\lambda_n)_{1\le n\le m}$, $\bar{\chi}$ (and possibly $(\tau_n)_{1\le n\le m}$). One of the main difficulties with respect to the case already studied by [34] is that, in our context, the involved operators $(U_n)_{1\le n\le m}$ are no longer self-adjoint.

A preliminary result will be needed:

**Lemma 4.1.** Let $m \in \mathbb{N} \setminus \{0\}$ be the total number of layers. For every layer indices $n \in \{1, \dots, m\}$ and $i \in \{1, \dots, n\}$, the norm of $U_n \circ \dots \circ U_i$ is equal to $\sqrt{a_{i,n}}$ with

$$a_{i,n} = \frac{1}{2} \sup_p \left( \beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2 + \sqrt{(\beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2)^2 - 4\beta_{i,n,p}^2 \eta_{i,n}^2} \right), \quad (27)$$

where $p$ indexes the eigenspaces of $T^*T$.

*Proof.* Thanks to expressions (12), (18), (19), and (20), we can calculate the norm of $U_n \circ \dots \circ U_i$. For every $z = (x, b)$, $U_n \circ \dots \circ U_i z = (W_{i,n}x + \tilde{W}_{i,n}b, \eta_{i,n}b)$ and

$$\begin{aligned} \|U_n \circ \dots \circ U_i z\|^2 &= \|W_{i,n}x + \tilde{W}_{i,n}b\|^2 + \eta_{i,n}^2 \|b\|^2 \\ &= \|W_{i,n}x\|^2 + 2\langle W_{i,n}x, \tilde{W}_{i,n}b \rangle + \|\tilde{W}_{i,n}b\|^2 + \eta_{i,n}^2 \|b\|^2. \end{aligned}$$

Let $(\beta_{i,n,p}, v_p)_{p \in \mathbb{N}}$ defined by (22) and $(\tilde{\beta}_{i,n,p}, v_p)_{p \in \mathbb{N}}$ defined by (23) be the respective eigensystems of $W_{i,n}$ and $\tilde{W}_{i,n}$. Let us decompose $(x, b)$ on the basis of eigenvectors $(v_p)_p$ of $T^*T$, as

$$\left\{ \begin{array}{l} x = \sum_p \xi_p v_p, \\ b = \sum_p \zeta_p v_p. \end{array} \right. \qquad (28)$$

In the following, we will assume that $\mathcal{X}$ is a real Hilbert space. A similar proof can be followed for complex Hilbert spaces. We have then

$$\|U_n \circ \dots \circ U_i z\|^2 = \sum_p \beta_{i,n,p}^2 \xi_p^2 + 2 \sum_p \beta_{i,n,p} \tilde{\beta}_{i,n,p} \xi_p \zeta_p + \sum_p (\tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2) \zeta_p^2.$$

By definition of the operator norm,

$$\|U_n \circ \dots \circ U_i\|^2 = \sup_{\|x\|^2 + \|b\|^2 = 1} \left( \sum_p \beta_{i,n,p}^2 \xi_p^2 + (\eta_{i,n}^2 + \tilde{\beta}_{i,n,p}^2) \zeta_p^2 + 2\beta_{i,n,p} \tilde{\beta}_{i,n,p} \xi_p \zeta_p \right).$$

Note that, for every integer $p$ and $\omega_p = (\xi_p, \zeta_p) \in \mathbb{R}^2$,

$$\beta_{i,n,p}^2 \xi_p^2 + (\eta_{i,n}^2 + \tilde{\beta}_{i,n,p}^2) \zeta_p^2 + 2\beta_{i,n,p} \tilde{\beta}_{i,n,p} \xi_p \zeta_p = \omega_p^\top A_{i,n,p} \omega_p \quad (29)$$
where

$$A_{i,n,p} = \begin{pmatrix} \beta_{i,n,p}^2 & \beta_{i,n,p} \tilde{\beta}_{i,n,p} \\ \tilde{\beta}_{i,n,p} \beta_{i,n,p} & \eta_{i,n}^2 + \tilde{\beta}_{i,n,p}^2 \end{pmatrix}.$$

Hence,

$$\|U_n \circ \cdots \circ U_i\|^2 = \sup_{z=(\omega_p)_p, \|z\|=1} \sum_p \omega_p^\top A_{i,n,p} \omega_p. \quad (30)$$

For every index $p$, let $\nu_{i,n,p}$ be the maximum eigenvalue of the symmetric positive semidefinite matrix $A_{i,n,p}$. We have

$$\sum_p \omega_p^\top A_{i,n,p} \omega_p \le \sum_p \nu_{i,n,p} \|\omega_p\|^2 \le \left( \sup_p \nu_{i,n,p} \right) \sum_p \|\omega_p\|^2 = \left( \sup_p \nu_{i,n,p} \right) \|z\|^2.$$

Altogether with (30), this shows that

$$\|U_n \circ \cdots \circ U_i\|^2 \le \sup_p \nu_{i,n,p}.$$

In addition, from the definition of the supremum, for every $\epsilon > 0$, there exists $p^*$ such that

$$\sup_p \nu_{i,n,p} - \epsilon < \nu_{i,n,p^*}.$$

By setting $z^* = (\omega_p^*)_p$ where $\omega_p^* = 0$ if $p \ne p^*$, and $\omega_{p^*}^*$ is a unit norm eigenvector associated with eigenvalue $\nu_{i,n,p^*}$,

$$\sum_p \omega_p^\top A_{i,n,p} \omega_p = \omega_{p^*}^\top A_{i,n,p^*} \omega_{p^*} = \nu_{i,n,p^*} > \sup_p \nu_{i,n,p} - \epsilon.$$

Since $\epsilon$ can be chosen arbitrarily small, we deduce that

$$\|U_n \circ \cdots \circ U_i\|^2 = \sup_p \nu_{i,n,p} \quad (31)$$

We finally note that the two eigenvalues of matrix $A_{i,n,p}$ are the roots of the characteristic polynomial

$$(\forall \nu \in \mathbb{R}) \det(A_{i,n,p} - \nu \mathbb{1}_2) = (\beta_{i,n,p}^2 - \nu)(\tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2 - \nu) - \beta_{i,n,p}^2 \tilde{\beta}_{i,n,p}^2 \\ = \nu^2 - (\beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2)\nu + \beta_{i,n,p}^2 \eta_{i,n}^2.$$

The discriminant of this second-order polynomial reads

$$\begin{aligned} \Delta_{i,n,p} &= (\beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2)^2 - 4\beta_{i,n,p}^2 \eta_{i,n}^2 \\ &= (\beta_{i,n,p}^2 - \tilde{\beta}_{i,n,p}^2 - \eta_{i,n}^2)^2 + 4\beta_{i,n,p}^2 \tilde{\beta}_{i,n,p}^2 \ge 0. \end{aligned}$$
Therefore, for every integer p,

$$ \nu_{i,n,p} = \frac{1}{2} \left( \beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2 + \sqrt{(\beta_{i,n,p}^2 + \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2)^2 - 4\beta_{i,n,p}^2\eta_{i,n}^2} \right) . \quad (32) $$

It then follows from (31) that

$$ \|U_n \circ \dots \circ U_i\|^2 = a_{i,n}. \qquad \square $$

**Remark 4.1.** The previous result assumes that the product space $\mathcal{X} \times \mathcal{X}$ is equipped with the standard norm. Since we may be interested in the stability w.r.t. variations of the observed data y, and $b_0 = T^*y$, it might be more insightful to consider instead the following weighted norm:

$$ (\forall z = (x, b) \in \mathcal{X} \times \mathcal{X}) \quad \|z\| = \sqrt{\|x\|^2 + \|b\|^2_{\bar{\omega}}}, \quad (33) $$

where b is decomposed as in (28),

$$ \|b\|_{\bar{\omega}}^2 = \sum_p \frac{|\zeta_p|^2}{\bar{\omega}_p}, \quad (34) $$

and $\bar{\omega} = (\bar{\omega}_p)_p$ is such that $\inf_p \bar{\omega}_p > 0$ and $\sup_p \bar{\omega}_p < +\infty$. The latter weights aim to compensate the effect of $T^*$ on the observed data. This space renormalization is possible since the operator $Q_n$ remains firmly nonexpansive after this norm change, because of its specific structure.

As shown in Appendix A, the expression of $a_{i,n}$ for $n \in \{1, \dots, m\}$ and $i \in \{1, \dots, n\}$ in Lemma 4.1 is then modified as follows:

$$ a_{i,n} = \frac{1}{2} \sup_p \left( \beta_{i,n,p}^2 + \bar{\omega}_p \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2 + \sqrt{(\beta_{i,n,p}^2 + \bar{\omega}_p \tilde{\beta}_{i,n,p}^2 + \eta_{i,n}^2)^2 - 4\beta_{i,n,p}^2 \eta_{i,n}^2} \right). \quad (35) $$

We will now quantify the Lipschitz regularity of the network.

**Proposition 4.2.** Let $m \in \mathbb{N} \setminus \{0\}$. For every $n \in \{1, \dots, m\}$ and $i \in \{1, \dots, n\}$, let $a_{i,n}$ be given by (27). Set $\theta_0 = 1$ and define $(\theta_n)_{1 \le n \le m}$ recursively by

$$ (\forall n \in \{1, \dots, m\}) \quad \theta_n = \sum_{i=1}^{n} \theta_{i-1} \sqrt{a_{i,n}}. $$

Then $\theta_m/2^{m-1}$ is a Lipschitz constant of the virtual network (12). In addition,

$$ \sqrt{a_{1,m}} \le \frac{\theta_m}{2^{m-1}} \le \left( \prod_{n=1}^{m} a_{n,n} \right)^{1/2}. \quad (36) $$
Proof. According to Proposition 4.11, if $\theta_m$ is given by (25), then $\theta_m/2^{m-1}$ is a Lipschitz constant of the virtual network (12). On the other hand, it follows from [40, Lemma 3.3] that $\theta_m$ can be calculated recursively as

$$ (\forall n \in \{1, \dots, m\}) \quad \theta_n = \sum_{i=1}^{n} \theta_{i-1} \|U_n \circ \dots \circ U_i\|, $$

with $\theta_0 = 1$. Finally, Lemma 4.1 allows us to substitute $(\sqrt{a_{i,n}})_{1 \le i \le n}$ for $(\|U_n \circ \dots \circ U_i\|)_{1 \le i \le n}$ in the above expression. In addition, according to [41, Proposition 4.3(i)],

$$ \|U_m \circ \dots \circ U_1\| \le \frac{\theta_m}{2^{m-1}} \le \prod_{n=1}^{m} \|U_n\|. $$

It follows from Lemma 4.1, that $\|U_m \circ \dots \circ U_1\| = \sqrt{a_{1,m}}$ and

$$ (\forall n \in \{1, \dots, m\}) \quad \|U_n\| = \sqrt{a_{n,n}}, $$

which allows us to deduce inequality (36). □

**Remark 4.2.** Proposition 4.2 shows that the complexity for computing the proposed Lipschitz bound is quadratic as a function of the number of layers $m$ (more precisely, $O(m(m+1)/2)$).

We can relate the bounds provided on the Lipschitz constant in the previous proposition to simpler ones.

**Corollary 4.1.** Assume that $\mathcal{X} \times \mathcal{X}$ is equipped with the norm defined by (33). Then

$$ \sup_p \max \left\{ \left| \prod_{n=1}^{m} \beta_p^{(n)} \right|, \sqrt{\varpi_p \tilde{\beta}_{1,m,p}^2 + \eta_{1,m}^2} \right\} \le \frac{\theta_m}{2^{m-1}} \\ \le \prod_{n=1}^{m} \left( \sup_p \left( (\beta_p^{(n)})^2 + \varpi_p \left( \frac{\lambda_n}{1 + \lambda_n \chi_n} \right)^2 + \eta_n^2 \right) \right)^{1/2} . \quad (37) $$

Proof. For every $n \in \{1, \dots, m\}$ and integer $p$,

$$ \sqrt{(\beta_{n,n,p}^2 + \varpi_p \tilde{\beta}_{n,n,p}^2 + \eta_{n,n}^2)^2 - 4\beta_{n,n,p}^2 \eta_{n,n}^2} \le \beta_{n,n,p}^2 + \varpi_p \tilde{\beta}_{n,n,p}^2 + \eta_{n,n}^2. $$

We deduce from (35) that

$$ \begin{aligned} a_{n,n} &\le \sup_p \beta_{n,n,p}^2 + \varpi_p \tilde{\beta}_{n,n,p}^2 + \eta_{n,n}^2 \\ &= \sup_p (\beta_p^{(n)})^2 + \varpi_p \left( \frac{\lambda_n}{1 + \lambda_n \chi_n} \right)^2 + \eta_n^2, \end{aligned} $$
where the last equality follows from (20) and (24). The majorization in (37) is then a consequence of the upper bound in (36).

On the other hand, by rewriting (35) in a slightly different way,

$$ 
\begin{aligned}
a_{1,m} &= \frac{1}{2} \sup_p \left( \beta_{1,n,p}^2 + \varpi_p \tilde{\beta}_{1,m,p}^2 + \eta_{1,m}^2 \right. \\
&\qquad \left. + \sqrt{(\beta_{1,m,p}^2 - \varpi_p \tilde{\beta}_{1,m,p}^2 - \eta_{1,m}^2)^2 + 4\varpi_p \beta_{1,m,p}^2 \tilde{\beta}_{1,m,p}^2} \right), \\
&\ge \frac{1}{2} \sup_p \left( \beta_{1,n,p}^2 + \varpi_p \tilde{\beta}_{1,m,p}^2 + \eta_{1,m}^2 + |\beta_{1,m,p}^2 - \varpi_p \tilde{\beta}_{1,m,p}^2 - \eta_{1,m}^2| \right) \\
&= \sup_p \max \left\{ \beta_{1,n,p}^2, \varpi_p \tilde{\beta}_{1,m,p}^2 + \eta_{1,m}^2 \right\}.
\end{aligned}
\quad \square $$

**Remark 4.3.** As shown in [34], if we are just interested in the sensitivity w.r.t. $x_0$ of the unfolded algorithm by assuming that $b_0$ is unperturbed, formulas in Proposition 4.2 remain valid by setting

$$ (\forall n \in \{1, \dots, m\})(\forall i \in \{1, \dots, n\}) \quad a_{i,n} = \sup_p |\beta_{i,n,p}|. $$

This can be viewed as a limit case of (35) where $\sup_p \varpi_p \to 0$ and, for every $n \in \{1, \dots, m\}$, $\eta_n = 0$.

We will next provide conditions ensuring that the virtual network is an averaged operator.

**Proposition 4.3.** Assume that $\mathcal{X} \times \mathcal{X}$ is equipped with the norm defined by (33). Let $m \in \mathbb{N} \setminus \{0\}$. Let $a_{1,m}$ be defined by (35) and $\theta_m$ be defined in Proposition 4.2. Let $\alpha \in [1/2, 1]$. Define

$$ b_{\alpha} = \frac{1}{2} \sup_p \left( (\beta_{1,m,p} - \gamma_{\alpha})^2 + (\eta_{1,m} - \gamma_{\alpha})^2 + \varpi_p \tilde{\beta}_{1,m,p}^2 \right. \\
\left. + \sqrt{((\beta_{1,m,p} - \gamma_{\alpha})^2 + (\eta_{1,m} - \gamma_{\alpha})^2 + \varpi_p \tilde{\beta}_{1,m,p}^2)^2 - 4(\beta_{1,m,p} - \gamma_{\alpha})^2 (\eta_{1,m} - \gamma_{\alpha})^2} \right), \quad (38) $$

with $\gamma_{\alpha} = 2^m(1 - \alpha)$. Then virtual network (12) is $\alpha$-averaged if

$$ \sqrt{b_{\alpha}} - \sqrt{a_{1,m}} \le 2^m \alpha - 2\theta_m. \quad (39) $$

*Proof.* Let us calculate the operator norms of $U$ and $U - \gamma_{\alpha}1$, where $U$ is given by (17).

**Norm of U.**

Applying Lemma 4.1 when $i = 1$ and $n = m$ yields

$$ \|U\|^2 = a_{1,m}. $$
Norm of $U - \gamma_\alpha \mathbf{1}$.

We follow the same reasoning as in the proof of Lemma 4.1 (see also Appendix A). We have

$$ \|U - \gamma_{\alpha} \mathbf{1}\|^2 = \sup_{z=(\omega_p)_p, \|z\|=1} \sum_p \omega_p^\top B_p \omega_p, $$

where $B_p$ is the symmetric positive semidefinite matrix given by

$$ B_p = \begin{pmatrix} (\beta_{1,m,p} - \gamma_\alpha)^2 & \sqrt{\overline{\omega_p}}(\beta_{1,m,p} - \gamma_\alpha)\tilde{\beta}_{1,m,p} \\ \sqrt{\overline{\omega_p}}(\beta_{1,m,p} - \gamma_\alpha)\tilde{\beta}_{1,m,p} & (\eta_{1,m} - \gamma_\alpha)^2 + \overline{\omega_p}\tilde{\beta}_{1,m,p}^2 \end{pmatrix}. $$

By definition of the spectral norm,

$$ \|U - 2^m(1 - \alpha)\mathbf{1}\|^2 = \sup_p \nu_p, \qquad (40) $$

where, for every integer $p$, $\nu_p$ is the maximum eigenvalue of $B_p$. The two eigenvalues of this matrix are the roots of the polynomial

$$ (\forall \nu \in \mathbb{R}) \quad \det(B_p - \nu \mathbf{1}_2) = \nu^2 - ((\beta_{1,m,p} - \gamma_\alpha)^2 + (\eta_{1,m} - \gamma_\alpha)^2 + \overline{\omega_p}\tilde{\beta}_{1,m,p}^2) \nu + (\beta_{1,m,p} - \gamma_\alpha)^2(\eta_{1,m} - \gamma_\alpha)^2. \qquad (41) $$

Solving the corresponding second-order equation leads to

$$ \nu_p = \frac{1}{2} \left( (\beta_{1,m,p} - \gamma_\alpha)^2 + (\eta_{1,m} - \gamma_\alpha)^2 + \overline{\omega_p} \tilde{\beta}_{1,m,p}^2 + \sqrt{((\beta_{1,m,p} - \gamma_\alpha)^2 + (\eta_{1,m} - \gamma_\alpha)^2 + \overline{\omega_p} \tilde{\beta}_{1,m,p}^2)^2 - 4(\beta_{1,m,p} - \gamma_\alpha)^2(\eta_{1,m} - \gamma_\alpha)^2} \right). $$

Then, it follows from (40) that $\|U - \gamma_\alpha \mathbf{1}\|^2 = b_\alpha$.

Conclusion of the proof.

Based on the previous calculations, Condition (39) is equivalent to (26). In addition, note that, as shown in Appendix A, for every $n \in \{1, \dots, m\}$, $Q_n$ in (12) is firmly nonexpansive. By applying now Proposition 4.12, we deduce that, when Condition (39) holds, virtual network (12) is $\alpha$-averaged. $\square$

**Remark 4.4.** Condition (39) just provides a sufficient condition for the averagedness of virtual network (12).

We conclude this subsection by a result emphasizing the interest of the leakage factors:

**Proposition 4.4.** Let $S$ be the virtual Model (12) without leakage factor, i.e., for every $n \in \{1, \dots, m\}$, $\eta_n = 1$. Assume that there exists $(x, b, b') \in \mathcal{X}^3$ such that the output $x_m$ of the unfolded algorithm associated with $(x, b)$ is distinct from its output $x'_m$ associated with $(x, b')$. Then the Lipschitz constant of $S$ is greater than 1.
Proof. If $\eta_{1,m} = 1$ then, for every $n \in \{1, \dots, m\}$,

$$ z_n = \begin{pmatrix} x_n \\ b_n \end{pmatrix}, \quad z_n = Q_n(U_n z_{n-1}), \quad \text{with} \quad \left\{ \begin{array}{l} Q_n = \begin{pmatrix} R_n \\ \mathbb{1} \end{pmatrix}, \\ U_n = \begin{pmatrix} W_n & V_n \\ 0 & \mathbb{1} \end{pmatrix}. \end{array} \right. \qquad (42) $$

Suppose that the virtual network $S$ is nonexpansive. Thus, for every $(x, b) \in \mathcal{X}^2$ and $(x', b') \in \mathcal{X}^2$, $\|S(x, b) - S(x', b')\| \le \|x - x'\|^2 + \|b - b'\|^2_{\bar{\omega}}$. Let $x_m$ (resp. $x'_m$) denote the first component in $\mathcal{X}$ of $S(x, b)$ (resp. $S(x', b')$). Since

$$ \begin{aligned} \|S(x, b) - S(x, b')\|^2 &= \left\| \begin{pmatrix} x_m - x'_m \\ b - b' \end{pmatrix} \right\|^2 \\ &= \|x_m - x'_m\|^2 + \|b - b'\|^2_{\bar{\omega}} \\ &\le \|x - x'\|^2 + \|b - b'\|^2_{\bar{\omega}}, \end{aligned} $$

we deduce that

$$ \|x_m - x'_m\| \le \|x - x'\|. $$

Consequently, when $x = x'$, we have $x_m = x'_m$ whatever the choice of $b \ne b'$. This contradicts our standing assumption. $\square$

Note that the assumption made on $S$ is relevant from the applicative viewpoint. It is typically satisfied if we have two different ground truth data $(\bar{x}, \bar{x}') \in \mathcal{X}^2$ and we observe $y = T\bar{x} + w$ and $y' = T\bar{x}' + w$ with $b = T^*y \ne T^*y' = b'$, while the same initialization $x$ is used in the iterative recovery process. Then, the condition $x_m \ne x'_m$ just means that the unfolded network is sensitive to the observed data, as can be expected from a model trained to perform suitable restoration.

The above result shows that the virtual network without leakage factors cannot be averaged since any averaged operator is nonexpansive.

## 4.3 Link with the original neural network – direct approach

In this subsection we go back to our initial model defined by (13). For simplicity, we assume that, for every $p, \bar{\omega}_p = 1$. We consider two distinct inputs $z = (x, b)$ and $z' = (x', b')$ in $\mathcal{X} \times \mathcal{X}$. Let $z_m = (x_m, b_m)$ and $z'_m = (x'_m, b'_m)$ be the outputs of the $m$-th layer of virtual network (12) associated with inputs $z$ and $z'$, respectively. Then,

$$ \begin{aligned} \|z'_m - z_m\|^2 &= \|x'_m - x_m\|^2 + \|b'_m - b_m\|^2 \\ &= \|x'_m - x_m\|^2 + \eta_{1,m}^2 \|b' - b\|^2, \end{aligned} $$

and, by applying Proposition 4.2,

$$ \|z'_m - z_m\|^2 \le \frac{\theta_m^2}{2^2(m-1)} (\|x' - x\|^2 + \|b' - b\|^2). $$
Then, the following inequality allows us to quantify the Lipschitz properties of the neural network (6) with an error on $b_0$:

$$ \|x'_{m} - x_{m}\|^{2} \leq \frac{\theta_{m}^{2}}{2^{2(m-1)}} \|x' - x\|^{2} + \left( \frac{\theta_{m}^{2}}{2^{2(m-1)}} - \eta_{1,m}^{2} \right) \|b' - b\|^{2}. $$

Two particular cases are of interest:

* If the network is initialized with a fixed signal, say $x_0 = x'_0 = 0$, then

$$ \|x'_{m} - x_{m}\|^{2} \leq \left( \frac{\theta_{m}^{2}}{2^{2(m-1)}} - \eta_{1,m}^{2} \right) \|b' - b\|^{2}. $$

So, a Lipschitz constant with respect to the input data $T^*y$ is

$$ \vartheta_m = \sqrt{\frac{\theta_m^2}{2^{2(m-1)}} - \eta_{1,m}^2}. \qquad (43) $$

This shows, in particular, that the Lipschitz constant of the virtual network $\theta_m/2^{m-1}$ cannot be lower than $\eta_{1,m}$. This result is consistent with the observation made at the end of Section 4.2.

* On the other hand, if the initialization is dependent on the observed image, i.e. $x_0 = b$ and $x'_0 = b'$,

$$ \|x'_{m} - x_{m}\|^{2} \leq \left( \frac{\theta_{m}^{2}}{2^{2m-3}} - \eta_{1,m}^{2} \right) \|b' - b\|^{2}. $$

So a higher Lipschitz constant value w.r.t. to the input data is obtained:

$$ \vartheta_m = \sqrt{\frac{\theta_m^2}{2^{2m-3}} - \eta_{1,m}^2}. \qquad (44) $$

## 4.4 Use of a semi-norm

Proposition 4.4 suggests that we may need a finer strategy to evaluate the nonex-pansiveness properties of Model (6). On the product space $\mathcal{X} \times \mathcal{X}$, we define the semi-norm which takes only into account the first component of the vectors:

$$ z = (x, b) \mapsto |z| = \|x\|. \qquad (45) $$

Let $L: \mathcal{X} \times \mathcal{X} \to \mathcal{X} \times \mathcal{X}$ be any bounded linear operator and, for every $z \in \mathcal{X} \times \mathcal{X}$, let $Lz = ((Lz)_x, (Lz)_b)$. We define the associated operator semi-norm

$$ |L| = \sup_{\|z\|=1} \|(Lz)_x\|. \qquad (46) $$

The following preliminary result will be useful subsequently.
**Lemma 4.2.** Assume that $\mathcal{X} \times \mathcal{X}$ is equipped with the norm defined by (33). Let $m \in \mathbb{N} \setminus \{0\}$. For every $n \in \{1, \dots, m\}$ and $i \in \{1, \dots, n\}$, the seminorm $|U_n \circ \dots \circ U_i|$ is equal to $\sqrt{\bar{a}_{i,n}}$ with

$$ \bar{a}_{i,n} = \sup_p \left( \beta_{i,n,p}^2 + \varpi_p \tilde{\beta}_{i,n,p}^2 \right) \le a_{i,n}. \qquad (47) $$

*Proof.* The inequality $\bar{a}_{i,n} \le a_{i,n}$ follows from the fact that $|U_n \circ \dots \circ U_i| \le \|U_n \circ \dots \circ U_i\|$.

The seminorm of $U_n \circ \dots \circ U_i$ is the same as the norm of $U_n \circ \dots \circ U_i$ where $\eta_n$ has been set to 0. The result thus follows from Lemma 4.1 where $\eta_{i,n} = 0$. □

### 4.4.1 Network with input $(x_0, b_0)$ and output $x_m$

The following result for the quantification of the Lipschitz regularity is an offspring of Proposition 4.2.

**Proposition 4.5.** Let $m \in \mathbb{N} \setminus \{0\}$. For every $i \in \{1, \dots, n\}$ and $n \in \{1, \dots, m-1\}$, let $a_{i,n}$ be defined by (27) and let $\bar{a}_{i,m}$ be given by (47). Set $\theta_0 = 1$ and define

$$ (\forall n \in \{1, \dots, m-1\}) \quad \theta_n = \sum_{i=1}^{n} \theta_{i-1} \sqrt{\bar{a}_{i,n}}, \qquad (48) $$

$$ \bar{\theta}_m = \sum_{i=1}^{m} \theta_{i-1} \sqrt{\bar{a}_{i,m}}. \qquad (49) $$

Then the network (13) with input $(x_0, b_0)$ and output $x_m$ is $\bar{\theta}_m/2^{m-1}$-Lipschitz. In addition,

$$ \sqrt{\bar{a}_{1,m}} \le \frac{\bar{\theta}_m}{2^{m-1}} \le \sqrt{\bar{a}_{m,m}} \left( \prod_{n=1}^{m-1} a_{n,n} \right)^{1/2}. \qquad (50) $$

*Proof.* Network (13) can be expressed as $R_m \circ \bar{U}_m \circ Q_{m-1} \circ U_{m-1} \circ \dots \circ Q_1 \circ U_1$ where

$$ \bar{U}_m = D_x \circ U_m \qquad (51) $$

and $D_x$ is the decimation operator

$$ D_x = [1 \ 0]. \qquad (52) $$

This network has the same Lipschitz properties as the network in (12) with $\eta_m = 0$. The result can thus be deduced from Proposition 4.2 by setting $\eta_{1,m} = 0$. □
Note that the resulting estimate of the Lipschitz constant is smaller than the one obtained for the VNN. Indeed, according to Lemma 4.2,

$$ \bar{\theta}_m = \sum_{i=1}^{m} \theta_{i-1} \sqrt{\bar{a}_{i,m}} \le \sum_{i=1}^{m} \theta_{i-1} \sqrt{\bar{a}_{i,m}} = \theta_m. \quad (53) $$

### 4.4.2 Network with a single input and a single output

Another possibility for investigating Lipschitz properties, including averagedness, consists of defining a network from $\mathcal{X}$ to $\mathcal{X}$. We will focus on specific networks of the form

$$ R_m \circ \bar{U}_m \circ Q_{m-1} \circ U_{m-1} \cdots \circ Q_1 \circ \hat{U}_1, \quad (54) $$

where $\bar{U}_m$ is given by (51) and

$$ \hat{U}_1 = U_1 \begin{bmatrix} F \\ 1 \end{bmatrix}. \quad (55) $$

Hereabove, $F$ is a bounded linear operator from $\mathcal{X}$ to $\mathcal{X}$. We will make the following additional assumption:

Operator $F$ can be diagonalized in basis $(v_p)_p$, its eigenvalues being denoted by $(\phi_p)_p$.

This assumption encompasses the three following cases:

1. The first one assumes that $x_0 = 0$ in (6), by setting $F = 0$.

2. The second one assumes that $x_0 = b_0$ in (6), by setting $F = 1$ and, for every $p$, $\phi_p = 1$.

3. The third one assumes that $x_0 = Fb_0 = FT^*y$ where $F$ is a pre-filter, while $T$ and $(D_n)_{1 \le n \le m}$ also are filtering operations. For instance, $F$ could perform a rough restoration of the sought signal from the observed one $b_0$ [47], possibly by ensuring that $FT^*$ is a Wiener filter. In this case, $(\phi_p)_p$ corresponds to the complex-valued frequency response of this pre-filter at discrete frequencies indexed by $p$.

Similarly to Lemma 4.1, we obtain the following result.

**Lemma 4.3.** Assume that $\mathcal{X} \times \mathcal{X}'$ is equipped with the norm defined by (33) and the input space is $(\mathcal{X}', \|\cdot\|_{\varpi})$. Let $m \in \mathbb{N} \setminus \{0, 1\}$. For every $n \in \{1, \dots, m-1\}$, the norm of $U_n \circ \dots \circ U_2 \circ \hat{U}_1$ is equal to $\sqrt{\hat{a}_{1,n}}$ with

$$ \hat{a}_{1,n} = \sup_p \left( \varpi_p |\beta_{1,n,p} \phi_p + \tilde{\beta}_{1,n,p} |^2 \right) + \eta_{1,n}^2 \quad (56) $$

and the norm of $\bar{U}_m \circ U_{m-1} \circ \dots \circ U_2 \circ \hat{U}_1$ is equal to $\sqrt{\hat{a}_{1,m}}$ with

$$ \hat{a}_{1,m} = \sup_p \varpi_p |\beta_{1,m,p} \phi_p + \tilde{\beta}_{1,m,p} |^2. \quad (57) $$
Proof. For every $n \in \{1, \dots, m-1\}$,

$$
\begin{align}
\|U_n \circ \dots \circ U_2 \circ \hat{U}_1 b_0\|^2 &= \|W_{1,n} F b_0 + \tilde{W}_{1,n} b_0\|^2 + \eta_{1,n}^2 \|b_0\|_{\bar{\omega}}^2 \nonumber \\
&= \sum_p (|\beta_{1,n,p} \phi_p + \tilde{\beta}_{1,n,p}|^2 + \bar{\omega}_p^{-1} \eta_{1,n}^2) |\zeta_p|^2 \tag{58} \\
&\le \sup_p (\bar{\omega}_p |\beta_{1,n,p} \phi_p + \tilde{\beta}_{1,n,p}|^2 + \eta_{1,n}^2) \|b_0\|_{\bar{\omega}}^2, \nonumber
\end{align}
$$

where $b_0$ has been decomposed as

$$b_0 = \sum_p \zeta_p v_p.$$

It can be deduced that

$$\|U_n \circ \dots \circ U_2 \circ \hat{U}_1\|^2 = \hat{a}_{1,n}.$$

Equation (58) remains valid when $n = m$, by setting $\eta_m = 0$, i.e. $\eta_{1,m} = 0$. This leads to the expression of $\hat{a}_{1,m}$. $\square$

By proceeding similarly to the proof of Proposition 4.3, we deduce the following stability result.

**Proposition 4.6.** Let $m \in \mathbb{N} \setminus \{0, 1\}$. For every $i \in \{2, \dots, n\}$ and $n \in \{1, \dots, m-1\}$, let $a_{i,n}$ be defined by (27) and let $\bar{a}_{i,m}$ be given by (47). For every $n \in \{1, \dots, m\}$, let $\hat{a}_{1,n}$ be defined by (56) and (57). Define $(\hat{\theta}_n)_{1 \le n \le m}$ recursively by

$$(\forall n \in \{1, \dots, m-1\}) \quad \hat{\theta}_n = \sqrt{\hat{a}_{1,n}} + \sum_{i=2}^{n} \hat{\theta}_{i-1} \sqrt{a_{i,n}}, \quad (59)$$

$$\hat{\theta}_m = \sqrt{\hat{a}_{1,m}} + \sum_{i=2}^{m} \hat{\theta}_{i-1} \sqrt{\bar{a}_{i,m}}, \quad \hat{\theta}_0 = 1. \quad (60)$$

Then network (54) is $\hat{\theta}_m/2^{m-1}$-Lipschitz. In addition,

$$\sqrt{\hat{a}_{1,m}} \le \frac{\hat{\theta}_m}{2^{m-1}} \le \left( \bar{a}_{m,m} \left( \prod_{n=2}^{m-1} a_{n,n} \right) \hat{a}_{1,1} \right)^{1/2}. \quad (61)$$

As shown hereafter, the obtained Lipschitz constant is an increasing function of the pre-filter components $(\phi_p)_p$ when those are nonnegative and real-valued.

**Proposition 4.7.** Assume that

$$(\forall n \in \{1, \dots, m\}) \quad \sup_p \lambda_n(\beta_{T,p} + \beta_{D_n,p}) \le 1. \quad (62)$$

Let us denote by $\hat{\theta}_m(\phi)$ the Lipschitz parameter defined in the previous proposition for a given value of $\phi = (\phi_p)_p$. If $\phi = (\phi_p)_p$ and $\phi' = (\phi'_p)$ are two sequences of
nonnegative reals, then

$$ (\forall p) \phi_p \ge \phi'_p \Rightarrow \hat{\theta}_m(\phi) \ge \hat{\theta}_m(\phi'). $$

*Proof.* The real values of $(a_{i,n})_{1 \le i \le n, 1 \le n \le m-1}$, $(\bar{a}_{i,m})_{1 \le i \le n}$, and $(\hat{a}_{1,n})_{1 \le n \le m}$, defined respectively by (27), (47), and (56)-(57), are nonnegative. If Condition (62) holds, then the constants defined in (21), (22), and (23) are nonnegative. Thus, according to (56)-(57), for every $n \in \{1, \dots, m\}$, $\hat{a}_{1,n}$ is an increasing function of the components of the associated pre-filter $\phi$. It follows from (59) that the monotonicity property holds for $(\hat{\theta}_n)_{1 \le n \le m-1}$ and, from (60), that it extends to $\hat{\theta}_m$. $\square$

We deduce, in particular, that the Lipschitz constant is lower when initializing with $x_0 = 0$ (case i) rather than $x_0 = b_0$ (case ii)).

# 5 Numerical illustrations

## 5.1 Stationary case

We first consider the case when the parameters do not vary across the layers. In this stationary case, we have: for every $n \in \{1, \dots, m\}$, $D_n = D$, $\lambda_n = \lambda$, $\eta_n = \eta$, and $\chi_n = \chi$. Details concerning this scenario are provided in Appendix B.

We set $T$ as a two-dimensional $3 \times 3$ discrete uniform blur such that $\|T\| = 1$. In addition,

$$ D = \sqrt{\tau} \begin{bmatrix} \nabla_H \\ \nabla_V \end{bmatrix} $$

where $\nabla_H$ (resp. $\nabla_V$) is the horizontal (resp. vertical) discrete gradient operator and $\tau = 10^{-2}$. Since $\mathcal{X} = \mathbb{R}^{N_1 \times N_2}$ ($N_1 = N_2 = 256$), the basis $(v_p)_{1 \le p \le N_1 N_2}$ is associated with the 2D discrete Fourier transform. In our tests, we employ the weighted norm defined by

$$ (\forall p \in \{1, \dots, N_1 N_2\}) \quad \bar{\omega}_p = \beta_{T,p}^2 + \epsilon \quad (63) $$

with $\epsilon = 10^{-2}$. In the following experiments, the strong convexity modulus $\bar{\chi}$ is zero.

Figure 2 illustrates the computed values of $\theta_m/2^{m-1}$ for the considered VNN, based on Proposition 4.2, as functions of the stepsize $\lambda$ and the leakage factor $\eta$, for one layer ($m=1$) and fifteen layers ($m=15$). The left part highlights the regions where the Lipschitz constant is less than one, hence the VNN is averaged. However, we observe that achieving nonexpansiveness of the virtual network requires relatively restrictive parameter choices. Interestingly, small values of the leakage factor exhibit a stabilizing effect, as suggested by Proposition 4.4. Moreover, the region of nonexpansiveness is larger for $\theta_{15}/2^{14}$ compared to $\theta_1$, indicating that a multi-layer stability analysis yields more accurate results. This observation is further confirmed in Figure 3, which shows the difference between the separable estimate after 15 layers, given by $(\theta_1)^{15}$, and the computed value $\theta_{15}/2^{14}$.

Figure 4 illustrates a similar behavior for the Lipschitz constant $\bar{\theta}_{15}/2^{14}$, based on the use of a semi-norm on the VNN output, which has been computed by using
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Stationary case: Lipschitz constant of the VNN as a function of $\lambda$ and $\eta$. First row: $\theta_1$, second one: $\theta_{15}/2^{14}$. First column: numerical value, second one: white if $\theta_m/2^{m-1} \le 1$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Stationary case: Difference $(\theta_1)^{15} - \theta_{15}/2^{14}$ for the VNN in the stationary case as a function of $\lambda$ and $\eta$.</figcaption></figure>

**Proposition 4.5.** The bottom map displays the difference between $\theta_{15}/2^{14}$ and $\bar{\theta}_{15}/2^{14}$, which, as shown in (53), is always nonnegative.

Lastly, for the network with a single input and a single output, Figure 5 presents the Lipschitz constant $\hat{\theta}_{15}/2^{14}$, computed using Proposition 4.6, for two different operator choices, $F=0$ and $F=1$. We observe that selecting $F=0$ results in greater stability compared to $F=1$, which aligns with Proposition 4.7.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Stationary case: Lipschitz constant of the network with inputs $(x_0, b_0)$ and output $x_{15}$ as a function of $\lambda$ and $\eta$. Top left: $\tilde{\theta}_{15}/2^{14}$. Top right: white when $\theta_{15}/2^{14} < 1$. Bottom: Difference $\theta_{15} - \tilde{\theta}_{15}$.</figcaption></figure>

## 5.2 Nonstationary case

We now consider the scenario most commonly encountered in unrolled algorithms, where the parameter values vary across layers. In our experiments, we set

$$ (\forall n \in \{1, \dots, m\}) \quad D_n = \sqrt{\tau_n} \begin{bmatrix} \nabla_H \\ \nabla_V \end{bmatrix}, $$

and choose $T$ as a two-dimensional discrete 5 × 5 uniform blur such that $\|T\| = 1$. Parameters $\{\tau_n\}_{1 \le n \le m} \subset [0.099, 0.0249]$ and $\{\lambda_n\}_{1 \le n \le m} \subset [0.5128, 0.9585]$, typical of those learned from an image dataset, are used. Constant values for $(\eta_n)_{1 \le n \le m}$ and $(\chi_n)_{1 \le n \le m}$ have been chosen. They are set to 0.98 and $10^{-3}$, respectively, although setting the latter strong convexity moduli to zero results in only a slight reduction in the estimated Lipschitz regularity. The weighted norm is defined similarly to the previous example (see (63)).

Figure 6 illustrates the variations in the Lipschitz constant estimates $\theta_n/2^{n-1}$, $\bar{\theta}_n/2^{n-1}$, and $\hat{\theta}_n/2^{n-1}$ as a function of the layer index $n$. We also compare the proposed estimates with the lower and upper bounds derived in Propositions 4.2, 4.5, and 4.6. It can be observed that the proposed bounds are significantly tighter than the separable lower bounds. While the resulting values exceed 1, they remain within a satisfactory order of magnitude when compared to those encountered in sensitive networks [44].
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Stationary case: Lipschitz constant $\hat{\theta}_{15}/2^{14}$ as a function of $\lambda$ and $\eta$. First row: $F = 0$, second one: $F = 1$. First column: numerical value, second one: white if $\hat{\theta}_{15}/2^{14} \le 1$.</figcaption></figure>

# 6 Conclusion

Unrolling optimization algorithms into neural architectures offers significant advantages, including a reduced number of parameters to learn, efficient GPU implementation, and the ability to optimize the quality of the results using relevant loss functions. However, there is currently a lack of theoretical analysis for these structures.

In this work, by leveraging tools from fixed-point theory, we conducted an analysis of the robustness of an unrolled version of the forward-backward algorithm, extending the work in [34], where sensitivity to observed data was overlooked. Our analysis relies on Lipschitz bounds that are computationally efficient (with quadratic complexity) and tighter than the separable bounds commonly used in the neural network community [42].

We acknowledge that, although these bounds are rigorous, they are not optimal. Future work could focus on deriving tighter bounds to further refine the analysis. Additionally, exploring more general proximal algorithms beyond the one considered in this study would be an interesting direction. For example, we could consider the presence of a heavy ball/inertial term in the forward-backward algorithm. Finally, we believe the $\alpha$-averaging properties we established could inspire the design of recurrent networks in the spirit of [40].
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Nonstationary case. From top to bottom: $\theta_n/2^{n-1}$, $\bar{\theta}_n/2^{n-1}$, and $(\hat{\theta}_n/2)^{n-1}$ versus $n$ (in blue), and the corresponding upper bound (in orange) and lower bound (in red). For $(\hat{\theta}_n/2)^{n-1}$, $F = 1$.</figcaption></figure>

**Acknowledgements.** The work by J.-C. Pesquet was supported by the ANR Chair in Artificial Intelligence BRIDGEABLE. E. Chouzenoux acknowledges support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925.
# Appendix A Calculations related to Remark 4.1.

## Firm nonexpansiveness of $Q_n$.

For every $n \in \{1, \dots, m\}$, $z = (x, b) \in \mathcal{X}^2$, and $z' = (x', b') \in \mathcal{X}^2$,

$$ 
\begin{aligned}
& \|Q_n z - Q_n z'\|^2 + \|(1 - Q_n)z + (1 - Q_n)z'\|^2 \\
& = \|R_n x - R_n x'\|^2 + \|b - b'\|_{\bar{\omega}}^2 + \|(1 - R_n)x - (1 - R_n)x'\|^2 \\
& \le \|x - x'\|^2 + \|b - b'\|_{\bar{\omega}}^2 \\
& = \|z - z'\|^2,
\end{aligned}
$$

where the inequality follows from the firm nonexpansiveness of $R_n$.

## Proof of (35).

The introduction of the weighted norm in (33) induces the following changes in the proof of Lemma 4.1.

Using the same notation as in the proof, the norm of $U_n \circ \dots \circ U_i z$ is now given by

$$ \|U_n \circ \dots \circ U_i z\|^2 = \sum_p \left( \beta_{i,n,p}^2 \xi_p^2 + 2\beta_{i,n,p} \tilde{\beta}_{i,n,p} \xi_p \zeta_p + \left( \tilde{\beta}_{i,n,p}^2 + \frac{\eta_{i,n}^2}{\bar{\omega}_p} \right) \zeta_p^2 \right). $$

The above expression has to be maximized on $(x, b)$ subject to the constraint

$$ \|x\|^2 + \|b\|_{\bar{\omega}}^2 = \sum_p \left( \xi_p^2 + \frac{\zeta_p^2}{\bar{\omega}_p} \right) = 1. $$

This suggests to set, for every integer $p$, $\omega'_p = (\xi_p, \zeta_p / \sqrt{\bar{\omega}_p})$, which leads to

$$ \|U_n \circ \dots \circ U_i\|^2 = \sup_{z'=(\omega'_p)_p, \|z'\|=1} \sum_p (\omega'_p)^\top A'_{i,n,p} (\omega'_p), $$

where the supremum is calculated on the standard unit ball and

$$ A'_{i,n,p} = \left( \frac{\beta_{i,n,p}^2}{\sqrt{\bar{\omega}_p}} \frac{\sqrt{\bar{\omega}_p} \beta_{i,n,p} \tilde{\beta}_{i,n,p}}{\eta_{i,n}^2 + \bar{\omega}_p \tilde{\beta}_{i,n,p}^2} \right). $$

The rest of the derivation of (35) is similar to that in the proof of Lemma 4.1.

# Appendix B Stationary case.

This appendix shows that, in the stationary case, some of the quantities of interest can be easily expressed. This allows us to state more explicitly necessary/sufficient conditions for the VNN to be 1-Lipschitz.
In this case, for every $p$, $\beta_p^{(n)} = \beta_p$. In addition, $\eta_{1,m} = \eta^m$, $\beta_{1,m,p} = \beta_p^m$, and (23) simplifies to

$$ \tilde{\beta}_{1,m,p} = \frac{\lambda}{1 + \lambda\chi} \sum_{j=0}^{m-1} \beta_p^j \eta^{m-j-1}. $$

Then, it follows from (37) that necessary conditions for $\theta_m/2^{m-1} = 1$ are

$$ (\forall p) \quad \begin{cases} |\beta_p| \le 1 \\ \varpi_p \left( \frac{\lambda}{1 + \lambda\chi} \right)^2 \left( \sum_{j=0}^{m-1} \beta_p^j \eta^{m-j-1} \right)^2 + \eta^{2m} \le 1 \quad \Rightarrow \eta < 1. \end{cases} $$

The first condition is classical for the convergence of the proximal gradient algorithm and it is satisfied if and only if $\lambda(\sup_p(\beta_{T,p} + \beta_{D,p}) - \chi) \le 2$.

Furthermore, (37) yields the following sufficient condition for the VNN to be 1-Lipschitz:

$$ (\forall p) \quad \beta_p^2 + \varpi_p \left( \frac{\lambda}{1 + \lambda\chi} \right)^2 + \eta^2 \le 1, $$

which, by virtue of (21), is also equivalent to

$$ (\forall p) \quad (1 - \lambda(\beta_{T,p} + \beta_{D,p}))^2 + \varpi_p \lambda^2 + (1 + \lambda\chi)^2 (\eta^2 - 1) \le 0. \quad (B1) $$

As expected, the strong convexity of the regularizing function $g$ allows this condition to be more easily satisfied since $\eta^2 - 1 \le 0$. It can be observed that a small value of $\eta$ also has a favorable effect. For simplicity, let us focus on the case when $\chi = 0$, which provides the following sufficient condition for (B1):

$$ (\forall p) \quad ((\beta_{T,p} + \beta_{D,p})^2 + \varpi_p)\lambda^2 - 2(\beta_{T,p} + \beta_{D,p})\lambda + \eta^2 \le 0. $$

These inequalities can be satisfied only if

$$ (\forall p) \quad \begin{cases} (\beta_{T,p} + \beta_{D,p})^2 - \eta^2((\beta_{T,p} + \beta_{D,p})^2 + \varpi_p) \ge 0 \\ \lambda_{1,p}(\eta) \le \lambda \le \lambda_{2,p}(\eta), \end{cases} $$

where $\lambda_{1,p}(\eta)$ and $\lambda_{2,p}(\eta)$ are the smallest and largest roots of the quadratic polynomial $((\beta_{T,p} + \beta_{D,p})^2 + \varpi_p)\lambda^2 - 2(\beta_{T,p} + \beta_{D,p})\lambda + \eta^2$. The first inequalities are satisfied if and only if

$$ \eta \le \inf_p \frac{\beta_{T,p} + \beta_{D,p}}{\sqrt{(\beta_{T,p} + \beta_{D,p})^2 + \varpi_p}} $$

while the second require that $\sup_p \lambda_{1,p}(\eta) \le \inf_p \lambda_{2,p}(\eta)$.

## References

[1] Bertero, M., Boccacci, P., Desiderà, G., Vicidomini, G.: Image blurring with Poisson data: from cells to galaxies. Inverse Problems 25(12), 123006 (2009)
[2] Chouzenoux, E., Moussaoui, S., Idier, J.: Majorize-minimize linesearch for inversion methods involving barrier function optimization. Inverse Problems 28(6), 065011 (2012)

[3] Wu, Z., Sun, Y., Matlock, A., Liu, J., Tian, L., Kamilov, U.S.: SIMBA: Scalable inversion in optical tomography using deep denoising priors. IEEE Journal of Selected Topics in Signal Processing 14(6), 1163–1175 (2020)

[4] Arsenault, L.-F., Neuberg, R., Hannah, L.A., Millis, A.J.: Projected regression method for solving Fredholm integral equations arising in the analytic continuation problem of quantum physics. Inverse Problems 33(11), 115007 (2017)

[5] Cherni, A., Chouzenoux, E., Delsuc, M.-A.: PALMA, an improved algorithm for DOSY signal processing. The Analyst 142(5), 772–779 (2017)

[6] Yu, Y., Si, X., Hu, C., Zhang, J.: A review of recurrent neural networks: LSTM cells and network architectures. Neural computation 31(7), 1235–1270 (2019)

[7] Tang, W., Chouzenoux, E., Pesquet, J.-C., Krim, H.: Deep transform and metric learning network: Wedding deep dictionary learning and neural network. Neurocomputing 509, 244–256 (2022)

[8] Combettes, P.L., Wajs, V.R.: Signal recovery by proximal forward-backward splitting. Multiscale Modeling & Simulation 4(4), 1168–1200 (2005)

[9] Combettes, P.L., Pesquet, J.-C.: Proximal splitting methods in signal processing. In: al., H.H.B. (ed.) Fixed-point Algorithms for Inverse Problems in Science and Engineering. Springer Optimization and Its Applications, vol. 49, pp. 185–212. Springer, New York, NY (2011)

[10] Hegland, M.: Variable Hilbert scales and their interpolation inequalities with applications to Tikhonov regularization. Applicable Analysis 59(1-4), 207–223 (1995)

[11] Hofmann, B., Yamamoto, M.: Convergence rates for Tikhonov regularization based on range inclusions. Inverse Problems 21(3), 805 (2005)

[12] Åkesson, E.O., Daun, K.J.: Parameter selection methods for axisymmetric flame tomography through Tikhonov regularization. Applied optics 47(3), 407–416 (2008)

[13] Daun, K.J., Thomson, K.A., Liu, F., Smallwood, G.J.: Deconvolution of axisymmetric flame properties using Tikhonov regularization. Applied optics 45(19), 4638–4646 (2006)

[14] Engl, H.W.: Regularization of Inverse Problems. Kluwer Academic Publishers,
[15] Hansen, P.C.: The l-curve and its use in the numerical treatment of inverse problems (1999)

[16] Antholzer, S., Haltmeier, M., Schwab, J.: Deep learning for photoacoustic tomography from sparse data. Inverse Problems in Science and Engineering **27**(7), 987–1005 (2019)

[17] Kofler, A., Haltmeier, M., Kolbitsch, C., Kachelrieß, M., Dewey, M.: A u-nets cascade for sparse view computed tomography. In: International Workshop on Machine Learning for Medical Image Reconstruction, pp. 91–99 (2018). Springer

[18] Hammernik, K., Klatzer, T., Kobler, E., Recht, M.P., Sodickson, D.K., Pock, T., Knoll, F.: Learning a variational network for reconstruction of accelerated MRI data. Magnetic Resonance in Medicine **79**(6), 3055–3071 (2018)

[19] Adler, J., Öktem, O.: Solving ill-posed inverse problems using iterative deep neural networks. Inverse Problems **33**(12), 124007 (2017)

[20] Aggarwal, H.K., Mani, M.P., Jacob, M.: MoDL: Model-based deep learning architecture for inverse problems. IEEE Transactions on Medical Imaging **38**(2), 394–405 (2018)

[21] Meinhardt, T., Moller, M., Hazirbas, C., Cremers, D.: Learning proximal operators: Using denoising networks for regularizing inverse imaging problems. In: Proceedings of the IEEE International Conference on Computer Vision (ICCV 2017), Venice, Italy, pp. 1781–1790 (2017)

[22] Pesquet, J.-C., Repetti, A., Terris, M., Wiaux, Y.: Learning maximally monotone operators for image recovery. SIAM Journal on Imaging Sciences **14**(3) (2021)

[23] Hasannasab, M., Hertrich, J., Neumayer, S., Plonka, G., Setzer, S., Steidl, G.: Parseval proximal neural networks. Journal of Fourier Analysis and Applications **26**(4) (2020)

[24] Galinier, M., Prato, M., Chouzenoux, E., Pesquet, J.-C.: A hybrid interior point – deep learning approach for Poisson image deblurring. In: Proceedings of the IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP 2020), Espoo, Finland, pp. 1–6 (2020)

[25] Li, H., Schwab, J., Antholzer, S., Haltmeier, M.: NETT: Solving inverse problems with deep neural networks. Inverse Problems **36**(6), 065005 (2020)

[26] Bubba, T.A., Galinier, M., Lassas, M., Prato, M., Ratti, L., Siltanen, S.: Deep neural networks for inverse problems with pseudodifferential operators: an application to limited-angle tomography. SIAM Journal on Imaging Sciences **14**(2)
(2021)

[27] Gilton, D., Ongie, G., Willett, R.: Neumann networks for linear inverse problems in imaging. IEEE Transactions on Computational Imaging **6**, 328–343 (2019)

[28] Ryu, E., Liu, J., Wang, S., Chen, X., Wang, Z., Yin, W.: Plug-and-play methods provably converge with properly trained denoisers. In: Proceedings of the International Conference on Machine Learning (ICML 2019), Long Beach, CA, pp. 5546–5557 (2019)

[29] Sun, Y., Wohlberg, B., Kamilov, U.S.: An online plug-and-play algorithm for regularized image reconstruction. IEEE Transactions on Computational Imaging **5**(3), 395–408 (2019)

[30] Hurault, S., Leclaire, A., Papadakis, N.: Proximal denoiser for convergent plug-and-play optimization with nonconvex regularization. In: Proceedings of the International Conference on Machine Learning (ICML 2022), vol. 16. Baltimore, MD, pp. 9483–9505 (2022)

[31] Gregor, K., LeCun, Y.: Fast approximations of sparse coding. In: Proceedings of the International Conference on Machine Learning (ICML 2010), Haifa, Israel, pp. 399–406 (2010)

[32] Borgerding, M., Schniter, P.: Onsager-corrected deep learning for sparse linear inverse problems. In: Proceedings of the IEEE Global Conference on Signal and Information Processing (GlobalSIP 2016), Greater Washington, DC, pp. 227–231 (2016)

[33] Jin, K.H., McCann, M.T., Froustey, E., Unser, M.: Deep convolutional neural network for inverse problems in imaging. IEEE Transactions on Image Processing **26**(9), 4509–4522 (2017)

[34] Bertocchi, C., Chouzenoux, E., Corbineau, M.-C., Pesquet, J.-C., Prato, M.: Deep unfolding of a proximal interior point method for image restoration. Inverse Problems **36**(3), 034005 (2020)

[35] Natarajan, S.K., Caro, M.A.: Particle swarm based hyper-parameter optimization for machine learned interatomic potentials. arXiv preprint arXiv:2101.00049 (2020)

[36] Monga, V., Li, Y., Eldar, Y.C.: Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing. IEEE Signal Processing Magazine **38**(2), 18–44 (2021)

[37] Gharbi, M., Chouzenoux, E., Pesquet, J.-C.: An unrolled half-quadratic approach for sparse signal recovery in spectroscopy. Signal Processing **218**, 109369 (2024)
[38] Le, H.T.V., Repetti, A., Pustelnik, N.: Unfolded proximal neural networks for robust image Gaussian denoising. IEEE Transactions on Image Processing **33**, 4475–4487 (2024)

[39] Genzel, M., Macdonald, J., März, M.: Solving inverse problems with deep neural networks–robustness included? IEEE Transactions on Pattern Analysis and Machine Intelligence **45**(1), 1119–1134 (2023)

[40] Combettes, P.L., Pesquet, J.-C.: Deep neural network structures solving variational inequalities. Set-Valued and Variational Analysis, 1–28 (2020)

[41] Combettes, P.L., Pesquet, J.-C.: Lipschitz certificates for layered network structures driven by averaged activation operators. SIAM Journal on Mathematics of Data Science **2**(2), 529–557 (2020)

[42] Serrurier, M., Mamalet, F., González-Sanz, A., Boissin, T., Loubes, J.-M., Barrio, E.: Achieving robustness in classification using optimal transport with hinge regularization. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2021), virtual (2021)

[43] Cisse, M., Bojanowski, P., Grave, E., Dauphin, Y., Usunier, N.: Parseval networks: Improving robustness to adversarial examples. In: Proceedings of the International Conference on Machine Learning (ICML 2017), Sydney, Australia, pp. 854–863 (2017)

[44] Neacșu, A., Pesquet, J.-C., Burileanu, C.: EMG-based automatic gesture recognition using Lipschitz-regularized neural networks. ACM Transactions on Intelligent Systems and Technology **18**, 1–25 (2024)

[45] Neacșu, A., Pesquet, J.-C., Vasilescu, V., Burileanu, C.: ABBA neural networks: Coping with positivity, expressivity, and robustness. SIAM Journal on Mathematics of Data Science **6**, 649–678 (2024)

[46] Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces vol. 408. Springer, New York (2017)

[47] Terris, M., Abdulaziz, A., Dabbech, A., Jiang, M., Repetti, A., Pesquet, J.-C., Wiaux, Y.: Deep post processing for sparse image deconvolution. In: Signal Processing with Adaptive Sparse Structured Representations (SPARS), Toulouse, France, p. 2 (2019)