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

HAL
open science

# Computable Lipschitz Bounds for Deep Neural Networks

Moreno Pintore, Bruno Després

► To cite this version:

Moreno Pintore, Bruno Després. Computable Lipschitz Bounds for Deep Neural Networks. 2024. [https://arxiv.org/abs/2405.04756](https://arxiv.org/abs/2405.04756)

HAL Id: hal-04756410

https://inria.hal.science/hal-04756410v1

Preprint submitted on 29 Oct 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COMPUTABLE LIPSCHITZ BOUNDS FOR DEEP NEURAL NETWORKS

MORENO PINTORE\* AND BRUNO DESPRÉS†

**Abstract.** Deriving sharp and computable upper bounds of the Lipschitz constant of deep neural networks is crucial to formally guarantee the robustness of neural-network based models. We analyse three existing upper bounds written for the $l^2$ norm. We highlight the importance of working with the $l^1$ and $l^\infty$ norms and we propose two novel bounds for both feed-forward fully-connected neural networks and convolutional neural networks. We treat the technical difficulties related to convolutional neural networks with two different methods, called explicit and implicit. Several numerical tests empirically confirm the theoretical results, help to quantify the relationship between the presented bounds and establish the better accuracy of the new bounds. Four numerical tests are studied: two where the output is derived from an analytical closed form are proposed; another one with random matrices; and the last one for convolutional neural networks trained on the MNIST dataset. We observe that one of our bound is optimal in the sense that it is exact for the first test with the simplest analytical form and it is better than other bounds for the other tests.

**Key words.** Lipschitz constant, deep neural networks, upper bounds, robustness.

**MSC codes.** 26A16, 68Q17, 68T99, 92B20.

**1. Introduction.** Deep neural networks are more and more common in most scientific applications, even though they are still unstable in presence of specific small input perturbations. The most common examples of such instabilities are the adversarial attacks in the context of image classification [28]. We also remark that, for engineering applications, the stability of the deep neural networks approximating physical functions with respect to small perturbations has to be controlled and quantified. A recent study in this direction is [4]. That is why it is of supreme importance for the development of the discipline to establish methods which qualify and quantify the stability of functions or operators represented by deep neural networks. The interest in the stability of the neural networks is also partially inspired by some of the previous works of the manuscript authours, where the theory of deep neural networks is analyzed from different perspectives (see [2, 3, 11]).

In [31], Szegedy-Goodfellow-et-al analyse the relationship between the Lipschitz constant of deep neural networks and their stability properties [15], proposing an upper bound obtained by multiplying the Lipschitz constants of all the layers. This bound, denoted by $K_*$ in the following discussion, is very pessimistic and so it is meaningless when the number of layers grows (the deep regime). Nevertheless, $K_*$ has been extensively utilized during the training or in the architecture of deep neural networks to improve the stability of the model [9, 27, 17]. Sharper estimates have been obtained in Combettes-Pesquet [8] and in Scaman-Virmaux [32]. In the current manuscript, we are interested in further extending these two specific works by presenting new certified Lipschitz bounds in $l^p$ norms ($1 \le p \le \infty$) and their extensions to convolutional neural networks. We are also interested in measuring and comparing the efficiency of all these estimates on simple tests.

Even though it is not the subject of the current paper, we highlight that there exist some estimates of the Lipschitz constant based on semidefinite minimization problems [14] or on polynomial approximation approaches [23, 7] (see [18] for a compact overview of both approaches), as well as estimates of the local Lipschitz constant [33, 35, 25]. We also highlight that all these bounds can be employed in theoretical analyses such as [21, 20] to better characterize the neural network output.

In our opinion, our main results are as follows.

* We provide theoretical reasons why the $l^1$ norm and the $l^\infty$ norm are in some cases preferable in

<footnote>\*MEGAVOLT team, INRIA, F-75013 Paris, France – SCAI, Sorbonne Université, F-75005 Paris, France. (moreno.pintore@inria.fr).</footnote>

<footnote>†Sorbonne Université, Université Paris Cité, CNRS, INRIA, Laboratoire Jacques-Louis Lions, LJLL, F-75005 Paris, France, (bruno.despres@sorbonne-universite.fr).</footnote>
evaluating the Lipschitz bounds, rather than $l^2$-based norms which are more standard.

* We define new bounds $K_3$ and $K_4$. The bound $K_4$ is sharper than the $K_1$ constant (which comes from the Combettes-Pesquet work [8]) in the $l^1$ and $l^\infty$ norms.

* We extend in two different ways the discussed bounds to the case of convolutional neural networks.

* We test our bounds on four numerical tests of different origins which confirm the better accuracy provided by our new bound $K_4$. In particular, we observe a surprising behavior for our second test problem for which we know the exact value of the Lipschitz constant. In this case, the numerical bound $K_4$ is exact as shown in Table 2.

The paper is structured as follows. In Section 2, we focus on fully-connected feed-forward neural networks. In particular, for this kind of neural networks, we summarize and prove existing Lipschitz upper bounds in Section 2.1 and we prove new bounds in Section 2.2. Convolutional neural networks are analysed in Section 3, where two approaches to bound the Lipschitz constants of networks including max-pooling layers are discussed in Sections 3.1 and 3.3. Numerical results in strong agreement with the theoretical ones are provided in Section 4, where we consider a feed-forward fully-connected neural network with random weights [8], with weights chosen to efficiently approximate polynomial functions and convolutional neural networks. In Sections 4.2 and 4.3 we propose two benchmarks where the output of a deep neural network is known in a simple closed form, because we believe that this is crucial in order to improve the theoretical understanding of deep neural networks with an arbitrary number of layers but in simplified scenarios. Section 5 includes conclusion remarks and future perspectives.

**Notation** For any vector $x \in \mathbb{R}^a$ with arbitrary size $a \ge 1$, we will denote by $\|x\|_{l^p(\mathbb{R}^a)} = (\sum_{i=1}^a |x_i|^p)^{1/p}$ its $l^p$ norm. Popular norms [19] are the $l^1$ norm $\|x\|_{l^1(\mathbb{R}^a)} = \sum_{i=1}^a |x_i|$, the $l^2$ norm $\|x\|_{l^2(\mathbb{R}^a)} = \sum_{i=1}^a x_i^2$ and the $l^\infty$ norm $\|x\|_{l^\infty(\mathbb{R}^a)} = \max_{i=1}^a |x_i|$. It yields the induced norm for matrices $M \in \mathcal{M}_{ab}(\mathbb{R}) = \mathbb{R}^{a \times b}$ where the output size $a \ge 1$ and the input size $b \ge 1$ are arbitrary. The induced norm is

$$ (1.1) \qquad \|M\|_{l^p(\mathcal{M}_{ab}(\mathbb{R}))} = \max_{x \neq 0} \frac{\|Mx\|_{l^p(\mathbb{R}^a)}}{\|x\|_{l^p(\mathbb{R}^b)}}. $$

It will be denoted by the simpler form $\|M\| = \|M\|_{l^p(\mathcal{M}_{ab}(\mathbb{R}))}$ if there is no ambiguity.

**2. Fully-connected feed-forward neural networks.** A generic fully-connected feed-forward neural network function $f$ [16, 11] is written under the form

$$ (2.1) \qquad f = f_{\ell} \circ S_{\ell} \circ f_{\ell-1} \circ S_{\ell-1} \circ \dots \circ S_1 \circ f_0. $$

Feed-forward neural networks are mostly used for regression purposes, which means that $f$ is usually close to a smooth objective function [11]. Here the functions $f_i$, $i = 1, \dots, \ell$, denote affine functions with varying input and output dimensions ($a_0, a_1, \dots, a_{\ell+1}) \in \mathbb{N}^{\ell+2}$ with $a_i > 0$, for $i = 1, \dots, \ell$. More precisely, $f_i(x_i) = W_i x_i + b_i \in \mathbb{R}^{a_{i+1}}$ for all $x_i \in \mathbb{R}^{a_i}$. The parameters of the affine function $f_i$ are called the weights $W_i \in \mathcal{M}_{a_{i+1}, a_i}(\mathbb{R})$ and the biases $b_i \in \mathbb{R}^{a_{i+1}}$. The intermediate functions $S_i$, $i = 1, \dots, \ell$, are nonlinear activation functions [30] such that $0 \le S'_i(x) \le 1$ a.e. $x \in \mathbb{R}$. The ReLU function $R(x) = \max(0, x)$ is a popular activation function that we will use in our numerical tests. The activation functions $S_i$ are all Lipschitz by assumption. Note that in practice, that is in the softwares that are used in machine learning, the derivative of all activation functions are defined everywhere. For example a standard choice is $R'(0) = 0$ (see [1, 29, 5]). Any activation function is generalized component wise to vectors of arbitrary dimension. The index $i$ in $S_i$ explains that one can change the activation functions from one layer to the other. The number $a_i$ is referred to as the number of neurons of the $i$-th layer. The integer $\ell$ is the number of layers of the neural network function (2.1).
The regularity of $f$ is

$$f \in C^0(\mathbb{R}^{a_0})^{a_{\ell+1}} \cap \text{Lip}(\mathbb{R}^{a_0})^{a_{\ell+1}}$$

since it is the composition of continuous and globally Lipschitz functions. The Rademacher's Theorem states that $f$ is differentiable almost everywhere. The Lipschitz constant of the function $f$ that we consider is

$$L = \sup_{x \in \mathbb{R}^{a_0}} \|\nabla f(x)\|_{lp}(M_{a_{\ell+1}, a_0}(\mathbb{R})) = \left( \|\nabla f\|_{lp}(M_{a_{\ell+1}, a_0}(\mathbb{R})) \right)_{L^{\infty}(\mathbb{R}^{a_0})}.$$

Let us denote by $D_r = D_r(x)$ the square diagonal matrix $D_r = S_r' \circ f_{r-1} \circ S_{r-1} \circ \dots \circ S_1 \circ f_0(x)$. With our assumptions, the diagonal coefficients of $D_r$ are between 0 and 1. We define the set

$$ (2.2) \qquad \mathcal{D}_r = \{D_r \in \mathcal{M}_{a_r, a_r}(\mathbb{R}) \mid D_r \text{ is a diagonal matrix with coefficients between 0 and 1}\}. $$

One can thus write $D_r \in \mathcal{D}_r$. The chain rule shows that the gradient of $f$ is the product of matrices

$$ (2.3) \qquad \nabla f(x) = W_{\ell} D_{\ell}(x) W_{\ell-1} D_{\ell-1}(x) \dots D_1(x) W_0, $$

where only the matrices $D_r$, $r = 1, \dots, \ell$, depend on $x$.

For convenience, we note $\mathcal{D} = \mathcal{D}_{\ell} \times \mathcal{D}_{\ell-1} \times \dots \times \mathcal{D}_1$, $D = D(x) = (D_{\ell}(x), D_{\ell-1}(x), \dots, D_1(x)) \in \mathcal{D}$ and $W = (W_{\ell}, W_{\ell-1}, \dots, W_0)$. As in [31], we introduce the quantity

$$ (2.4) \qquad K = K(W) = \max_{D \in \mathcal{D}} \|W_{\ell} D_{\ell} W_{\ell-1} D_{\ell-1} \dots D_1 W_0\|, $$

so that one clearly has

$$ L \le K. $$

**Remark 2.1 (Complexity).** Set $\text{Ext}(\mathcal{D}_r) = \{D_r \in \mathcal{D}_r \mid$ the diagonal coefficients are equal to 0 or 1}$ which is the finite ensemble of exterior points of the convex set $\mathcal{D}_r$. In (2.4) one notices that $K$ is a convex function of $D_r$, so the maximal value $K$ is reached at the external points. Therefore one also has

$$ K = \max_{D \in \text{Ext}(\mathcal{D})} \|W_{\ell} D_{\ell} W_{\ell-1} D_{\ell-1} \dots D_1 W_0\| $$

where $\text{Ext}(\mathcal{D}) = \text{Ext}(\mathcal{D}_{\ell}) \times \dots \times \text{Ext}(\mathcal{D}_0)$. Since $\text{Ext}(\mathcal{D})$ is finite, it is possible to calculate $K$ by exhaustion of the possibilities. Nevertheless the number of cases is the number of corners of the hypercube $\mathcal{D} \in \mathbb{R}^{a_{\ell}+\dots+a_0}$. So the complexity of this calculation is $O(2^{a_{\ell}+\dots+a_0})$ where $a_{\ell} + \dots + a_0$ is the total number of neurons. That is, the complexity of the calculation of $K$ is exponential with respect to the total number of neurons. In many reasonable cases, this cost makes this direct calculation of $K$ just impossible.

Our goal in this work is to examine some bounds on $K$ which are computable and sharp. This is what we call **certified bounds**.

**2.1. Previous works.** In this section we introduce some upper bounds for the Lipschitz constant of fully-connected feed-forward neural networks available in literature.

**DEFINITION 2.2 (Worst bound).** $K_{\star} := \prod_{r=0}^{\ell} \|W_r\|$.

The first bound $K \le K_{\star}$ is evident from the sub-multiplicativity of the norm [8, 32, 11]. However, it will be clear in the numerical Section that this bound is not sharp and is actually the worst one.

The next result is an excerpt from [8]. For convenience, we begin with a definition.
**DEFINITION 2.3.** Given $0 \le s < t \le \ell + 1$, define the matrices $W_{(t,s)} = W_{t-1}W_{t-2}\dots W_{s+1}W_s$.

**THEOREM 2.4** (Combettes-Pesquet bound). One has the bound $K \le K_1$ where

$$ (2.5) \qquad K_1 = \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \|W_{(\ell+1,r_n)}\| \|W_{(r_n,r_{n-1})}\| \dots \|W_{(r_2,r_1)}\| \|W_{(r_1,0)}\| $$

*Remark 2.5.* a) Necessarily $n \le \ell$ in the formula (2.5), and b) the number of terms in the sum is $2^\ell$.

*Remark 2.6.* We prove below in Proposition 2.7 that $K_1 \le K_\star$. But before proving in the inequality, we stress that there are two intuitive reasons to argue that $K_1$ could be much smaller than $K_\star$, that is $K_1 \ll K_\star$, in some cases. The first one is that the decomposition for $n = \ell$ gives back one term which is equal to $\frac{1}{2^\ell} K_\star$. It results in a drastic reduction of the influence of the worst bound in the deep regime. The second one is that for all $n < \ell$, there is a least one matrix $W_{(r',r)}$ with $r' \ge r + 2$. Then this term is the product of at least two matrices, for which one can expect some cancellations and thus a strong reduction in the numerical value of the norm of the product.

*Proof.* The technical part of the proposed proof is a simplification of the one in [8]. We also show that the bound (2.5) holds for all induced norms, while the seminal reference [8] considers only the $l^2$ norm.

One expresses the matrices $D_r \in \mathcal{M}_{a_r, a_r}(\mathbb{R})$ in (2.4) as

$$ (2.6) \qquad D_r = \frac{1}{2}(I_r + Z_r) $$

where $I_r \in \mathcal{M}_{a_r, a_r}(\mathbb{R})$ is the identity matrix and $Z_r = 2D_r - I_r$ is such that $\|Z_r\| \le 1$. We denote by $\mathcal{Z}_r$ the set of matrices $\mathcal{Z}_r = \{Z_r | \exists D_r \in \mathcal{D}_r : Z_r = 2D_r - I_r\}$. This is generalized as $\mathcal{Z} = \mathcal{Z}_\ell \times \dots \times \mathcal{Z}_1$ and we introduce the notation

$$ Z = (Z_\ell, \dots, Z_1) \in \mathcal{Z}. $$

To explain the interest of the matrices $Z_r$, let us consider the following example with 4 hidden layers

$$ (2.7) \qquad f = f_4 \circ S_4 \circ f_3 \circ S_3 \circ f_2 \circ S_2 \circ f_1 \circ S_1 \circ f_0. $$

The gradient formula (2.3) becomes

$$ (2.8) \qquad \nabla f = \frac{1}{2^4} W_4(I + Z_4)W_3(I + Z_3)W_2(I + Z_2)W_1(I + Z_1)W_0 $$

where the identity matrices have different sizes in the general case. Expansion of all terms yields

$$ (2.9) \qquad \nabla f = \frac{1}{2^4} (W_4 W_3 W_2 W_1 W_0 + \dots + W_4 Z_4 W_3 Z_3 W_2 Z_2 W_1 Z_1 W_0) $$

where in the parenthesis, the first term is the product of all weight matrices, the last term is the product of all weight matrices combined with the $Z_r$, and the middle term ... represents all other possible chunks of weight matrices $W_r$ separated with matrices $Z_r$. Then a triangular inequality yields the result because $\|Z_r\| = 1$ for all $r$.
The general case is treated as follows.

$$
\begin{align*}
K &= \max_{D \in \mathcal{D}} \|W_{\ell} D_{\ell} W_{\ell-1} D_{\ell-1} \dots D_1 W_0\| \\
&= \max_{Z \in \mathcal{Z}} \left\| W_{\ell} \left( \frac{1}{2}(I_{\ell} + Z_{\ell}) \right) W_{\ell-1} \left( \frac{1}{2}(I_{\ell-1} + Z_{\ell-1}) \right) \dots \left( \frac{1}{2}(I_1 + Z_1) \right) W_0 \right\| \\
&\le \frac{1}{2^{\ell}} \max_{Z \in \mathcal{Z}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)} Z_{r_n} W_{(r_n,r_{n-1})} \dots W_{(r_2,r_1)} Z_{r_1} W_{(r_1,0)} \| \\
&\le \frac{1}{2^{\ell}} \max_{Z \in \mathcal{Z}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)} \| \| Z_{r_n} \| \| W_{(r_n,r_{n-1})} \| \dots \| W_{(r_2,r_1)} \| \| Z_{r_1} \| \| W_{(r_1,0)} \| \\
&\le \frac{1}{2^{\ell}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)} \| \| W_{(r_n,r_{n-1})} \| \dots \| W_{(r_2,r_1)} \| \| W_{(r_1,0)} \| = K_1,
\end{align*}
$$

where we have removed the dependence on Z by using the bound $||Z_r|| \le 1$. □

PROPOSITION 2.7. One has the bound $K_1 \le K_*$.

Proof. One has

$$
(2.11) \qquad \|W_{(t,s)}\| = \|W_{t-1}W_{t-2}\dots W_{s+1}W_s\| \le \|W_{t-1}\|\|W_{t-2}\|\dots \|W_{s+1}\|\|W_s\|
$$

because of the norm submultiplicativity. We now bound all the norms involved in the definition of $K_1$ as in
2.11 obtaining

$$
\begin{align*}
K_1 &= \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)} \| \| W_{(r_n,r_{n-1})} \| \dots \| W_{(r_2,r_1)} \| \| W_{(r_1,0)} \| \\
&\le \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_\ell \| \| W_{\ell-1} \| \dots \| W_1 \| \| W_0 \| = \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} K_* = K_*.
\end{align*}
$$

The last equality holds because the sum includes $2^\ell$ identical terms. □

The next result is an extension to generic $l^p$ matrix norms of the Virmaux-Scaman upper bound [32].
One starts from the generic decomposition of $W_i \in \mathcal{M}_{a_{i+1},a_i}(\mathbb{R})$ written as

$$
W_i = U_i \Sigma_i V_i^T.
$$

This decomposition is a singular value decomposition [6] by choosing $U_i \in \mathcal{M}_{a_{i+1},a_{i+1}}(\mathbb{R})$ and $V_i \in \mathcal{M}_{a_i,a_i}(\mathbb{R})$ as orthogonal matrices, and by imposing that $\Sigma$ is a diagonal matrix with positive entries. Other choices are of course possible.

THEOREM 2.8 (Virmaux-Scaman bound). One has the bound $K \le K_2$ where

$$
(2.12) \qquad K_2 = \prod_{i=0}^{\ell-1} \max_{D_{i+1} \in \mathcal{D}_{i+1}} \| \tilde{\Sigma}_{i+1} V_{i+1}^T D_{i+1} U_i \tilde{\Sigma}_i \|,
$$

where $\tilde{\Sigma}_0 = \Sigma_0 V_0^T$, $\tilde{\Sigma}_\ell = U_\ell \Sigma_\ell$ and $\tilde{\Sigma}_i = \Sigma_i^{1/2}$ for all $i \notin \{0, \ell\}$.
Proof. The theorem can be proved in the following way:

$$
\begin{align*}
K &= \max_{D \in \mathcal{D}} \| (U_\ell \Sigma_\ell V_\ell^T) D_\ell (U_{\ell-1} \Sigma_{\ell-1} V_{\ell-1}^T) D_{\ell-1} \dots D_2 (U_1 \Sigma_1 V_1^T) D_1 (U_0 \Sigma_0 V_0^T) \| \\
&= \max_{D \in \mathcal{D}} \| \tilde{\Sigma}_\ell V_\ell^T D_\ell U_{\ell-1} \tilde{\Sigma}_{\ell-1} \tilde{\Sigma}_{\ell-1} V_{\ell-1}^T D_{\ell-1} \dots D_2 U_1 \tilde{\Sigma}_1 \tilde{\Sigma}_1 V_1^T D_1 U_0 \tilde{\Sigma}_0 \| \\
&\le \prod_{i=0}^{\ell-1} \max_{D_{i+1} \in \mathcal{D}_{i+1}} \| \tilde{\Sigma}_{i+1} V_{i+1}^T D_{i+1} U_i \tilde{\Sigma}_i \|.
\end{align*}
$$

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

The prescription from [32] is to take $\tilde{\Sigma}_0 = \Sigma_0$ and $\tilde{\Sigma}_\ell = \Sigma_\ell$, when the $l^2$ norm is used, because $V_0^T$ and $U_\ell$ are orthogonal matrices. Various similar bounds can be obtained with other matrix factorizations.

**Remark 2.9** (Complexity of the Virmaux-Scaman bound). By comparison with Remark 2.1, one can calculate the Virmaux-Scaman bound by varying all $D_{i+1} \in \text{Ext}(\mathcal{D}_{i+1})$ independently one of the other. Therefore a direct calculation is possible with a cost $O(2^{a_\ell} + \dots + 2^{a_0})$. Clearly, this is much less than the complexity of $K$. However, this is still impracticable if the number of neurons on one layer is reasonably large (which is the case for standard neural networks). For this reason, in [32] the authors propose to approximate $K_2$ instead of computing it exactly.

**2.2. New certified bounds.** In this section, we introduce two new upper bounds for the Lipschitz constant of fully-connected feed-forward neural networks. Since we will need either one of the equality expressed in (2.13), the symbol $\| \cdot \|$ will denote in this Section either $\| \cdot \|_1$ or $\| \cdot \|_\infty$ (and only one of these two norms).

**Definition 2.10** (Element-wise absolute value of a matrix). Given a matrix $A \in \mathcal{M}_{m,n}(\mathbb{R})$ with entries $\{A_{ij}\}_{ij}$, we denote by $A^{\text{abs}} \in \mathcal{M}_{m,n}(\mathbb{R})$ the matrix obtained by applying the absolute value to each entry of $A$, i.e. $(A^{\text{abs}})_{ij} = |A_{ij}|, \forall i, j$.

From the well known [19] identities $\|A\|_1 = \max_{1 \le i \le n} \sum_{j=1}^m |A_{ij}|$ and $\|A\|_\infty = \max_{1 \le i \le m} \sum_{j=1}^n |A_{ij}|$, one obtains the equalities

$$ (2.13) \qquad \|A\|_1 = \|A^{\text{abs}}\|_1 \text{ and } \|A\|_{\infty} = \|A^{\text{abs}}\|_{\infty}. $$

For other norms, then the norm of a matrix may be different from the norm of its absolute value. Consider for example $A = \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}$ together with the euclidian norm $\|A\|_2$. Then it is easy to check that $\|A\|_2 = \sqrt{2} < \|A^{\text{abs}}\|_2 = 2$. We introduce two evident lemmas used later in the rest of the Section.

**LEMMA 2.11.** Let $A_1, A_2, \dots, A_N$ be $N$ matrices which are compatible in the sense that the product $A_N A_{N-1} \dots A_2 A_1 \in \mathcal{M}_{m,n}(\mathbb{R})$ is well defined. Then one has

$$ |(A_N A_{N-1} \dots A_2 A_1)_{i,j}| \le (A_N^{\text{abs}} A_{N-1}^{\text{abs}} \dots A_2^{\text{abs}} A_1^{\text{abs}})_{i,j} \quad \text{for all } 1 \le i \le m \text{ and } 1 \le j \le n $$

$$ \text{and } \|A_N A_{N-1} \dots A_2 A_1\| \le \|A_N^{\text{abs}} A_{N-1}^{\text{abs}} \dots A_2^{\text{abs}} A_1^{\text{abs}}\|. $$

Proof. The lemma is proved by expressing the entries of the product of matrices in terms of the entries
$a_{i,j}^r$ of the initial matrices $A_r$, $r = 1, \dots, N$. One checks that

$$
\begin{align*}
(A_N A_{N-1} \dots A_2 A_1)_{i,j}^{\text{abs}} &= \left| \sum_{k_1, k_2, \dots, k_{N-1}} a_{i,k_{N-1}}^N a_{k_{N-1},k_{N-2}}^{N-1} \dots a_{k_2,k_1}^2 a_{k_1,j}^1 \right| \\
&\le \sum_{k_1, k_2, \dots, k_{N-1}} \left| a_{i,k_{N-1}}^N a_{k_{N-1},k_{N-2}}^{N-1} \dots a_{k_2,k_1}^2 a_{k_1,j}^1 \right| \\
&= \sum_{k_1, k_2, \dots, k_{N-1}} \left| a_{i,k_{N-1}}^N \right| \left| a_{k_{N-1},k_{N-2}}^{N-1} \right| \dots \left| a_{k_2,k_1}^2 \right| \left| a_{k_1,j}^1 \right| \\
&= (A_N^{\text{abs}} A_{N-1}^{\text{abs}} \dots A_2^{\text{abs}} A_1^{\text{abs}})_{i,j}.
\end{align*}
$$

The definitions of the $l^1$ norm and of the $l^\infty$ norm yield the second inequality. $\square$

LEMMA 2.12. Take matrices $A, C \in \mathcal{M}_{m,n}(\mathbb{R})$ and $B, D \in \mathcal{M}_{n,p}(\mathbb{R})$ be such that $0 \le A_{i,j} \le C_{i,j}$ and $0 \le B_{i,j} \le D_{i,j}$ for all values of $i$ and $j$. Then one has $(AB)_{i,j} \le (CD)_{i,j}, \forall i = 1, \dots, m, j = 1, \dots, p, \text{ and} \\ \|AB\| \le \|CD\|$. $\square$

Proof. Indeed one has $(AB)_{i,j} = \sum_{k=1}^n A_{i,k} B_{k,j} \le \sum_{k=1}^n C_{i,k} D_{k,j} = (CD)_{i,j}$ for all $i = 1, \dots, m$ and $j = 1, \dots, p$. The definitions of the $l^1$ and $l^\infty$ norms yield the second inequality. $\square$

THEOREM 2.13. One has the bound $K \le K_3$ where $K_3 = \|W_\ell^{\text{abs}} W_{\ell-1}^{\text{abs}} \dots W_0^{\text{abs}}\|$. Moreover, this new bound is sharper than the worst bound, that is $K_3 \le K_\star$.

Proof. The proof is made for the $l^\infty$ norm. One has

$$
\begin{align*}
K &= \max_{D \in \mathcal{D}} \|W_\ell D_\ell W_{\ell-1} D_{\ell-1} \dots D_1 W_0\|_{\infty} \\
&= \max_{D \in \mathcal{D}} \max_i \sum_j \left| \sum_{k_\ell, \dots, k_1} w_{i,k_\ell}^\ell d_{k_\ell}^\ell w_{k_\ell, k_{\ell-1}}^{\ell-1} \dots d_{k_1}^1 w_{k_1,j}^0 \right| \\
&\le \max_{D \in \mathcal{D}} \max_i \sum_j \sum_{k_\ell, \dots, k_1} \left| w_{i,k_\ell}^\ell d_{k_\ell}^\ell w_{k_\ell, k_{\ell-1}}^{\ell-1} \dots d_{k_1}^1 w_{k_1,j}^0 \right| \\
&\le \max_i \sum_j \sum_{k_\ell, \dots, k_1} \left| w_{i,k_\ell}^\ell \right| \left| w_{k_\ell, k_{\ell-1}}^{\ell-1} \right| \dots \left| w_{k_1,j}^0 \right| = \|W_\ell^{\text{abs}} W_{\ell-1}^{\text{abs}} \dots W_0^{\text{abs}}\|_{\infty}.
\end{align*}
$$

The proof for the $l^1$ norm is derived by exchanging $i$ and $j$ in (2.14). It yields the first part of the claim.

The second part is obtained as follows. One has $K_3 = \|W_\ell^{\text{abs}} W_{\ell-1}^{\text{abs}} \dots W_0^{\text{abs}}\| \le \|W_\ell^{\text{abs}}\| \|W_{\ell-1}^{\text{abs}}\| \dots \|W_0^{\text{abs}}\|$ by sub-multiplicativity of the norm. Using (2.13) one has $\|W_\ell^{\text{abs}}\| \|W_{\ell-1}^{\text{abs}}\| \dots \|W_0^{\text{abs}}\| = \|W_\ell\| \|W_{\ell-1}\| \dots \|W_0\|$, so $K_3 \le K_\star$. $\square$

Next, we combine the Combettes-Pesquet technique with the bound $K_3$ and our choice of norm (either $\| \cdot \|_1$ or $\| \cdot \|_\infty$). It defines a new bound

$$
(2.15) \qquad K_4 = \frac{1}{2\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \|.
$$

THEOREM 2.14. One has the bounds

$$
(2.16) \qquad K \le K_4 \le K_1.
$$
Proof. We start by proving $K \le K_4$. The third line in (2.10) rewrites as

$$ (2.17) \qquad K \le \frac{1}{2^\ell} \max_{Z \in \mathcal{Z}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)} Z_{r_n} W_{(r_n,r_{n-1})} \dots W_{(r_2,r_1)} Z_{r_1} W_{(r_1,0)} \|. $$

Lemma 2.11 shows that

$$ (2.18) \qquad \| W_{(\ell+1,r_n)} Z_{r_n} W_{(r_n,r_{n-1})} \dots W_{(r_2,r_1)} Z_{r_1} W_{(r_1,0)} \| \le \| W_{(\ell+1,r_n)}^{\text{abs}} Z_{r_n}^{\text{abs}} W_{(r_n,r_{n-1})}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} Z_{r_1}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \|. $$

Generically, one has $(Z_r^{\text{abs}})_{ii} \le 1$ for all $i$ and $(Z_r^{\text{abs}})_{ij} = 0$ for all $i \ne j$. That is the matrix $Z_r$ is element-wise dominated by the identity matrix of the same size. Lemma 2.12 used recursively yields that

$$ (2.19) \qquad \| W_{(\ell+1,r_n)}^{\text{abs}} Z_{r_n}^{\text{abs}} W_{(r_n,r_{n-1})}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} Z_{r_1}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \| \le \| W_{(\ell+1,r_n)}^{\text{abs}} W_{(r_n,r_{n-1})}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \|. $$

The combination of (2.17–2.19) yields the claim.

The inequality $K_4 \le K_1$ comes from

$$ \begin{aligned} K_4 &\le \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)}^{\text{abs}} \| \dots \| W_{(r_2,r_1)}^{\text{abs}} \| \| W_{(r_1,0)}^{\text{abs}} \| \\ &\le \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)}^{\text{abs}} \| \dots \| W_{(r_2,r_1)}^{\text{abs}} \| \| W_{(r_1,0)}^{\text{abs}} \| = K_1 \end{aligned} $$

thanks to the norm submultiplicativity and (2.13). □

PROPOSITION 2.15. One has the bound $K_4 \le K_3$.

Proof. Consider the generic term in the sum of $K_4$

$$ (2.20) \qquad \| W_{(\ell+1,r_n)}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \| $$

for suitable indices $1 \le r_1 < \dots < r_n \le \ell$. Then consider the generic matrix $W_{r_{i+1},r_i}^{\text{abs}}$ in (2.20), $i = 0, \dots, n$, where $r_0 = 0$ and $r_{n+1} = \ell + 1$. For all admissible indices $i, j$, the following inequality holds as a consequence of Lemma 2.11

$$ (2.21) \qquad \left( W_{r_{k+1},r_k}^{\text{abs}} \right)_{i,j} = \left| \left( W_{r_{k+1},r_k} \right)_{i,j} \right| \le \left( W_{r_{k+1}-1}^{\text{abs}} W_{r_{k+1}-2}^{\text{abs}} \dots W_{r_{k+1}}^{\text{abs}} W_{r_k}^{\text{abs}} \right)_{i,j}. $$

We can recursively use Lemma 2.12 to claim that

$$ (2.22) \qquad \| W_{(\ell+1,r_n)}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \| \le \| W_{\ell}^{\text{abs}} \dots W_1^{\text{abs}} W_0^{\text{abs}} \|. $$

Therefore

$$ \begin{aligned} K_4 &= \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{(\ell+1,r_n)}^{\text{abs}} \dots W_{(r_2,r_1)}^{\text{abs}} W_{(r_1,0)}^{\text{abs}} \| \\ &\le \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \| W_{\ell}^{\text{abs}} \dots W_1^{\text{abs}} W_0^{\text{abs}} \| = \frac{1}{2^\ell} 2^\ell \| W_{\ell}^{\text{abs}} \dots W_1^{\text{abs}} W_0^{\text{abs}} \| = K_3, \end{aligned} $$

where we use the definition of $K_3$ and the fact that the sum comprises $2^\ell$ terms that can be bounded by the same quantity independent from the subscripts $r_1, \dots, r_n$. □
In summary, we derived the series of inequalities

$$L \le K \le K_4 \le \min(K_1, K_3) \le \max(K_1, K_3) \le K_*.$$

With this respect $K_4$ is a better upper bound than the others. How much better? It will be evaluated in the numerical Section with basic numerical experiments.

**3. Convolutional neural networks.** A CNN is a specialized kind of neural network [16, 11] for processing data that have a known, grid-like topology. CNNs are tremendously successful in practical applications. It is therefore appealing to examine under what conditions the previous material, developed for fully-connected feed-forward neural networks, can be generalized to CNNs. A generic CNN function [26] is written under the form

$$(3.1) \qquad g = g_{\ell} \circ T_{\ell} \circ g_{\ell-1} \circ T_{\ell-1} \circ \cdots \circ T_1 \circ g_0.$$

This structure formally seems very similar to the one of feed-forward neural networks (2.1) but there are important differences which are summarized in the next Remarks.

*Remark 3.1.* In (2.1) the $f_r$ are linear and have ability to change the dimension (of the underlying space), while the $S_r$ are non linear activation functions and do not change the dimension. In (3.1) the notation cannot be as simple. We choose arbitrarily that the $g_r$ have ability to change the dimension and that the $T_r$ do not change the dimension. It means that some $g_r$ can be non linear while some $T_r$ can be linear.

*Remark 3.2.* The main example that motivates the classification made in above Remark 3.1 is the max-pooling function detailed below. Indeed the max-pooling function is non linear and changes the dimension so it cannot fit in the simple structure (2.1). Since a max-pooling function $g_r$ is, to the best of our knowledge, rarely followed by a non linear function, it is necessary to allow the functions $T_{r+1}$ to be equal to the identity if it is needed, that is $T_{r+1} = I$. Then it makes $T_{r+1}$ a linear function. In summary all $g_r$ and all $T_r$ can be linear or non linear in our modeling of convolutional neural networks.

*Remark 3.3.* As max-pooling layers, also average-pooling layers are rarely followed by non linear activation functions. In this case, the operator $g_{r+1} \circ T_r \circ g_r$, composed by an average-pooling layer $g_r$, the identity operator $T_r$ and another linear operator $g_{r+1}$, is linear and can be exactly represented by a single linear layer which matrix is the product between the matrices associated with $g_{r+1}$ and $g_r$. Since this merging operation is always possible and preserve the structure (3.1), in the following we always assume that all triplets of consecutive linear operators are already merged, i.e. if $g_{r+1}$ and $g_r$ are linear, then $T_r$ is non linear.

*Remark 3.4.* It is common for CNN [16] to finish with another non linear function $T_{\ell+1}$ called softmax function [11], because it is very effective for classification tasks. We just disregard $T_{\ell+1}$ for the simplicity of notation. It brings no additional difficulty to keep it because its Lipschitz constant is naturally equal to 1.

*Remark 3.5.* CNN are mainly used for classification purposes, contrary to feed-forward neural networks which are mainly used for regression purposes. For classification, it is highly possible that one tries to approximate (by tuning the coefficient of the networks) a function with low regularity. In this case, it is not clear what the Lipschitz constant of the function modeled by (3.1) should be. Nevertheless it is easy to imagine simple examples which are relevant for classification where the objective function is smooth (in particular if a thresholding post-processing mechanism is added). In this case smoothness could as well bring some stability in the approximation/training process. That is why we think that it is a relevant mathematical question to inquire about the Lipschitz constant of convolutional neural networks.

Let us now describe more precisely the different functions we have in mind.
The functions $g_i$, $i = 0, \dots, \ell$ can be either linear functions denoted as $g_i^{\text{lin}}(x) = W_i x + b_i$ (as for a feed-forward neural networks), or convolutional functions denoted by $g_i^{\text{conv}}$, or average-pooling functions denoted by $g_i^{\text{avg-pool}}$ (a pooling function is a particular averaging function). These functions being all linear, they do not bring conceptually new material with respect to feed-forward neural networks. Another type of non linear function $g_i^{\text{max-pool}}$, named max-pooling function, will be introduced below. All these different functions $g_i$ have the ability to change the dimension of the data: that is $g_i : \mathbb{R}^{a_i} \to \mathbb{R}^{a_{i+1}}$ where $a_{i+1} \neq a_i$ is a possibility. The functions $T_i$, $i = 1, \dots, \ell$ can be any activation function or the identity function. All these different functions $T_i$ do not change the dimension of the data: that is $T_i : \mathbb{R}^{a_i} \to \mathbb{R}^{a_{i+1}}$. As in the previous Section, the non linear activation functions $T_i$ are such that $0 \le T_i'(x) \le 1$ for almost all $x \in \mathbb{R}$ and are applied component wise.

A linear function $g_i^{\text{lin}} : \mathbb{R}^{a_i} \to \mathbb{R}^{a_{i+1}}$ is written as before:

$$g_i^{\text{lin}}(x) = W_i x + b_i.$$

A convolutional function $g_i^{\text{conv}} : \mathbb{R}^{a_i} \to \mathbb{R}^{a_{i+1}}$ operates the convolution between the input and a convolution matrix $K_i^{\text{conv}}$, with the addition of a bias

$$g_i^{\text{conv}}(x) = K_i^{\text{conv}} * \bar{x} + \bar{b}_i.$$

Since convolution operators are linear operators, there exists a matrix $W_i = W_i^{\text{conv}}$ such that

$$g_i^{\text{conv}}(x) = W_i x + b_i,$$

where $x$ and $b_i$ are vectors obtained by serializing/reindexing the vectors $\bar{x}$ and $\bar{b}_i$. With this notation, $W_i^{\text{conv}}$ is a double circulant matrix (this can be considered as a practical definition of a convolution operator in our context). The correspondance between $x$ and $\bar{x}$ is explained in more details in the following example. Let us take a vertical vector

$$(3.2) \qquad x = [x_{11} \ x_{12} \ x_{13} \ x_{21} \ x_{22} \ x_{23} \ x_{31} \ x_{32} \ x_{33}]^T.$$

A re-indexation allows to write

$$(3.3) \qquad \bar{x} = \begin{bmatrix} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{bmatrix}.$$

An average-pooling function $g_i^{\text{avg-pool}} : \mathbb{R}^{a_i} \to \mathbb{R}^{a_{i+1}}$ calculates the average of $x$ on patches which correspond to a certain multidimensional structure of the data stored in $x$. Since it is a linear operator, there exist a matrix $W_i^{\text{avg-pool}}$ such that $g_i^{\text{avg-pool}}(x) = W_i^{\text{avg-pool}} x$ as for convolutional layers. An average-pooling function with filter size (2,2) and stride 1 applied to the matrix correspond to the averaging of all $2 \times 2$ sub-matrix possible in (3.3). It yields

$$(3.4) \qquad g_i^{\text{avg-pool}}(\bar{x}) = \begin{bmatrix} \text{avg}(x_{11}, x_{12}, x_{21}, x_{22}) & \text{avg}(x_{12}, x_{13}, x_{22}, x_{23}) \\ \text{avg}(x_{21}, x_{22}, x_{31}, x_{32}) & \text{avg}(x_{22}, x_{23}, x_{32}, x_{33}) \end{bmatrix}$$

where $\text{avg}(a, b, c, d) = \frac{1}{4}(a + b + c + d)$. The very same function can of course be defined for the vector $x$ in
(3.2). One obtains

$$ (3.5) \qquad g_i^{\text{avg-pool}}(x) = \begin{bmatrix} \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} \end{bmatrix} \begin{bmatrix} x_{11} \\ x_{12} \\ x_{13} \\ x_{21} \\ x_{22} \\ x_{23} \\ x_{31} \\ x_{32} \\ x_{33} \end{bmatrix}. $$

Since the difference between (3.4) and (3.5) is just a matter of re-indexation, these two functions are the same and they are identified. The same re-indexation and identification process is used for convolution functions [16, 11].

The max-pooling operation of function is the most original one since it introduces a non linearity. Using the notation (3.3), it is written as

$$ (3.6) \qquad g_i^{\max, \text{pool}}(\bar{x}) = \begin{bmatrix} \max(x_{11}, x_{12}, x_{21}, x_{22}) & \max(x_{12}, x_{13}, x_{22}, x_{23}) \\ \max(x_{21}, x_{22}, x_{31}, x_{32}) & \max(x_{22}, x_{23}, x_{32}, x_{33}) \end{bmatrix}. $$

By construction, a CNN (3.1) is a continuous and globally Lipschitz

$$ g \in C^0(\mathbb{R}^{a_0})^{a_{\ell+1}} \cap \text{Lip}(\mathbb{R}^{a_0})^{a_{\ell+1}} $$

and the Rademacher Theorem is still valid. The rest of the Section is devoted to generalize to CNNs the bounds $K, K_*, K_1, K_3$ and $K_4$ already developed for fully-connected feed-forward neural networks. There is no additional difficulty for linear functions and for classical activation functions, that is for $g_i^{\text{lin}}$, $g_i^{\text{conv}}$, $g_i^{\text{avg-pool}}$ and $T_i$. The technical difference is for the max-pooling functions $g_i^{\text{max-pool}}$. We distinguish two approaches called the explicit approach and the implicit approach.

To be consistent with Section 2.2, we restrict our analysis to the norms $\| \cdot \|_1$ and $\| \cdot \|_{\infty}$ and we denote by $\| \cdot \|$ any of these two norms. We also highlight that we do not consider the bound $K_2$ since the matrices associated with convolutional neural networks are too large and it would be impossible to compute such bound as explained in Remark 2.9.

**3.1. Explicit approach.** The first method that we consider for the decomposition of a max-pooling function $g_i^{\text{max-pool}}$ is based on the functional identity

$$ (3.7) \qquad \begin{aligned} \max(x_1, x_2) &= \frac{1}{2}(x_1 + x_2) + \frac{1}{2}|x_1 - x_2| \\ &= \frac{1}{2}[1 \quad 1]\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + \frac{1}{2}\text{abs}\left([1 \quad -1]\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}\right) \\ &= \frac{1}{2}[1 \quad 1]\begin{bmatrix} x_1 \\ x_2 \end{bmatrix} + \frac{1}{2}z[1 \quad -1]\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \end{aligned} $$

where $z = \pm 1$. This decomposition is the generalization to max-pooling layers of the decomposition (2.6) which was considered for an activation function in the Combettes-Pesquet approach [8]. We call it explicit because the linear matrices $[1 \quad 1]$ and $[1 \quad -1]$ are explicit. The non linear part depends on just one coefficient $z$.
An interesting property is linked to this decomposition. Let us consider $u(x_1, x_2) := \max(x_1, x_2)$. We also consider $v(x_1, x_2) = x_1 + x_2$ which corresponds to the first part of the decomposition (3.7) and $w(x_1, x_2) = x_1 - x_2$ which corresponds to the second part of the decomposition. We write $u = \frac{1}{2}v + \frac{1}{2}zw$ with $z = \pm 1$. The gradient of $u$ is (almost everywhere) $\nabla u = [1 \ 0]$ or $\nabla u = [0 \ 1]$. Also $\nabla v = [1 \ 1]$ and $\nabla w = [1 \ -1]$. So one has the relation for the gradients (almost everywhere with respect to $x$)

$$ (3.8) \qquad \nabla u = \frac{1}{2} \nabla v + \frac{1}{2} z \nabla w $$

with $z = \pm 1$. Comparing with (3.7), one obtains the identity $u(x) = \nabla u(x) \cdot x$ which is the Euler identity for homogeneous function of degree one (the function $u(x_1, x_2) := \max(x_1, x_2)$ is indeed homogeneous of degree one). Then the triangular inequality yields (once again almost everywhere with respect to $x$)

$$ (3.9) \qquad \|\nabla u\|_{l^p} \le \frac{1}{2} \|\nabla v\|_{l^p} + \frac{1}{2} \|\nabla w\|_{l^p} $$

where the $l^p$ norm is the induced norm for operators, $1 \le p \le \infty$.

**LEMMA 3.6.** The $l^1$ norm is optimal for the triangular inequality (3.9), in the sense that it is an equality for $p = 1$ and a strict inequality for $1 < p \le \infty$.

*Proof.* Clearly $L(u) = (\|\nabla u\|_{l^p})_{L^\infty(\mathbb{R}^2)} = 1$. The $l^p$ norm of the operator $\nabla v = [1 \ 1]$ is evaluated with (1.1). It yields

$$ L(v) = \|\nabla v\|_{l^p} = \max_{(a,b) \ne 0} \frac{|a+b|}{(|a|^p + |b|^p)^{\frac{1}{p}}} \le \max_{(a,b) \ne 0} \frac{(|a|^p + |b|^p)^{\frac{1}{p}} (2)^{\frac{1}{q}}}{(|a|^p + |b|^p)^{\frac{1}{p}}} = 2^{\frac{1}{q}} $$

where we used a Hölder inequality and $q$ is conjugate to $p$, that is $\frac{1}{p} + \frac{1}{q} = 1$. This is actually an equality since the Hölder inequality is optimal, so $L(v) = 2^{\frac{1}{q}}$. For the same reason $L(w) = 2^{\frac{1}{q}}$. Then the triangular inequality $L(u) \le \frac{1}{2}L(v) + \frac{1}{2}L(w)$ reduces to $1 \le 2^{\frac{1}{q}} \iff 1 \le 2^{\frac{p-1}{p}}$. This inequality is an equality only for $p = 1$ which yields the claim. $\square$

This method can be applied recursively to compute a max-pooling function of any structure. A first application is as follows.

**LEMMA 3.7.** The representation (3.7) composed $n$ times yields the function $\max(x_1, x_2, \dots, x_n)$. A bound on the $l^p$ norm of the gradient is $2^{\frac{(p-1)(n-1)}{p}}$.

*Proof.* Write $\max(x_1, x_2, \dots, x_n) = \max(x_1, \max(x_2, \dots, x_n))$ and iterate Lemma 3.6. $\square$

A second application is max-pooling combined with convolution. Instead of presenting a complicated theory which will bring very little in terms of ideas, let us consider just one example which operates on either (3.2) or (3.3) with kernel size (2, 1) and stride 1 (row-wise max-pooling). So the pooling is made with every pairs in (3.3) on the same **row**. One writes

$$ (3.10) \qquad g_{\text{row}}^{\max, \text{pool}}(x) = \frac{1}{2} M_{\text{row}}^{+} x + \frac{1}{2} Z_{\text{row}}^{\max, \text{pool}}(x) M_{\text{row}}^{-} x $$
where $x \in \mathbb{R}^9$ is given in (3.2) and

$$ M_{\text{row}}^{+} = \begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \end{bmatrix}, \quad M_{\text{row}}^{-} = \begin{bmatrix} 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -1 \end{bmatrix} $$

and $Z_{\text{row}}^{\text{max-pool}}(x) \in \mathcal{M}_{6,6}(\mathbb{R})$ is a diagonal matrix with diagonal elements 1 or $-1$. The result is a vector of size 6.

Analogously, a max-pooling layer is performed on **columns** with kernel size (1, 2) and stride 1 (column-wise max-pooling) applied on the resulting vector $x \in \mathbb{R}^6$. One writes

$$ (3.11) \qquad g_{\text{col}}^{\text{max-pool}}(x) = \frac{1}{2} M_{\text{col}}^{+} x + \frac{1}{2} Z_{\text{col}}^{\text{max-pool}}(x) M_{\text{col}}^{-} x $$

where $x = [x_{11} \ x_{12} \ x_{21} \ x_{22} \ x_{31} \ x_{32}]^T \in \mathbb{R}^6$ and

$$ M_{\text{col}}^{+} = \begin{bmatrix} 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 1 \end{bmatrix}, \quad M_{\text{col}}^{-} = \begin{bmatrix} 1 & 0 & -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & 0 & -1 \end{bmatrix} $$

where $Z_{\text{col}}^{\text{max-pool}}(x) \in \mathcal{M}_{4,4}(\mathbb{R})$ is a diagonal matrix with elements 1 and $-1$. The composition yields a max-pooling function $g^{\text{max-pool}}: \mathbb{R}^9 \to \mathbb{R}^4$ (with kernel size (2, 2))
(3.12)

$$ g^{\text{max-pool}}(x) = g_{\text{col}}^{\text{max-pool}} \circ g_{\text{row}}^{\text{max-pool}}(x) = \left( \frac{1}{2} M_{\text{col}}^{+} + \frac{1}{2} Z_{\text{col}}^{\text{max-pool}}(x) M_{\text{col}}^{-} \right) \left( \frac{1}{2} M_{\text{row}}^{+} + \frac{1}{2} Z_{\text{row}}^{\text{max-pool}}(x) M_{\text{row}}^{-} \right) x $$

The gradient of the function is (almost everywhere in $x$)

$$ (3.13) \qquad \nabla g^{\text{max-pool}}(x) = \left( \frac{1}{2} M_{\text{col}}^{+} + \frac{1}{2} Z_{\text{col}}^{\text{max-pool}}(x) M_{\text{col}}^{-} \right) \left( \frac{1}{2} M_{\text{row}}^{+} + \frac{1}{2} Z_{\text{row}}^{\text{max-pool}}(x) M_{\text{row}}^{-} \right). $$

As for the example (3.8), we observe that the identity $g^{\text{max-pool}}(x) = \nabla g^{\text{max-pool}}(x) \cdot x$ holds, which comes from the fact that $g^{\text{max-pool}}(x)$ is homogeneous of degree one with respect to $x$.

To simplify the notations, we will develop the theory for (2, 2) max-pooling kernels described by (3.13). The general case will be treated in Remarks.

**Remark 3.8.** Max-pooling functions with arbitrary kernel size $(m, n)$ can be expressed as the compositions of $m-1$ row-wise max-pooling layers and $n-1$ column-wise max-pooling layers. A justification is Lemma 3.7. It can be extended to any number of dimensions. It constructs max-pooling functions with a control of the Lipschitz constant with respect to the $l^p$ norm.
Remark 3.9. In view of Lemma 3.6, it is appealing to use systematically the $l^1$ norm to get free of extra multiplicative constants for $l^p$ norms with $p > 1$. This may probably become a sensitive issue if the recursive structure is used for kernel sizes $(m, n)$ with either $m \gg 1$ or $n \gg 1$.

Let us consider $D_r(x) = T'_r \circ g_{r-1} \circ T_{r-1} \circ \dots \circ T_1 \circ g_0(x) \in \mathcal{D}_r$, which is a square diagonal matrix. In the case the function $g_r$ is linear we note $W_r = \nabla g_r$. If $g_r$ is a max-pooling function with kernel (2, 2) under the form (3.12) we note $Z_{\text{row},r}^{\text{max-pool}}(x)$ and $Z_{\text{col},r}^{\text{max-pool}}(x)$ the two diagonal matrices that show up. We do note consider kernels larger than (2, 2) for the simplicity of notation but they are immediate to treat using Remark 3.8. The main point is that these matrices $Z_{\text{row},r}^{\text{max-pool}}(x)$ and $Z_{\text{col},r}^{\text{max-pool}}(x)$ play a role similar to the matrices $D_r(x)$, so they must be treated as well with the Combettes-Pesquet trick. We note

$$ \{ Z_{\text{col},r}^{\text{max-pool}}(x), Z_{\text{row},r}^{\text{max-pool}}(x) \} \in \mathcal{Z}_r^{\text{max-pool}} $$

where $\mathcal{Z}_r^{\text{max-pool}}$ is the ensemble of all pairs of diagonal matrices of convenient size used for max-pooling layers with diagonal coefficients equal to $\pm 1$. Finally, we note

$$ Z_{\text{col}}(x) = \left( \{ Z_{\text{col},r}^{\text{max-pool}}(x), Z_{\text{row},r}^{\text{max-pool}}(x) \} \right)_r \in \mathcal{Z}^{\text{max-pool}} = \Pi_r \mathcal{Z}_r^{\text{max-pool}} $$

where the indices are restricted of course to max-pooling layers. We also denote by $W$ the collection of all remaining weight matrices for layers representing linear operators, that is $W = (\dots, W_r, \dots)$ where $\ell \ge r \ge 0$ and $g_r$ is not a max-pooling layer. The ensemble of all these weights is the vectorial spaces $\mathcal{W}$ such that $W \in \mathcal{W}$.

We now have the set of notations needed to study the Lipschitz constant of the CNN function (3.1). The chain rule yields

$$ (3.14) \qquad \nabla g = Y_{\ell}(x) D_{\ell}(x) Y_{\ell-1}(x) D_{\ell-1}(x) \dots D_{1}(x) Y_{0}(x) $$

where $D(x) = (D_{\ell}(x), \dots, D_0(x)) \in \mathcal{D}$ and $Y_r(x)$ is a matrix which is as follows. Either $Y_r(x) = W_r$, which is the layer's matrix of weights and it is constant with respect to $x$, or

$$ (3.15) \qquad Y_r(x) = \left( \frac{1}{2} M_{\text{col},r}^{+} + \frac{1}{2} Z_{\text{col},r}^{\text{max-pool}}(x) M_{\text{col},r}^{-} \right) \left( \frac{1}{2} M_{\text{row},r}^{+} + \frac{1}{2} Z_{\text{row},r}^{\text{max-pool}}(x) M_{\text{row},r}^{-} \right) $$

which depends on $x$ through the matrices $(Z_{\text{col},r}^{\text{max-pool}}(x), Z_{\text{row},r}^{\text{max-pool}}(x)) \in \mathcal{Z}_r^{\text{max-pool}}$. From the chain rule identity (3.14) one defines the Lipschitz constant

$$ L = \left( \| \nabla g(x) \|_{l^p(\mathcal{M}_{a_{\ell+1}, a_0}(\mathbb{R}))} \right)_{L^{\infty}(\mathbb{R}^{a_0})} $$

For given weight matrices $W \in \mathcal{W}$, one has by definition $L \le K$ where

$$ (3.16) \qquad K^{\text{expl}} = K^{\text{expl}}(W) = \max_{(D,Z) \in \mathcal{D} \times \mathcal{Z}^{\text{max-pool}}} \| Y_{\ell} D_{\ell} Y_{\ell-1} D_{\ell-1} \dots D_{1} Y_{0} \| $$

where either $Y_r = W_r$ or $Y_r = \left( \frac{1}{2} M_{\text{col},r}^{+} + \frac{1}{2} Z_{\text{col},r}^{\text{max-pool}} M_{\text{col},r}^{-} \right) \left( \frac{1}{2} M_{\text{row},r}^{+} + \frac{1}{2} Z_{\text{row},r}^{\text{max-pool}} M_{\text{row},r}^{-} \right)$, where the degrees of freedom are $(Z_{\text{col},r}^{\text{max-pool}}, Z_{\text{row},r}^{\text{max-pool}}) \in \mathcal{Z}_r^{\text{max-pool}}$ because the other matrices are given. This first bound is referred to as the explicit one since it is based on the explicit representation (3.7).
**3.2. Explicit CNN bounds.** The comparison between the expression (3.16) for CNNs and the expression (2.4) for feed-forward dense networks shows that the only difference is in the additional set of matrices $Z^{\max-pool}$. Apart from the inflation of notational issues due to the explicit matrices $M_{\text{col},r}^{+}$, $M_{\text{col},r}^{-}$, $M_{\text{row},r}^{+}$ and $M_{\text{row},r}^{-}$ (for relevant max-pooling layers $r$), the results are fundamentally similar to the ones in Sections 2.1 and 2.2. The main idea is explained in the following example

$$ (3.17) \qquad g = g_4 \circ T_4 \circ g_3 \circ T_3 \circ g_2 \circ T_2 \circ g_1 \circ T_1 \circ g_0. $$

We make the same assumptions as for the example (2.7) except that $g_2$ is now a max-pooling layer with kernel size (2, 2). As explained in Remark 3.2, the operator $T_3$ is now the identity. With the notation (3.15), the gradient formula (2.3) becomes

$$ (3.18) \qquad \nabla g = \frac{1}{2^5} W_4(I + Z_4)W_3 \left( (M_{\text{col},2}^{+} + Z_{\text{col},2}^{\text{max-pool}}M_{\text{col},2}^{-})(M_{\text{row},2}^{+} + Z_{\text{row},2}^{\text{max-pool}}M_{\text{row},2}^{-}) \right) (I + Z_2)W_1(I + Z_1)W_0 $$

where the identity matrices may have different sizes as in the general case. All matrices of type $Z$ are square matrices with $\pm 1$ on the diagonal, so their norm is equal to 1. Also one sees that a matrix $Z_{\text{col},2}^{\text{max-pool}}$ or $Z_{\text{row},2}^{\text{max-pool}}$ is necessarily on the left of a matrix $M_{\text{col},2}^{-}$ or $M_{\text{row},2}^{-}$.

In order to generalize this example, we define two matrices $W^{+}$ and $W^{-}$ for every linear, convolution, average-pooling or max-pooling layer. For linear layers we set $W^{+} = W^{-} = W^{\text{lin}}$. For convolutional layers we set $W^{+} = W^{-} = W^{\text{conv}}$. For average-pooling layers we set $W^{+} = W^{-} = W^{\text{avg-pool}}$. For row-wise max-pooling layers (see eq. (3.15)) we set $W^{+} = M_{\text{row}}^{+}$, $W^{-} = M_{\text{row}}^{-}$ and for column-wise max-pooling layers we set $W^{+} = M_{\text{col}}^{+}$, $W^{-} = M_{\text{col}}^{-}$. Note that max-pooling layers are fundamentally considered as two different non-linear operations, which explains why we make a distinction between **row**-like and **column**-like matrices. We set

$$ \ell' = \ell + \text{number of max-pooling layers.} $$

and we express each two-dimensional max-pooling layer as a combination of one-dimensional max-pooling layers, we obtain a neural network with $\ell' \ge \ell$ layers. We define

$$ (3.19) \qquad W_{[t,s]} = W_{t-1}^{-} W_{t-2}^{+} W_{t-3}^{+} \dots W_{s+2}^{+} W_{s+1}^{+} W_{s}^{+} $$

where the ordering is the natural one for a neural network with $\ell'$ layers.

Note that $W_{[t,s]}$ and $W_{(t,s)}$ coincide for fully connected neural networks since $W^{+} = W^{-}$ for linear layers. We consider

$$ K_{*}^{\text{expl}} := \prod_{i=0}^{\ell'} \|W_{i}^{+}\|. $$

**LEMMA 3.10.** *One has $K_{*}^{\text{expl}} \le K_{*}^{\text{expl}}$ (restricted to $l^1$ and $l^{\infty}$ norms).*

*Proof.* Consider the typical matrices of a max-pooling layer (3.15). The key property is $\|M_{\text{col},r}^{+}\| = \|M_{\text{col},r}^{-}\|$ for the $l^1$ and the $l^{\infty}$ norms. Since $\|Z_{\text{col},r}^{\text{max-pool}}\| = 1$ by construction, then $\|\frac{1}{2}M_{\text{col},r}^{+} + \frac{1}{2}Z_{\text{col},r}^{\text{max-pool}}(x)M_{\text{col},r}^{-}\| \le \|M_{\text{col},r}^{+}\| = \|W_{r}^{+}\|$. The rest of the proof is based on the sub-multiplicativity of norms. □

Now we define

$$ (3.20) \qquad K_1^{\text{expl}} = \frac{1}{2^{\ell'}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell'} \|W_{[\ell'+1, r_n]}\| \|W_{[r_n, r_{n-1}]}\| \dots \|W_{[r_2, r_1]}\| \|W_{[r_1, 0]}\|, $$
**THEOREM 3.11.** *One has $K^{\text{expl}} \le K_1^{\text{expl}} \le K_*^{\text{expl}}$.*

*Proof.* The proof of the first inequality is a direct generalization of the one of Theorem 2.4. It is based on a direct expansion of (3.14–3.15). The matrices $Z_r(x)$ (coming from $D_r(x)$), the matrices $Z_{\text{col},r}^{\text{max-pool}}$ and the matrices $Z_{\text{row},r}^{\text{max-pool}}$ divide all terms in the expansion between chunks of matrices. One can check that our notations are such that a matrix $W^{-}$ is always just before a matrix $Z$ (of any of the previous types). The rest is a matter of direct calculus.

The second inequality is proved as in Proposition 2.7, exploiting the fact that $\|W_r^+\| = \|W_r^-\|$ for one-dimensional max-pooling layers. □

**Remark 3.12** (Extension to general max-pooling layers of size $(m, n)$). Such results hold for the convolutive structure (3.15) which corresponds to the (2,2) convolutive kernel for which we gave the details of the construction in (3.10–3.12). As explained in Remark 3.8, it is sufficient to compose with additional **row-wise** and **column-wise** matrices to model two-dimensional max-pooling layers of any size. It is easy to generalize to max-pooling kernels in higher dimensions $(m, n, o, p...)$). The number $\ell'$ of terms in the formulas must be changed of course. The general formula is the sum of the numbers of terms for a feed-forward fully-connected neural networks, plus all the additional layers needed to represent the CNN.

Next we define

$$K_3^{\text{expl}} := \left\| \prod_{i=0}^{\ell'} (W_i^+)^{\text{abs}} \right\| \le K_*^{\text{expl}}.$$

**LEMMA 3.13.** *One has*

$$ (3.21) \qquad K^{\text{expl}} \le K_3^{\text{expl}} \le K_*^{\text{expl}}. $$

*Proof.* For $* \in \{\text{row, col}\}$ one has

$$ \begin{align*} \left( \frac{1}{2} M_*^+ + \frac{1}{2} Z_* M_*^- \right)_{ij}^{\text{abs}} &\le \frac{1}{2} (M_*^+)^{\text{abs}}_{ij} + \frac{1}{2} [Z_*^{\text{abs}} (M_*^-)^{\text{abs}}]_{ij} \\ &= \frac{1}{2} (M_*^+)^{\text{abs}}_{ij} + \frac{1}{2} [I (M_*^+)^{\text{abs}}]_{ij} \\ &= (M_*^+)^{\text{abs}}_{ij}. \end{align*} $$

Then, the inequalities in (3.21) follow from theorem 2.13. □

We define the following quantity

$$K_4^{\text{expl}} := \frac{1}{2^{\ell'}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell'} \left\| W_{[\ell'+1, r_n]}^{\text{abs}} W_{[r_n, r_{n-1}]}^{\text{abs}} \dots W_{[r_2, r_1]}^{\text{abs}} W_{[r_1, 0]}^{\text{abs}} \right\|.$$

**THEOREM 3.14.** *One has $K^{\text{expl}} \le K_4^{\text{expl}} \le K_1^{\text{expl}}$ and $K_4^{\text{expl}} \le K_3^{\text{expl}}$.*

*Proof.* Same proofs as in Theorem 2.14 and in Proposition 2.15. □

**3.3. Implicit approach.** The interest of what we call the implicit approach is that it simplifies the implementation for the max-pooling functions, because it is a one step technique for the modeling of $\max(x_1, \dots, x_n)$. It avoids the recursive technique used in Lemma 3.7. Moreover, this leads to bounds that are computationally more efficient because the number of terms in the sum in $K_1$ and $K_4$ grows exponentially with the number of layers.
Consider the set $\mathcal{X} = \{x_1, x_2, \dots, x_n\}$, $n \ge 2$. Let $x_M$ be the maximum element in $\mathcal{X}$ and let $\mathcal{X}^{-M} = \mathcal{X} \setminus \{x_M\}$ be the set of the remaining elements. Then, $\max(x_1, \dots, x_n)$ can be computed as

$$ (3.22) \qquad \max(x_1, \dots, x_n) = \frac{1}{2}(x_1 + x_2 + \dots + x_n) + \frac{1}{2} \left( x_M - \sum_{x_i \in \mathcal{X}^{-M}} x_i \right). $$

We consider two matrices $1 = (1, \dots, 1) \in \mathbb{R}^{1 \times n}$ and $Z_M = (-1, \dots, -1, +1, -1, \dots, -1) \in \mathbb{R}^{1 \times n}$. The matrix $Z_M$ is made only of coefficients -1 except one (at position $M$) which is equal to +1. Next, we modify the notations at the beginning of Section 3.1. Let us denote $u(x_1, \dots, x_n) := \max(x_1, \dots, x_n)$, $v(x_1, \dots, x_n) = 1 \times$ and $w(x_1, \dots, x_n) = Z_M x$. One has $u = \frac{1}{2}v + \frac{1}{2}w$. Then, almost everywhere, $\nabla u = \frac{1}{2}1 + \frac{1}{2}Z_M$, $\nabla v = 1$ and $\nabla w = Z_M$. So one has the relation for the gradients (almost everywhere with respect to $x$) $\nabla u = \frac{1}{2}\nabla v + \frac{1}{2}\nabla w$. We call this approach the *implicit one* because $Z_M$ depends on $M$.

The generalization is as follows, but once again, instead of a complicated theory, we explain how to use this decomposition for the example (3.2)-(3.3) for the max-polling function (3.6). One simply has to modify the matrices $1$ and $Z_M$. One notices that a max-pooling function like (3.6) satisfies two identities (almost everywhere), which are

$$ (3.23) \qquad g(x) = \frac{1}{2}1x + \frac{1}{2}Z_M x $$

and

$$ \nabla g(x) = \frac{1}{2}1 + \frac{1}{2}Z_M. $$

Here one has

$$ 1 = \begin{bmatrix} 1 & 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 0 & 1 & 1 \end{bmatrix}. $$

Concerning the matrix $Z_M$ there is at most 16 different possibilities (so $M < 16$) since all four lines of $M$ display coefficients equal to 0 or -1, except that one -1 per line is changed into +1. The number of different possibilities is not equal to 16 because there are some redundancies between lines. An example is

$$ Z_M = \begin{bmatrix} 1 & -1 & 0 & -1 & -1 & 0 & 0 & 0 & 0 \\ 0 & -1 & -1 & 0 & 1 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & -1 & 0 & -1 & 1 & 0 \\ 0 & 0 & 0 & 0 & -1 & -1 & 0 & 1 & -1 \end{bmatrix}. $$

A trivial property, important for further developments, holds for such matrices.

**LEMMA 3.15.** *One has $Z_M^{\text{abs}} = 1$ which does not depend on $M$. Similarly, one has $\|Z_M\| = \|1\|$.*

*Proof.* Consider the matrices $1$ and $Z_M$ of the above example. In this case, the lemma is trivial. It is the same in the general case. □

It is now simple to implement this representation of the max-pooling functions and their gradients in (3.14). Instead of (3.15) one takes for max-pooling functions only $Y_r(x) = \frac{1}{2}1_r + \frac{1}{2}Z_{M,r}$. Then one sets

$$ (3.24) \qquad K^{\text{impl}} = K^{\text{impl}}(W) = \max \|Y_r D_r Y_{r-1} D_{r-1} \dots D_1 Y_0\| $$
where either $Y_r = W_r$ or $Y_r = \frac{1}{2}1_r + \frac{1}{2}Z_{M,r}$. The maximum is taken over all possible matrices $D_r$ and all possible matrices $Z_{M,r}$.

The theoretical simplification offered by the implicit approach is clear since only one matrix $Z_{M,r}$ (with nevertheless possibly different $M$) is enough to represent the gradient of max-pooling functions for any kernels, instead of the recursive representations explained in Remark 3.8.

*Remark 3.16.* In the implicit approach, a max-pooling layer $g_r$ is expressed as the sum of a term $\frac{1}{2}1_r$ independent from the input and a term $\frac{1}{2}Z_{M,r}$ which depends on the input. A similar expansion, as the sum of a term $\frac{1}{2}M_{*,r}^+$ independent from the input and a term $\frac{1}{2}Z_{*,r}^{\max-pool}M_{*,r}^-$ depending on the input, is introduced in Section 3.1 for the one-dimensional max-pooling layer $g_*^{\max-pool}$, $* \in \{\text{row, col}\}$. The key difference between these two expansions is that the dependency on the input in the explicit approach is represented by a diagonal matrix with diagonal entries equal to 1 or -1, which is multiplied by a fixed known matrix. In the implicit approach, instead, the term depending on the input is a single rectangular matrix with norm greater than 1. This difference leads to the following alternative upper bounds.

**3.4. Implicit CNN bounds.** Since the implicit approach has a structure which is similar to the explicit one, we just write the different bounds and terms without further explanations.

Consider once again the example (3.17) where the max-pooling function $g_2$ is modeled with the implicit approach (3.23). The gradient is represented as

$$ \nabla g = \frac{1}{2^4} W_4(I + Z_4)W_3(1_2 + Z_{M,2}) (I + Z_2)W_1(I + Z_1)W_0. $$

The norm of the gradient is obtained after full expansion and use of Lemma 3.15.

To model a more general two-dimensional max-pooling layer $g_r$ we denote by $W_r^+ = 1_r$ and $W_r^- = Z_{M,r}$. We introduce the matrices

$$ (3.25) \qquad W_{\{t,s\}} = \begin{cases} W_{t-2}^+ W_{t-3}^+ \dots W_{s+1}^+ W_s^+ & \text{if } g_{t-1} \text{ is a max-pooling layer,} \\ W_{t-1}^- W_{t-2}^+ W_{t-3}^+ \dots W_{s+1}^+ W_s^+ & \text{otherwise,} \end{cases} $$

To have a compact notation, we set

$$ (3.26) \qquad R_t = \begin{cases} W_{t-1}^{-} & \text{if } g_{t-1} \text{ is a max-pooling layer,} \\ I_{t-1} & \text{otherwise,} \end{cases} $$

where $I_{t-1}$ is an identity matrix with as many rows as $W_{t-1}^-$. To prove the following bounds, it is sufficient to extend the bounds in Section 3.2 by using, for each max-pooling layer $g_r$, the quantities $(W_r^-)^{\text{abs}}$ and $\|W_r^-\|$ instead of $W_r^-$.

The following bound holds:

$$ K^{\text{impl}} \le K_*^{\text{impl}} := \prod_{i=0}^{\ell} \|W_i^+\|. $$

One also has

$$ (3.27) \qquad K^{\text{impl}} \le K_1^{\text{impl}} = \frac{1}{2^\ell} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \|W_{\{\ell+1, r_n\}}\| \|R_{r_n}\| \|W_{\{r_n, r_{n-1}\}}\| \|R_{r_{n-1}}\| \dots \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \|R_{r_2}\| \|W_{\{r_2, r_1\}}\| \|R_{r_1}\| \|W_{\{r_1, 0\}}\|. $$
One has

$$ (3.28) \qquad K^{\text{impl}} \le K_3^{\text{impl}} := \left\| \prod_{i=0}^{\ell} (W_i^+)^{\text{abs}} \right\| \le K_*^{\text{impl}}. $$

In view of our general results, the most interesting quantity is

$$ K_4^{\text{impl}} := \frac{1}{2^{\ell'}} \sum_{1 \le r_1 < r_2 < \dots < r_n \le \ell} \left\| W_{\{\ell+1, r_n\}}^{\text{abs}} R_{r_n}^{\text{abs}} W_{\{r_n, r_{n-1}\}}^{\text{abs}} R_{r_{n-1}}^{\text{abs}} \dots R_{r_2}^{\text{abs}} W_{\{r_2, r_1\}}^{\text{abs}} R_{r_1}^{\text{abs}} W_{\{r_1, 0\}}^{\text{abs}} \right\|. $$

The following inequalities hold

$$ K^{\text{impl}} \le K_4^{\text{impl}} \le K_1^{\text{impl}} \quad \text{and} \quad K_4^{\text{impl}} \le K_3^{\text{impl}}. $$

**4. Numerical results.** We provide and discuss some numerical experiments to illustrate the theoretical results presented in the previous sections. The numerical results correspond to a variety of very different test problems that come from different origins. The first test problem is the evaluation of the Lipschitz constant for a fully-connected neural network with three layers with random matrices proposed in [8]. The second test problem is particular in the sense that one has a reference solution [34, 10, 11] for the neural network approximation of $x \mapsto x^2$. It allows for sound and simple numerical comparisons, see also some developments in [12, 13]. The third test problem is a new reference solution that we propose for the function $(x, y) \mapsto xy$. The fourth problem is the now standard CNN for the MNIST problem [24]. All the numerical tests have been implemented in Tensorflow [1].

**4.1. Neural networks with random weights.** In this section, we repeat the numerical test with random matrices and weights proposed in Combettes-Pesquet [8, Example 2.1]. One considers a neural network with the following architecture

$$ (4.1) \qquad f = f_3 \circ R \circ f_2 \circ R \circ f_1. $$

Here $R: \mathbb{R} \to \mathbb{R}$ denotes the ReLU activation function, that is $R(x) = \max(0, x)$. The matrices associated with the linear layers are $W_1 \in \mathcal{M}_{10,8}(\mathbb{R})$, $W_2 \in \mathcal{M}_{6,10}(\mathbb{R})$ and $W_3 \in \mathcal{M}_{3,6}(\mathbb{R})$, their entries are i.i.d. realizations of the normal distribution $\mathcal{N}(0, 1)$. As in [8], we initialize 1000 neural networks and analyse the behaviour of the bounds. If the $l^2$ norm is used, the results presented in [8] are obtained (not reproduced here). To be coherent with the theory developed in Section 2.2 with a more general approach, we prefer to present the results computed with respect to the $l^\infty$ norm.

Note that the smaller is a bound, the closer it is to the true Lipschitz constant. The same behaviour holds for the normalized bounds shown in Table 1, where each bound is divided by $K_*$. In particular, when such ratio is close to 1, it means that the considered bound is not significantly better than the naive bound $K_*$. On the other hand, when it is close to $K/K_*$, it means that the bound is almost as sharp as the ideal (but usually incomputable) bound $K$.

The obtained statistics are reported in Table 1. Due to the small dimensionality of the considered network, it is possible to explicitly compute the ideal bound $K$ in (2.4). Coherently with the theory, we observe the inequalities $K \le K_4 \le K_3 \le K_*$ and $K_4 \le K_1 \le K_*$. Moreover $K_1$ and $K_3$ are the bounds with the lowest and highest variances (respectively). The smaller value of the ratio $K_\dagger/K_*$, $\dagger \in \{1, 2, 3, 4\}$, is obtained using the bound $K_2$.
<table><caption>**Table 1** Statistics over 1000 realizations of the network in (4.1).</caption><thead><tr><th>Statistic</th><th>K/K<sub>*</sub></th><th>K<sub>1</sub>/K<sub>*</sub></th><th>K<sub>2</sub>/K<sub>*</sub></th><th>K<sub>3</sub>/K<sub>*</sub></th><th>K<sub>4</sub>/K<sub>*</sub></th></tr></thead><tbody><tr><td>Maximum</td><td>0.2772</td><td>0.5786</td><td>0.6789</td><td>0.8023</td><td>0.4608</td></tr><tr><td>Average</td><td>0.1422</td><td>0.4539</td><td>0.3256</td><td>0.5461</td><td>0.2875</td></tr><tr><td>Minimum</td><td>0.0595</td><td>0.3703</td><td>0.1597</td><td>0.2897</td><td>0.1604</td></tr><tr><td>Standard deviation</td><td>0.0343</td><td>0.0350</td><td>0.0685</td><td>0.0813</td><td>0.0483</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Graphical representation of $g_1$, $g_2$ and $g_3$.</figcaption></figure>

**4.2. Approximation of $x^2$.** We present a neural network with an arbitrary number of layers which output is known in a simple closed form over their entire domain [34, 10, 11]. Understanding these types of models, where multiple layers and nonlinearities are present but the typical uncontrolled oscillations of deep neural networks are absent, is crucial in order to improve the deep learning mathematical theory. Let us consider the one-dimensional hat function $g: [0, 1] \to [0, 1]$

$$g(x) = \begin{cases} 2x, & x \in [0, 0.5), \\ 2(1-x), & x \in [0.5, 1]. \end{cases}$$

The function $g_r: [0, 1] \to [0, 1]$ obtained by composing $r$ times the function $g$ with itself

$$g_r(x) = \underbrace{g \circ g \circ \dots \circ g \circ g}_{r \text{ times}}.$$

The functions $g_1$, $g_2$ and $g_3$ are represented in Figure 1. The series $\sum_{r=0}^{\infty} \frac{g_r(x)}{4^r}$ converges to the function $x - x^2$ (two different proves are proposed in [11]). It is then possible to express $x^2$ as

$$(4.2) \qquad x^2 = x - \sum_{r=0}^{\infty} \frac{g_r(x)}{4^r}.$$

Since the function $g$ can be exactly expressed as a small ReLU network and the series on the right hand side of (4.2) converges exponentially fast, it is possible to efficiently approximate the squaring function with a
deep ReLU network in the following way.

* Construct a ReLU network which output coincides with the function $g$. Different choices are available. For example, if we use ReLU networks with two layers and two neurons in the internal one, $g$ can be expressed as

$$ (4.3) \qquad g(x) = 1 - R(2x - 1) - R(1 - 2x), \quad x \in [0, 1] $$

or

$$ (4.4) \qquad g(x) = R(2x) - R(4x - 2), \quad x \in [0, 1]. $$

The function in (4.3) can be represented as a two-layers ReLU neural network $g(x) = f_1 \circ R \circ f_0(x)$ with weight matrices and vectors:

$$ (4.5) \qquad W_0 = \begin{bmatrix} 2 \\ -2 \end{bmatrix}, \quad b_0 = \begin{bmatrix} -1 \\ 1 \end{bmatrix}, \quad W_1 = \begin{bmatrix} -1 & -1 \end{bmatrix}, \quad b_1 = \begin{bmatrix} 1 \end{bmatrix}. $$

Analogously, the function in (4.4) can be represented by a neural network with the same architecture but with the following weights:

$$ (4.6) \qquad W_0 = \begin{bmatrix} 2 \\ 4 \end{bmatrix}, \quad b_0 = \begin{bmatrix} 0 \\ -2 \end{bmatrix}, \quad W_1 = \begin{bmatrix} 1 & -1 \end{bmatrix}, \quad b_1 = \begin{bmatrix} 0 \end{bmatrix}. $$

* Concatenate the networks representing the function $g$ and merge each pair of consecutive linear layers into a single new linear layer. For the sake of clarity, we assume that the same representation of the function $g$ (e.g. the one in (4.3) or (4.4)) is adopted on every layer. However such a constraint is not necessary. Then the function $g_3$ can be computed as

$$ \begin{aligned} g_3 &= g \circ g \circ g \\ &= (f_1 \circ R \circ f_0) \circ (f_1 \circ R \circ f_0) \circ (f_1 \circ R \circ f_0) \\ &= f_1 \circ R \circ f_{0,1} \circ R \circ f_{0,1} \circ R \circ f_0, \end{aligned} $$

where the new layer $f_{0,1}$ is naturally defined as:

$$ \begin{aligned} f_{0,1}(x) &= f_0 \circ f_1(x) \\ &= W_0(W_1x + b_1) + b_0 \\ &= W_0W_1x + W_0b_1 + b_0 \\ &= W_{0,1}x + b_{0,1}, \end{aligned} $$

with $W_{0,1} = W_0W_1$ and $b_{0,1} = W_0b_1 + b_0$.

* Add one neuron to each layer to evaluate and add all the terms in the series in equation (4.2). It yields a collation channel [10]. The specific rows and columns that have to be added to the weight matrices depend on the chosen representation of the function $g$. The resulting neural network associated with the function $x - \sum_{r=0}^{3} \frac{g_r(x)}{4^r}$ is represented in Figure 2.

We now evaluate the bounds presented in Section 2 on this specific networks. The obtained bounds are shown in Figure 3. In the left subplot, the function $g$ is represented as in (4.3). It is possible to observe that $K_*$, $K_1$ and $K_2$ grow exponentially with the number of layers, whereas $K_3$ grows only linearly and
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Graphical representation of the network representing the function $x - \sum_{r=0}^{3} \frac{g_r(x)}{4^r}$. The blue dots and edges are associated with the construction of the function $g$. The red ones are used to store and sum $x$ and $-g_i/4^i$, $i = 1, 2, 3$. The weight associated with the dashed lines is 0.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Lipschitz bounds for networks approximating the function $x^2$. The function $g$ is represented as in (4.3) (left) or as in (4.4) (right).</figcaption></figure>

$K_4$ is exact, as highlighted in Table 2 for the first networks. In such table, the exact Lipschitz constant $L$ is computable due to the specific network structure. Instead, if the representation in (4.4) is chosen, all bounds grow exponentially, apart from $K_4$ which grows linearly. In both cases, coherently with the theory, the following inequalities hold: $K_4 \le K_1 \le K_*$, $K_4 \le K_3 \le K_*$.

Denoting by $K_\dagger^\ell$ (for $\dagger \in \{*, 1, 2, 3, 4\}$) the evaluation of one of the proposed bounds for a network with $\ell$ layers, we define the rate of growth of $K_\dagger$ as $G_\dagger^\ell = \frac{K_\dagger^{\ell+2} - K_\dagger^{\ell+1}}{K_\dagger^{\ell+1} - K_\dagger^\ell}$. This quantity is useful to understand the behavior of the growths shown in Fig 3. In fact, $K_\dagger$ grows exponentially when $G_\dagger^\ell$ is constant and larger than 1, $K_\dagger$ grows linearly when $G_\dagger^\ell = 1$, and $K_\dagger$ converges exponentially to a constant value when $G_\dagger^\ell$ is constant and smaller than 1. The behaviour of $G_\dagger^\ell$, coherent with the growths in Figure 3, is shown in Figure 4.

**4.3. Approximation of xy.** We propose a way to construct a ReLU neural network to approximate a specific polynomial function of two variables. This method is a completely different alternative to the polarization method of Yarostky for which we refer to [34, 11]. To the best of our knowledge, this construction is fully original with respect to the literature.

Let us consider the two-dimensional domain $[-1, 1]^2$ and the mesh $\mathcal{T}$ shown in Figure 5. We denote by
<table><caption>**Table 2** Exact Lipshitz constants and upper bounds for networks approximating the function $x^2$. The function $g$ is represented as in (4.3).</caption><thead><tr><th>ℓ</th><th>L</th><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th></tr></thead><tbody><tr><td>1</td><td>1.5</td><td>3.0</td><td>2.0</td><td>1.5</td><td>2.0</td><td>1.5</td></tr><tr><td>2</td><td>1.75</td><td>9.0</td><td>3.53125</td><td>3.64531</td><td>3.0</td><td>1.75</td></tr><tr><td>3</td><td>1.875</td><td>33.0</td><td>6.92187</td><td>6.45925</td><td>4.0</td><td>1.875</td></tr><tr><td>4</td><td>1.93875</td><td>129.0</td><td>14.42877</td><td>12.45882</td><td>5.0</td><td>1.93875</td></tr><tr><td>5</td><td>1.96875</td><td>513.0</td><td>30.75637</td><td>24.68184</td><td>6.0</td><td>1.96875</td></tr><tr><td>6</td><td>1.984375</td><td>2049.0</td><td>66.00227</td><td>49.24501</td><td>7.0</td><td>1.984375</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** Rates of growth for networks approximating the function $x^2$. The function $g$ is represented as in (4.3) (left) or as in (4.4) (right).</figcaption></figure>

$\varphi_{*}, * \in \{\alpha, \beta, \gamma, \delta, A, B, \dots, H, I\}$ the $\mathbb{P}_1(\mathcal{T})$ finite element basis function with value 1 at the node * and 0 at the other ones. Such basis functions are used to construct the function $\Lambda: \mathbb{R}^2 \to \mathbb{R}^2$, $\Lambda = (\Lambda_1, \Lambda_2)$, where

$$ (4.7) \qquad \begin{aligned} \Lambda_1 &= \varphi_\alpha - \varphi_\beta + \varphi_\gamma - \varphi_\delta \\ \Lambda_2 &= \varphi_D - \varphi_A + \varphi_B - \varphi_C + \varphi_F - \varphi_E - \varphi_G + \varphi_H - \varphi_I, \end{aligned} $$

Denoting by $T$ the target function $T(x, y) = xy$, we define the function $e_0$ as $e_0 = T - \frac{1}{4}T \circ \Lambda$.

**PROPOSITION 4.1.** *The function $e_0$ is linear in each triangle of the mesh $\mathcal{T}$, i.e. $e_0 \in \mathbb{P}_1(\mathcal{T})$, and it interpolates $T$ in each node of the mesh $\mathcal{T}$.*

*Proof.* Let us consider the triangle $T_{(\alpha,C,D)}$ of vertices $\alpha, C$ and $D$. Due to the local nature of the finite element basis functions, only the basis functions $\varphi_\alpha, \varphi_C$ and $\varphi_D$ are non-zero on $T_{(\alpha,C,D)}$. Therefore, $\Lambda$ can be computed as $\Lambda = (\varphi_\alpha, \varphi_D - \varphi_C)$. In particular, $\Lambda|_{T_{(\alpha,C,D)}} = (2y, 2x - 1)$. A direct evaluation of the function $e_0$ in the triangle $T_{(\alpha,C,D)}$ shows that

$$ e_0(x, y) = T(x, y) - \frac{1}{4}T \circ \Lambda(x, y) = xy - \frac{1}{4}T(2y, 2x - 1) = xy - \frac{(2y)(2x - 1)}{4} = xy - \frac{4xy - 2y}{4} = \frac{y}{2}. $$

Moreover, we observe that: $e_0(0, 0) = 0 = T(0, 0)$, $e_0(\frac{1}{2}, \frac{1}{2}) = \frac{1}{4} = T(\frac{1}{2}, \frac{1}{2})$ and $e_0(1, 0) = 0 = T(1, 0)$, i.e. $e_0$
<figure><img src="image_8.png" /><figcaption>Fig. 5 Mesh used to construct the neural network approximating the function $xy$ on $[-1, 1]^2$.</figcaption></figure>

and $T$ coincide on the vertices of the triangle. Repeating the same computation in each triangle, it can be shown that $e_0$ is linear in each triangle and $e_0 = T$ in each node of $T$. $\square$

From the definition of $e_0$, we can express $T$ as

$$ (4.8) \qquad \begin{aligned} T &= e_0 + \frac{1}{4} T \circ \Lambda = e_0 + \frac{1}{4} \left( e_0 + \frac{1}{4} T \circ \Lambda \right) \circ \Lambda \\ &= e_0 + \frac{1}{4} e_0 \circ \Lambda + \frac{1}{4^2} T \circ \Lambda \circ \Lambda = e_0 + \frac{1}{4} e_0 \circ \Lambda + \frac{1}{4^2} \left( e_0 + \frac{1}{4} T \circ \Lambda \right) \circ \Lambda \circ \Lambda \\ &= e_0 + \frac{1}{4} e_0 \circ \Lambda + \frac{1}{4^2} e_0 \circ \Lambda \circ \Lambda + \frac{1}{4^3} T \circ \Lambda \circ \Lambda \circ \Lambda = \sum_{r=0}^{\infty} \frac{1}{4^r} e_0 \circ \underbrace{\Lambda \circ \Lambda \circ \dots \circ \Lambda \circ \Lambda}_{r \text{ times}} = \sum_{r=0}^{\infty} \frac{g_r}{4^r}, \end{aligned} $$

where

$$ (4.9) \qquad g_r = e_0 \circ \underbrace{\Lambda \circ \Lambda \circ \dots \circ \Lambda \circ \Lambda}_{r \text{ times}}. $$

The derived expression is similar to (4.2), i.e. the target function is expressed as an infinite series of terms that are combinations of piecewise linear functions. Therefore, if one is able to exactly represent the functions $e_0$ and $\Lambda$ as ReLU neural networks, then it is possible to efficiently represent the target function $T$.

Let us consider the reference square $[-1, 1]^2$ and let us subdivide it in four triangles by joining the vertices with the origin. Let $\hat{\varphi}$ be the piece-wise linear function with value 1 at the origin and 0 at the vertices and which is linear in each triangle. This is our reference basis function and can be constructed, for example, as:

$$ (4.10) \qquad \hat{\varphi}(x, y) = R[2 - R(x + y) - R(x - y) - R(1 - x)], $$

or

$$ (4.11) \qquad \hat{\varphi}(x, y) = R \left[ 1 - R \left( \frac{x}{2} + \frac{y}{2} \right) - R \left( \frac{x}{2} - \frac{y}{2} \right) - R \left( -\frac{x}{2} + \frac{y}{2} \right) - R \left( -\frac{x}{2} - \frac{y}{2} \right) \right]. $$
Such functions can be expressed as ReLU neural networks with architecture $\tilde{f}_2 \circ R \circ \tilde{f}_1 \circ R \circ \tilde{f}_0$. We point out that adding the external operator $R$ (and thus the linear layer $\tilde{f}_2$) is not necessary when considering $\hat{\varphi}$ on the reference square $[-1, 1]^2$, but it is important in the neural network construction to obtain local basis functions. Using suitable linear change of variables $\mathcal{L}_*$, we construct the local basis function $\varphi_*$, $* \in \{\alpha, \beta, \gamma, \delta, B, C, D, F, H\}$ as $\varphi_* = \hat{\varphi} \circ \mathcal{L}_*$. Denoting by $\tilde{f}_0^*$ the new linear layer defined as $\tilde{f}_0^* = \tilde{f}_0 \circ \mathcal{L}_*$, each basis function $\phi_*$ can be exactly computed as $\tilde{f}_2 \circ R \circ \tilde{f}_1 \circ R \circ \tilde{f}_0^*$. Lastly, we compute $\varphi_A, \varphi_E, \varphi_G$ and $\varphi_I$ as

$$ \varphi_A(x, y) = R(R(x + y - 1)), \qquad \varphi_E(x, y) = R(R(x - y + 1)), \\ \varphi_G(x, y) = R(R(-x - y - 1)), \qquad \varphi_I(x, y) = R(R(-x + y + 1)), $$

where the double application of the ReLU operator is needed to obtain neural networks as deep as the previous ones.

Let $f_\varphi = f_2 \circ R \circ f_1 \circ R \circ f_0$ be the ReLU neural network with 2 inputs and 13 outputs obtained combining the neural networks of the 13 basis functions $\varphi_*$, $* \in \{\alpha, \beta, \gamma, \delta, A, B, \dots, H, I\}$. Let $f_3 : \mathbb{R}^{13} \to \mathbb{R}^2$ be the linear operator mapping the vector $[\varphi_\alpha, \varphi_\beta, \varphi_\gamma, \varphi_\delta, \varphi_A, \varphi_B, \dots, \varphi_H, \varphi_I]$ to the vector $[\Lambda_1, \Lambda_2]$ according to (4.7) and let $f_{3,2}$ be the linear operator $f_3 \circ f_2$. A ReLU neural network exactly representing the operator $\Lambda$ is thus $f_\Lambda = f_{3,2} \circ R \circ f_1 \circ R \circ f_0$. Analogously, let $f_4 : \mathbb{R}^{13} \to \mathbb{R}$ be the linear operator mapping the vector $[\varphi_\alpha, \varphi_\beta, \varphi_\gamma, \varphi_\delta, \varphi_A, \varphi_B, \dots, \varphi_H, \varphi_I]$ to the function $e_0$ and let $f_{4,2}$ be the linear operator $f_4 \circ f_2$. Then, a ReLU neural network exactly representing the function $e_0$ is $f_{e_0} = f_{4,2} \circ R \circ f_1 \circ R \circ f_0$.

Substituting the derived neural networks in equation (4.9), $g_r$ can be thus exactly represented as:

$$ (4.12) \qquad \begin{aligned} g_r &= e_0 \circ \underbrace{\Lambda \circ \Lambda \circ \dots \circ \Lambda \circ \Lambda}_{r \text{ times}} = f_{e_0} \circ \underbrace{f_\Lambda \circ f_\Lambda \circ \dots \circ f_\Lambda \circ f_\Lambda}_{r \text{ times}} \\ &= (f_{4,2} \circ R \circ f_1 \circ R \circ f_0) \circ \underbrace{(f_{3,2} \circ R \circ f_1 \circ R \circ f_0) \circ \dots \circ (f_{3,2} \circ R \circ f_1 \circ R \circ f_0)}_{r \text{ times}}. \end{aligned} $$

As done previously, we merge the consecutive linear layers $f_0$ and $f_{3,2}$ into a new linear layer $f_{0,3,2}$. We denote the resulting ReLU neural network by $f_{g_r}$. Finally, as in Section 4.2, the target function $T$ is approximated by truncating the series (4.8) and by enlarging the resulting neural network to store and sum all the intermediate terms involved in such a series. Adding a single neuron to each layer in Section 4.2 was enough because the target function and the approximation were always positive. Instead, since in this new case they can be both positive and negative in the domain, we need to add two neurons to propagate the information by representing the identity operator $I$ as $I(x) = R(x) - R(-x)$.

We are now able to test the bounds discussed in Section 2 on this new benchmark neural network. In Figure 6, such bounds are evaluated on neural networks of different length and constructed using the representations (4.10) and (4.11). We point out that, when (4.10) is used, the resulting neural network is a series of alternating layers of width 33 and 15, whereas when (4.11) is used, the widths of the layers are 42 and 15. In both cases, it is possible to observe that all bounds grow exponentially with respect to the number of layers and that, as expect from the theory, $K_4 \le K_1 \le K_*$ and $K_4 \le K_3 \le K_*$. Note that it is not possible to evaluate the $K_2$ bound because the layers are already too large.

In Figure 7 we show the rates of growths of the bounds shown in Figure 6. Coherently with the exponential growth of such bounds, all the rates of growths are positive and constant. We also highlight that the bounds and the rates of growth are smaller when (4.11) is used, even if it leads to a larger network. The improved performance is probably due to the symmetry of such a representation. Indeed, in Section 4.2, we observed that the bounds are sharper and the rates of growth smaller when a symmetric representation of the function $g$ is chosen.
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Lipschitz bounds for networks approximating the function $xy$. The function $\hat{\varphi}$ is represented as in (4.10) (left) or as in (4.11) (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Rates of growth for networks approximating the function $xy$. The function $\hat{\varphi}$ is represented as in (4.10) (left) or as in (4.11) (right).</figcaption></figure>

**4.4. Convolutional neural networks.** In this section, we test the bounds proposed in Section 3 on convolutional neural networks trained on the MNIST dataset. To compare the performance of the discussed bounds we train three different convolutional neural networks with the following architectures:

* Model A: Convolutional layer (5 filters) - Average-pooling layer - Convolutional layer (10 filters) - Max-pooling layer - Dense layer (20 neurons) - Dense layer (10 neurons);

* Model B: Convolutional layer (5 filters) - Max-pooling layer - Convolutional layer (10 filters) - Max-pooling layer - Dense layer (20 neurons) - Dense layer (10 neurons);

* Model C: Convolutional layer (10 filters) - Max-pooling layer - Convolutional layer (20 filters) - Max-pooling layer - Dense layer (40 neurons) - Dense layer (10 neurons).

Note that from Model A to Model B we only change the first pooling layer from an *Average-pooling layer* to a *Max-pooling layer*, whereas from Model B to Model C we only double the number of filters and neurons in the
inner layers.

We adopt the ReLU activation function everywhere, apart from the last layer where we use the softmax activation function. We always use a glorot normal initialization and we train the networks with the ADAM optimizer [22], with standard values of the hyperparameters, for 100 epochs. Each network is trained three times with different $l^2$ regularization parameters $\lambda_{reg}$ and we evaluate all the available bounds. The numerical results are shown in Tables 3, 4 and 5 for $\lambda_{reg} = 0$, $\lambda_{reg} = 1e-3$ and $\lambda_{reg} = 1e-2$ respectively.

Coherently with the theory presented in Section 3, the inequalities $K_4 \le K_1 \le K_*$ and $K_4 \le K_3 \le K_*$ hold for both the implicit and the explicit approach. The numerical results show that, in the considered cases, the implicit approach bounds are all sharper than the ones obtained with the explicit approach, both in the $l^1$ norm and in the $l^\infty$ norm. For example, note that the value of $K_4$ computed with the explicit approach is always greater than the value of $K_*$ compute with the implicit one. Lastly, we highlight that higher values of $\lambda_{reg}$ can be used to obtain neural networks with lower Lipschitz bounds.

<table><caption>**Table 3** Lipschitz bounds on convolutional neural networks trained without regularization. The accuracy on the test set is always between 97% and 98%.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Approach</th><th colspan="4">l<sup>1</sup></th><th colspan="4">l<sup>∞</sup></th></tr><tr><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th></tr></thead><tbody><tr><td rowspan="2">Model A</td><td>Explicit</td><td>1.777e6</td><td>1.347e5</td><td>2.362e5</td><td>3.119e4</td><td>2.836e7</td><td>2.381e6</td><td>8.320e6</td><td>9.236e5</td></tr><tr><td>Implicit</td><td>2.172e3</td><td>4.061e2</td><td>8.016e2</td><td><b>2.349e2</b></td><td>7.602e4</td><td>1.268e4</td><td>3.364e4</td><td><b>8.197e3</b></td></tr><tr><td rowspan="2">Model B</td><td>Explicit</td><td>3.567e8</td><td>9.815e6</td><td>7.780e6</td><td>6.364e5</td><td>9.811e8</td><td>6.063e7</td><td>2.584e8</td><td>1.752e7</td></tr><tr><td>Implicit</td><td>8.916e3</td><td>1.759e3</td><td>3.190e3</td><td><b>9.946e2</b></td><td>3.925e5</td><td>6.421e4</td><td>1.286e5</td><td><b>3.328e4</b></td></tr><tr><td rowspan="2">Model C</td><td>Explicit</td><td>8.630e9</td><td>2.179e8</td><td>1.381e8</td><td>1.176e7</td><td>2.555e10</td><td>1.478e9</td><td>4.704e9</td><td>2.861e8</td></tr><tr><td>Implicit</td><td>1.348e4</td><td>2.488e3</td><td>3.790e3</td><td><b>1.216e3</b></td><td>6.388e5</td><td>1.007e5</td><td>1.602e5</td><td><b>3.832e4</b></td></tr></tbody></table>

<table><caption>**Table 4** Lipschitz bounds on convolutional neural networks trained with $l^2$ regularization with parameter $\lambda_{reg} = 1e-3$. The accuracy on the test set is always between 97% and 98%.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Approach</th><th colspan="4">l<sup>1</sup></th><th colspan="4">l<sup>∞</sup></th></tr><tr><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th></tr></thead><tbody><tr><td rowspan="2">Model A</td><td>Explicit</td><td>4.753e5</td><td>4.102e4</td><td>5.699e4</td><td>9.077e3</td><td>9.618e6</td><td>8.795e5</td><td>1.902e6</td><td>2.333e5</td></tr><tr><td>Implicit</td><td>7.542e2</td><td>1.638e2</td><td>2.552e2</td><td><b>9.340e1</b></td><td>2.428e4</td><td>4.569e3</td><td>9.621e3</td><td><b>2.777e3</b></td></tr><tr><td rowspan="2">Model B</td><td>Explicit</td><td>1.074e8</td><td>3.996e6</td><td>1.996e6</td><td>2.608e5</td><td>2.102e8</td><td>2.177e7</td><td>6.191e7</td><td>6.672e6</td></tr><tr><td>Implicit</td><td>2.686e3</td><td>6.127e2</td><td>9.377e2</td><td><b>3.614e2</b></td><td>8.409e4</td><td>1.651e4</td><td>3.592e4</td><td><b>1.089e4</b></td></tr><tr><td rowspan="2">Model C</td><td>Explicit</td><td>2.839e9</td><td>7.436e7</td><td>3.108e7</td><td>3.226e6</td><td>5.999e9</td><td>3.607e8</td><td>1.007e9</td><td>6.863e7</td></tr><tr><td>Implicit</td><td>4.436e3</td><td>9.379e2</td><td>1.014e3</td><td><b>4.224e2</b></td><td>1.500e5</td><td>2.482e4</td><td>4.233e4</td><td><b>1.131e4</b></td></tr></tbody></table>

**5. Conclusion.** In this work we analysed and compared various upper bounds of the Lipschitz constant of deep neural networks. We considered five different bounds: the naive one ($K_*$), the ones proposed in [8] and in [32] ($K_1$ and $K_2$, respectively), and two novel bounds ($K_3$ and $K_4$). For fully-connected feed-forward neural networks, $K_*$, $K_1$ and $K_2$ are upper bounds in any $l^p$ matrix norm. However, we mainly focused on the $l^1$ and $l^\infty$ norms to provide new bounds and prove that they are sharper than the existing ones with the same computational cost, i.e. $K_3 \le K_*$ and $K_4 \le K_1$.

Two types of generalizations of such bounds to convolutional neural networks involving convolutional layers, linear layers, max-pooling layers and average-pooling layers are presented. The theoretical inequalities obtained for fully-connected feed-forward neural networks are also extended to prove that the generalization of $K_4$ is still the sharper bound among the available ones.
<table><caption>**Table 5** Lipschitz bounds on convolutional neural networks trained with $l^2$ regularization with parameter $\lambda_{\text{reg}} = 1e-2$. The accuracy on the test set is always between 95% and 97%.</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Approach</th><th colspan="4">l<sup>1</sup></th><th colspan="4">l<sup>∞</sup></th></tr><tr><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th><th>K<sub>*</sub></th><th>K<sub>1</sub></th><th>K<sub>3</sub></th><th>K<sub>4</sub></th></tr></thead><tbody><tr><td rowspan="2">Model A</td><td>Explicit</td><td>1.058e5</td><td>1.107e4</td><td>9.377e3</td><td>2.222e3</td><td>1.565e6</td><td>1.734e5</td><td>2.993e5</td><td>4.756e4</td></tr><tr><td>Implicit</td><td>1.562e2</td><td>4.721e1</td><td>5.191e1</td><td><b>2.739e3</b></td><td>3.678e3</td><td>9.021e2</td><td>1.859e3</td><td><b>6.424e2</b></td></tr><tr><td rowspan="2">Model B</td><td>Explicit</td><td>1.712e7</td><td>7.353e5</td><td>1.599e5</td><td>3.626e4</td><td>2.878e7</td><td>3.171e6</td><td>4.590e6</td><td>6.804e5</td></tr><tr><td>Implicit</td><td>4.280e2</td><td>1.423e2</td><td>1.772e2</td><td><b>1.016e2</b></td><td>1.151e4</td><td>2.823e3</td><td>6.203e3</td><td><b>2.294e3</b></td></tr><tr><td rowspan="2">Model C</td><td>Explicit</td><td>4.950e8</td><td>1.910e7</td><td>3.120e6</td><td>7.327e5</td><td>8.241e8</td><td>7.749e7</td><td>9.068e7</td><td>1.170e7</td></tr><tr><td>Implicit</td><td>7.735e2</td><td>2.337e2</td><td>1.853e2</td><td><b>1.230e2</b></td><td>2.060e4</td><td>4.570e3</td><td>6.574e3</td><td><b>2.307e3</b></td></tr></tbody></table>

We provide numerical results on fully-connected feed-forward neural networks, with random weights or approximating specific polynomials, and on convolutional neural networks with different architectures. We also proposed two ways to construct neural networks converging to the functions $x \to x^2$ and $(x, y) \to xy$ exponentially with respect to the number of layers. This is important in order to improve the theoretical understanding of deep neural networks in a simplified scenarios, where multiple layers are present but the represented function is still known in closed form.

Future extension of this works include a theoretical analysis in other norms, the development of training strategies involving the proposed bounds to obtain more stable neural networks, and the generalization of the networks converging to $x \to x^2$ and $(x, y) \to xy$ to more general functions.

**Acknowledgements.** This work is funded by PEPER/IA.

## REFERENCES

[1] M. ABADI ET AL., *TensorFlow: Large-scale machine learning on heterogeneous systems*, 2015. Software available from `tensorflow.org`.

[2] S. BERRONE, C. CANUTO, AND M. PINTORE, *Solving PDEs by variational physics-informed neural networks: an a posteriori error analysis*, Annali dell'Università di Ferrara, 68 (2022), pp. 575-595.

[3] S. BERRONE, C. CANUTO, AND M. PINTORE, *Variational physics informed neural networks: the role of quadratures and test functions*, Journal of Scientific Computing, 92 (2022), pp. 1-27.

[4] M. BOUCHEREAU, P. CHARTIER, M. LEMOU, AND F. MÉHATS, *Machine learning methods for autonomous ordinary differential equations*, 2023.

[5] J. BRADBURY, R. FROSTIG, P. HAWKINS, M. J. JOHNSON, C. LEARY, D. MACLAURIN, G. NECULA, A. PASZKE, J. VANDERPLAS, S. WANDERMAN-MILNE, AND Q. ZHANG, *JAX: composable transformations of Python+NumPy programs*, http://github.com/google/jax, (2018).

[6] S. L. BRUNTON AND J. N. KUTZ, *Singular Value Decomposition (SVD)*, Cambridge University Press, 2019, pp. 3-46.

[7] T. CHEN, J. LASSERRE, V. MAGRON, AND E. PAUWELS, *Semialgebraic optimization for lipschitz constants of relu networks*, in Proceedings of the 34th International Conference on Neural Information Processing Systems, NIPS '20, 2020.

[8] P. L. COMBETTES AND J. C. PESQUET, *Lipschitz certificates for layered network structures driven by averaged activation operators*, SIAM Journal on Mathematics of Data Science, 2 (2020), pp. 529-557.

[9] N. COUELLAN, *The Coupling Effect of Lipschitz Regularization in Neural Networks*, SN Computer Science, 2 (2021).

[10] I. DAUBECHIES, R. DEVORE, S. FOUCART, B. HANIN, AND G. PETROVA, *Nonlinear approximation and (deep) relu networks*, Constructive Approximation, 55 (2022), pp. 127-172.

[11] B. DESPRÉS, *Neural Networks and Numerical Analysis*, De Gruyter, 2022.

[12] ——, *A convergent Deep Learning algorithm for approximation of polynomials*, Comptes Rendus. Mathématique, 361 (2023), pp. 1029-1040.

[13] B. DESPRÉS AND M. ANCELLIN, *A functional equation with polynomial solutions and application to neural networks*, Comptes Rendus. Mathématique, 358 (2021), pp. 1059-1072.

[14] M. FAZLYAB, A. ROBEY, H. HASSANI, M. MORARI, AND G. J. PAPPAS, *Efficient and accurate estimation of lipschitz constants for deep neural networks*, in Proceedings of the 33rd International Conference on Neural Information Processing
Systems, 2019.

[15] J. Y. FRANCESCHI, A. FAWZI, AND O. FAWZI, Robustness of classifiers to uniform $\ell_p$ and gaussian noise, in Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, vol. 84 of Proceedings of Machine Learning Research, 2018, pp. 1280–1288.

[16] I. GOODFELLOW, Y. BENGIO, AND A. COURVILLE, Deep Learning, MIT Press, 2016. http://www.deeplearningbook.org.

[17] H. GOUK, B. PFAHRINGER, AND M. J. CREE, Regularisation of neural networks by enforcing Lipschitz continuity, Machine Learning, 110 (2021).

[18] K. GUPTA, F. KAAKAI, J. PESQUET, AND F. D. MALLIAROS, Multivariate Lipschitz Analysis of the Stability of Neural Networks, Frontiers in Signal Processing, 2 (2022).

[19] N. J. HIGHAM, Matrix Norms, Philadelphia: Soc. Industrial and Appl. Math., 1996.

[20] C. W. HUANG, A. TOUATI, P. VINCENT, G. K. DZIUGAITE, A. LACOSTE, AND A. C. COURVILLE, Stochastic neural network with kronecker flow, in International Conference on Artificial Intelligence and Statistics, 2019.

[21] D. KATSELIS, X. XIE, C. L. BECK, AND R. SRIKANT, On concentration inequalities for vector-valued lipschitz functions, Statistics & Probability Letters, 173 (2021), p. 109071.

[22] D. P. KINGMA AND J. BA, Adam: a method for stochastic optimization, arXiv preprint arXiv:1412.6980, (2014).

[23] F. LATORRE, P. T. Y. ROLLAND, AND V. CEVHER, Lipschitz constant estimation of neural networks via sparse polynomial optimization, in 8th International Conference on Learning Representations, 2020.

[24] Y. LECUN, L. BOTTOU, Y. BENGIO, AND P. HAFFNER, Gradient-based learning applied to document recognition, Proceedings of the IEEE, 86 (1998), pp. 2278–2324.

[25] S. LEE, J. LEE, AND S. PARK, Lipschitz-certifiable training with a tight outer bound, in Advances in Neural Information Processing Systems, vol. 33, 2020, pp. 16891–16902.

[26] Z. LI, F. LIU, W. YANG, S. PENG, AND J. ZHOU, A survey of convolutional neural networks: analysis, applications, and prospects, IEEE transactions on neural networks and learning systems, 33 (2021), pp. 6999–7019.

[27] H. T. D. LIU, F. WILLIAMS, A. JACOBSON, S. FIDLER, AND O. LITANY, Learning smooth neural functions via lipschitz regularization, in ACM SIGGRAPH 2022 Conference Proceedings, 2022.

[28] A. MADRY, A. MAKELOV, L. SCHMIDT, D. TSIPRAS, AND A. VLADU, Towards deep learning models resistant to adversarial attacks, in 6th International Conference on Learning Representations, ICLR, 2018.

[29] A. PASZKE ET AL., Pytorch: An imperative style, high-performance deep learning library, in Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 8024–8035.

[30] S. SHARMA, S. SHARMA, AND A. ATHAIYA, Activation functions in neural networks, Towards Data Sci, 6 (2017), pp. 310–316.

[31] C. SZEGEDY, W. ZAREMBA, I. SUTSKEVER, J. BRUNA, D. ERHAN, I. GOODFELLOW, AND R. FERGUS, Intriguing properties of neural networks, arXiv preprint arXiv:1312.6199, (2013).

[32] A. VIRMAUX AND K. SCAMAN, Lipschitz regularity of deep neural networks: analysis and efficient estimation, in Advances in Neural Information Processing Systems, vol. 31, Curran Associates, Inc., 2018.

[33] L. WENG, H. ZHANG, H. CHEN, Z. SONG, C.-J. HSIEH, L. DANIEL, D. BONING, AND I. DHILLON, Towards fast computation of certified robustness for ReLU networks, in Proceedings of the 35th International Conference on Machine Learning, vol. 80 of Proceedings of Machine Learning Research, 2018, pp. 5276–5285.

[34] D. YAROTSKY, Error bounds for approximations with deep ReLU networks, Neural Networks, 94 (2017), pp. 103–114.

[35] H. ZHANG, T.-W. WENG, P.-Y. CHEN, C.-J. HSIEH, AND L. DANIEL, Efficient neural network robustness certification with general activation functions, in Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS'18, 2018, pp. 4944–4953.