Generalization bounds for neural ordinary differential equations and deep residual networks

Introduction

Neural ordinary differential equations (neural ODEs, Chen et al., 2018) are a flexible family of neural networks used in particular to model continuous-time phenomena. Along with variants such as neural stochastic differential equations (neural SDEs, Tzen and Raginsky, 2019) and neural controlled differential equations (Kidger et al., 2020), they have been used in diverse fields such as pharmokinetics (Lu et al., 2021; Qian et al., 2021), finance (Gierjatowicz et al., 2020), and transportation (Zhou et al., 2021). We refer to Massaroli et al. (2020) for a self-contained introduction to this class of models.

Despite their empirical success, the statistical properties of neural ODEs have not yet been fully investigated. What is more, neural ODEs can be thought of as the infinite-depth limit of (properly scaled) residual neural networks (He et al., 2016a), a connection made by, e.g., E (2017); Haber and Ruthotto (2017); Lu et al. (2017). Since standard measures of statistical complexity of neural networks grow with depth (see, e.g., Bartlett et al., 2019), it is unclear why infinite-depth models, including neural ODEs, should enjoy favorable generalization properties.

To better understand this phenomenon, our goal in this paper is to study the statistical properties of a class of time-dependent neural ODEs that write

Since the parameters $\theta_{i}$ belong to an infinite-dimensional space, in practice they need to be approximated in a finite-dimensional basis of functions. For example, the residual neural networks (3) can be seen as an approximation of the neural ODEs (1) on a piecewise-constant basis of function. But more complex choices are possible, such as B-splines (Yu et al., 2022). However, the formulation (4) is agnostic from the choice of finite-dimensional approximation. This more abstract point of view is fruitful to derive generalization bounds, for at least two reasons. First, the statistical properties of the parameterized ODEs (4) only depend on the characteristics of the functions $\theta_{i}$ and not on the specifics of the approximation scheme, so it is more natural and convenient to study them at the continuous level. Second, their properties can then be transferred to any specific discretization, such as the deep residual networks (3), resulting in generalization bounds for the latter.

Regarding the characteristics of the functions $\theta_{i}$, we make the structural assumption that they are Lipschitz-continuous and uniformly bounded. This is a natural assumption to ensure that the initial value problem (4) has a unique solution in the usual sense of the Picard-Lindelöf theorem (Arnold, 1992). Remarkably, this assumption on the parameters also enables us to obtain statistical guarantees despite the fact that we are working with an infinite-dimensional set of parameters.

We provide a generalization bound for the large class of parameterized ODEs (4), which include time-dependent and time-independent neural ODEs (1) and (2). To the best of our knowledge, this is the first available bound for neural ODEs in supervised learning. By leveraging on the connection between (time-dependent) neural ODEs and deep residual networks, our approach allows us to provide a depth-independent generalization bound for the class of deep residual networks (3). The bound is precisely compared with earlier results. Our bound depends in particular on the magnitude of the difference between successive weight matrices, which is, to our knowledge, a novel way of controlling the statistical complexity of neural networks. Numerical illustration is provided to show the relationship between this quantity and the generalization ability of neural networks.

Organization of the paper.

Section 2 presents additional related work. In Section 3, we specify our class of parameterized ODEs, before stating the generalization bound for this class and for neural ODEs as a corollary. The generalization bound for residual networks is presented in Section 4 and compared to other bounds, before some numerical illustration. Section 5 concludes the paper. The proof technique is discussed in the main paper, but the core of the proofs is relegated to the Appendix.

Related work

The fields of deep learning and dynamical systems have recently benefited from sustained cross-fertilization. On the one hand, a large line of work is aimed at modeling complex continuous-time phenomena by developing specialized neural architectures. This family includes neural ODEs, but also physics-informed neural networks (Raissi et al., 2019), neural operators (Li et al., 2021) and neural flows (Biloš et al., 2021). On the other hand, successful recent advances in deep learning, such as diffusion models, are theoretically supported by ideas from differential equations (Huang et al., 2021).

Generalization for continuous-time neural networks.

Obtaining statistical guarantees for continuous-time neural networks has been the topic of a few recent works. For example, Fermanian et al. (2021) consider recurrent neural networks (RNNs), a family of neural networks handling time series, which is therefore a different setup from our work that focuses on vector-valued inputs. These authors show that a class of continuous-time RNNs can be written as input-driven ODEs, which are then proved to belong to a family of kernel methods, which entails a generalization bound. Lim et al. (2021) also show a generalization bound for ODE-like RNNs, and argue that adding stochasticity (that is, replacing ODEs with SDEs) helps with generalization. Taking another point of view, Yin et al. (2021) tackle the separate (although related) question of generalization when doing transfer learning across multiple environments. They propose a neural ODE model and provide a generalization bound in the case of a linear activation function. Closer to our setting, Hanson and Raginsky (2022) show a generalization bound for parameterized ODEs for manifold learning, which applies in particular for neural ODEs. Their proof technique bears similarities with ours, but the model and task differ from our approach. In particular, they consider stacked time-independent parameterized ODEs, while we are interested in a time-dependent formulation. Furthermore, these authors do not discuss the connection with residual networks.

Lipschitz-based generalization bounds for deep neural networks.

From a high-level perspective, our proof technique is similar to previous works (Bartlett et al., 2017; Neyshabur et al., 2018) that show generalization bounds for deep neural networks, which scale at most polynomially with depth. More precisely, these authors show that the network satisfies some Lipschitz continuity property (either with respect to the input or to the parameters), then exploit results on the statistical complexity of Lipschitz function classes. Under stronger norm constraints, these bounds can even be made depth-independent (Golowich et al., 2018). However, their approach differs from ours insofar as we consider neural ODEs and the associated family of deep neural networks, whereas they are solely interested in finite-depth neural networks. As a consequence, their hypotheses on the class of neural networks differ from ours. Section 4 develops a more thorough comparison. Similar Lipschitz-based techniques have also been applied to obtain generalization bounds for deep equilibrium networks (Pabbaraju et al., 2021). Going beyond statistical guarantees, Béthune et al. (2022) study approximation and robustness properties of Lipschitz neural networks.

Generalization bounds for parameterized ODEs

We start by recalling the usual supervised learning setup and introduce some notation in Section 3.1, before presenting our parameterized ODE model and the associated generalization bound in Section 3.2. We then apply the bound to the specific case of time-invariant neural ODEs in Section 3.3.

2 Generalization bound

for some $R_{\Theta}>0$ and $K_{\Theta}\geqslant 0$. Then, for $\theta\in\Theta$, the following Proposition, which is a consequence of the Picard-Lindelöf Theorem, shows that the mapping $x\mapsto F_{\theta}(x)$ is well-defined.

An immediate consequence of Proposition 1 is that it is legitimate to consider $\mathcal{F}_{\Theta}=\{F_{\theta},\theta\in\Theta\}$ for our model class.

Note that we consider the case where the dynamics at time $t$ are linear with respect to the parameter $\theta_{i}(t)$. Nevertheless, we emphasize that the mapping $x\mapsto F_{\theta}(x)$ remains a highly non-linear function of each $\theta_{i}(t)$. To fix ideas, this setting can be seen as analogue to working with pre-activation residual networks instead of post-activation (see He et al., 2016b, for definitions of the terminology), which is a mild modification.

Statistical analysis

Since $\Theta$ is a subset of an infinite-dimensional space, complexity measures based on the number of parameters cannot be used. Instead, our approach is to resort to Lipschitz-based complexity measures. More precisely, to bound the complexity of our model class, we propose two building blocks: we first show that the model $F_{\theta}$ is Lipschitz-continuous with respect to its parameters $\theta$. This allows us to bound the complexity of the model class depending on the complexity of the parameter class. In a second step, we assess the complexity of the class of parameters itself.

The proof, given in the Appendix, makes extensive use of Grönwall’s inequality (Pachpatte and Ames, 1997), a standard tool to obtain estimates in the theory of ODEs, in order to bound the magnitude of the solution $H_{t}$ of (5).

The next step is to assess the magnitude of the covering number of $\Theta$. Recall that, for $\varepsilon>0$, the $\varepsilon$-covering number of a metric space is the number of balls of radius $\varepsilon$ needed to completely cover the space, with possible overlaps. More formally, considering a metric space $\mathscr{M}$ and denoting by $B(x,\varepsilon)$ the ball of radius $\varepsilon$ centered at $x\in\mathscr{M}$, the $\varepsilon$-covering number of $\mathscr{M}$ is equal to $\inf\{n\geqslant 1|\exists x_{1},\dots,x_{n}\in\mathscr{M},\mathscr{M}\subseteq\bigcup_{i=1}^{n}B(x_{i},\varepsilon)\}$.

For $\varepsilon>0$, let $\mathcal{N}(\varepsilon)$ be the $\varepsilon$-covering number of $\Theta$ endowed with the distance associated to the $(1,\infty)$-norm (6). Then

Proposition 3 is a consequence of a classical result, see, e.g., Kolmogorov and Tikhomirov (1959, example 3 of paragraph 2). A self-contained proof is given in the Appendix for completeness. We also refer to Gottlieb et al. (2017) for more general results on covering numbers of Lipschitz functions.

The two propositions above and an $\varepsilon$-net argument allow to prove the first main result of our paper (where we recall that the notations are defined in Section 3.1).

Consider the class of parameterized ODEs $\mathcal{F}_{\Theta}=\{F_{\theta},\theta\in\Theta\}$, where $F_{\theta}$ is given by (5) and $\Theta$ by (7). Let $\delta>0$.

Then, for $n\geqslant 9\max(m^{-2}R_{\Theta}^{-2},1)$, with probability at least $1-\delta$,

Three terms appear in our upper bound of $\mathscr{R}(\widehat{\theta}_{n})-\widehat{\mathscr{R}}_{n}(\widehat{\theta}_{n})$. The first and the third ones are classical (see, e.g. Bach, 2023, Sections 4.4 and 4.5). On the contrary, the second term is more surprising with its convergence rate in $\mathcal{O}(n^{-1/4})$. This slower convergence rate is due to the fact that the space of parameters is infinite-dimensional. In particular, for $K_{\Theta}=0$, corresponding to a finite-dimensional space of parameters, we recover the usual $\mathcal{O}(n^{-1/2})$ convergence rate, however at the cost of considering a much more restrictive class of models. Finally, it is noteworthy that the dimensionality appearing in the bound is not the input dimension $d$ but the number of mappings $m$.

Note that this result is general and may be applied in a number of contexts that go beyond deep learning, as long as the instantaneous dependence of the ODE dynamics to the parameters is linear. One such example is the predator-prey model, describing the evolution of two populations of animals, which reads $dx_{t}=x_{t}(\alpha-\beta y_{t})dt$ and $dy_{t}=-y_{t}(\gamma-\delta x_{t})dt$, where $x_{t}$ and $y_{t}$ are real-valued variables and $\alpha$, $\beta$, $\gamma$ and $\delta$ are model parameters. This ODE falls into the framework of this section, if one were to estimate the parameters by empirical risk minimization. We refer to Deuflhard and Röblitz (2015, section 3) for other examples of parameterized biological ODE dynamics and methods for parameter identification.

Nevertheless, for the sake of brevity, we focus on applications of this result to deep learning, and more precisely to neural ODEs, which is the topic of the next section.

3 Application to neural ODEs

As explained in Section 1, parameterized ODEs include both time-dependent and time-independent neural ODEs. Since the time-independent model is more common in practice, we develop this case here and leave the time-dependent case to the reader. We thus consider the following neural ODE:

where $\sigma_{ij}(x)=\sigma(x_{j})e_{i}$. Each $\sigma_{ij}$ is itself $K_{\sigma}$-Lipschitz, hence we fall in the framework of Section 3.2. In other words, the functions $f_{i}$ of our general parameterized ODE model form a shallow neural network with pre-activation. Denote by $\|W\|_{1,1}$ the sum of the absolute values of the elements of $W$. We consider the following set of parameters, which echoes the set $\Theta$ of Section 3.2:

for some $R_{\mathcal{W}}>0$. We can then state the following result as a consequence of Theorem 1.

Consider the class of neural ODEs $\mathcal{F}_{\mathcal{W}}=\{F_{W},W\in\mathcal{W}\}$, where $F_{W}$ is given by (8) and $\mathcal{W}$ by (9). Let $\delta>0$.

Then, for $n\geqslant 9R_{\mathcal{W}}^{-1}\max(d^{-4}R_{\mathcal{W}}^{-1},1)$, with probability at least $1-\delta$,

Note that the term in $\mathcal{O}(n^{-1/4})$ from Theorem 1 is now absent. Since we consider a time-independent model, we are left with the other two terms, recovering a standard $\mathcal{O}(n^{-1/2})$ convergence rate.

Generalization bounds for deep residual networks

As highlighted in Section 1, there is a strong connection between neural ODEs and discrete residual neural networks. The previous study of the continuous case in Section 3 paves the way for deriving a generalization bound in the discrete setting of residual neural networks, which is of great interest given the pervasiveness of this architecture in modern deep learning.

We begin by presenting our model and result in Section 4.1, before detailing the comparison of our approach with other papers in Section 4.2 and giving some numerical illustration in Section 4.3.

We consider the following class of deep residual networks:

Also denoting $\|\cdot\|_{\infty}$ the element-wise maximum norm for a matrix, we consider the class of matrices

for some $R_{\mathcal{W}}>0$ and $K_{\mathcal{W}}\geqslant 0$, which is a discrete analogue of the set $\Theta$ defined by (7).

In particular, the upper bound on the difference between successive weight matrices is to our knowledge a novel way of constraining the parameters of a neural network. It corresponds to the discretization of the Lipschitz continuity of the parameters introduced in (7). By analogy, we refer to it as a constraint on the Lipschitz constant of the weights. Note that, for standard initialization schemes, the difference between two successive matrices is of the order $\mathcal{O}(1)$ and not $\mathcal{O}(1/L)$, or, in other words, $K_{\mathcal{W}}$ scales as $\mathcal{O}(L)$. This dependence of $K_{\mathcal{W}}$ on $L$ can be lifted by adding correlations across layers at initialization. For instance, one can take, for $k\in\{1,\dots,L\}$ and $i,j\in\{1,\dots,d\}$, $\mathbf{W}_{k,i,j}=\frac{1}{\sqrt{d}}f_{i,j}(\frac{k}{L})$, where $f_{i,j}$ is a smooth function, for example a Gaussian process with the RBF kernel. Such a non-i.i.d. initialization scheme is necessary for the correspondence between deep residual networks and neural ODEs to hold (Marion et al., 2022). Furthermore, Sander et al. (2022) prove that, with this initialization scheme, the constraint on the Lipschitz constant also holds for the trained network, with $K_{\mathcal{W}}$ independent of $L$. Finally, we emphasize that the following developments also hold in the case where $K_{\mathcal{W}}$ depends on $L$ (see also Section 4.2 for a related discussion).

Statistical analysis.

At first sight, a reasonable strategy would be to bound the distance between the model (10) and its limit $L\to\infty$ that is a parameterized ODE, then apply Theorem 1. This strategy is straightforward, but comes at the cost of an additional $\mathcal{O}(1/L)$ term in the generalization bound, as a consequence of the discretization error between the discrete iterations (10) and their continuous limit. For example, we refer to Fermanian et al. (2021) where this strategy is used to prove a generalization bound for discrete RNNs and where this additional error term is incurred. We follow another way by mimicking all the proof with a finite $L$. This is a longer approach but it yields a sharper result since we avoid the $\mathcal{O}(1/L)$ discretization error. The proof structure is similar to Section 3: the following two Propositions are the discrete counterparts of Propositions 2 and 3.

Let $\mathcal{N}(\varepsilon)$ be the covering number of $\mathcal{W}$ endowed with the distance associated to the $(1,1,\infty)$-norm (11). Then

The proof of Proposition 4 is a discrete analogous of Proposition 2. On the other hand, Proposition 5 can be proven as a consequence of Proposition 3, by showing the existence of an injective isometry from $\mathcal{W}$ into a set of the form (7). Equipped with these two propositions, we are now ready to state the generalization bound for our class of residual neural networks.

Consider the class of neural networks $\mathcal{F}_{\mathcal{W}}=\{F_{\mathbf{W}},\mathbf{W}\in\mathcal{W}\}$, where $F_{\mathbf{W}}$ is given by (10) and $\mathcal{W}$ by (12). Let $\delta>0$.

Then, for $n\geqslant 9R_{\mathcal{W}}^{-1}\max(d^{-4}R_{\mathcal{W}}^{-1},1)$, with probability at least $1-\delta$,

We emphasize that this result is non-asymptotic and valid for any width $d$ and depth $L$. Furthermore, the depth $L$ does not appear in the upper bound (13). This should not surprise the reader since Theorem 1 can be seen as the deep limit $L\to\infty$ of this result, hence we expect that our bound remains finite when $L\to\infty$ (otherwise the bound of Theorem 1 would be infinite). However, $L$ appears as a scaling factor in the definition of the neural network (10) and of the class of parameters (12). This is crucial for the depth independence to hold, as we will comment further on in the next section.

Furthermore, the depth independence comes at the price of a $\mathcal{O}(n^{-1/4})$ convergence rate. Note that, by taking $K_{\mathcal{W}}=0$, we obtain a generalization bound for weight-tied neural networks with a faster convergence rate in $n$, since the term in $\mathcal{O}(n^{-1/4})$ vanishes.

2 Comparison with other bounds

As announced in Section 2, we now compare Theorem 2 with the results of Bartlett et al. (2017) and Golowich et al. (2018). Beginning by Bartlett et al. (2017), we first state a slightly weaker version of their result to match our notations and facilitate comparison.

where $\widehat{\mathscr{R}}_{n}(\mathbf{W})\leqslant n^{-1}\sum_{i=1}^{n}\mathbf{1}_{F_{\mathbf{W}}(x_{i})_{y_{i}}\leqslant\gamma+\max_{j\neq y_{i}}f(x_{i})_{j}}$ and $C$ is a universal constant.

Comparing (13) and (14), we see that our bound enjoys a better dependence on the depth $L$ but a worse dependence on the width $d$. Regarding the depth, our bound (13) does not depend on $L$, whereas the bound (14) scales as $\mathcal{O}(\sqrt{L})$. This comes from the fact that we consider a smaller set of parameters (12), by adding the constraint on the Lipschitz norm of the weights. This constraint allows us to control the complexity of our class of neural networks independently of depth, as long as $K_{\mathcal{W}}$ is independent of $L$. If $K_{\mathcal{W}}$ scales as $\mathcal{O}(L)$, which is the case for i.i.d. initialization schemes, our result also features a scaling in $\mathcal{O}(\sqrt{L})$. As for the width, Bartlett et al. (2017) achieve a better dependence by a subtle covering numbers argument that takes into account the geometry induced by matrix norms. Since our paper focuses on a depth-wise analysis by leveraging the similarity between residual networks and their infinite-depth counterpart, improving the scaling of our bound with width is left for future work. Finally, note that both bounds have a similar exponential dependence in $R_{\mathcal{W}}$.

As for Golowich et al. (2018), they consider non-residual neural networks of the form $x\mapsto M_{L}\sigma(M_{L-1}\sigma(\dots\sigma(M_{1}x))).$ These authors show that the generalization error of this class scales as

where $\Pi_{F}$ is an upper-bound on the product of the Frobenius norms $\prod_{k=1}^{L}\|M_{k}\|_{F}$ and $\pi_{S}$ is a lower-bound on the product of the spectral norms $\prod_{k=1}^{L}\|M_{k}\|$. Under the assumption that both $\Pi_{F}$ and $\nicefrac{{\Pi_{F}}}{{\pi_{S}}}$ are bounded independently of $L$, their bound is indeed depth-independent, similarly to ours. Interestingly, as ours, the bound presents a $\mathcal{O}(n^{-1/4})$ convergence rate instead of the more usual $\mathcal{O}(n^{-1/2})$. However, the assumption that $\Pi_{F}$ is bounded independently of $L$ does not hold in our residual setting, since we have $M_{k}=I+\frac{1}{L}W_{k}$ and thus we can lower-bound

In our setting, it is a totally different assumption, the constraint that two successive weight matrices should be close to one another, which allows us to derive depth-independent bounds.

3 Numerical illustration

The bound of Theorem 2 features two quantities that depend on the class of neural networks, namely $R_{\mathcal{W}}$ that bounds a norm of the weight matrices and $K_{\mathcal{W}}$ that bounds the maximum difference between two successive weight matrices, i.e. the Lipschitz constant of the weights. The first one belongs to the larger class of norm-based bounds that has been extensively studied (see, e.g., Neyshabur et al., 2015). We are therefore interested in getting a better understanding of the role of the second quantity, which is much less common, in the generalization ability of deep residual networks.

To this aim, we train deep residual networks (10) (of width $d=30$ and depth $L=1000$) on MNIST. We prepend the network with an initial weight matrix to project the data $x$ from dimension $768$ to dimension $30$, and similarly postpend it with another matrix to project the output $F_{\mathbf{W}}(x)$ into dimension $10$ (i.e. the number of classes in MNIST). Finally, we consider two training settings: either the initial and final matrices are trained, or they are fixed random projections. We use the initialization scheme outlined in Section 4.1. Further experimental details are postponed to the Appendix.

We report in Figure 1(a) the generalization gap of the trained networks, that is, the difference between the test and train errors (in terms of cross entropy loss), as a function of the maximum Lipschitz constant of the weights $\sup_{0\leqslant k\leqslant L-1}(\|W_{k+1}-W_{k}\|_{\infty})$. We observe a positive correlation between these two quantities. To further analyze the relationship between the Lipschitz constant of the weights and the generalization gap, we then add the penalization term \lambda\cdot\big{(}\sum_{k=0}^{L-1}\|W_{k+1}-W_{k}\|_{F}^{2}\big{)}^{1/2} to the loss, for some $\lambda\geqslant 0$. The obtained generalization gap is reported in Figure 1(b) as a function of $\lambda$. We observe that this penalization allows to reduce the generalization gap. These two observations go in support of the fact that a smaller Lipschitz constant improves the generalization power of deep residual networks, in accordance with Theorem 2.

However, note that we were not able to obtain an improvement on the test loss by adding the penalization term. This is not all too surprising since previous work has investigated a related penalization, in terms of the Lipschitz norm of the layer sequence $(H_{k})_{0\leqslant k\leqslant L}$, and was similarly not able to report any improvement on the test loss (Kelly et al., 2020).

Finally, the proposed penalization term slightly departs from the theory that involves $\sup_{0\leqslant k\leqslant L-1}(\|W_{k+1}-W_{k}\|_{\infty})$. This is because the maximum norm is too irregular to be used in practice since, at any one step of gradient descent, it only impacts the maximum weights and not the others. As an illustration, Figure 2 shows the generalization gap when penalizing with the maximum max-norm $\sup_{0\leqslant k\leqslant L-1}(\|W_{k+1}-W_{k}\|_{\infty})$ and the $L_{2}$ norm of the max-norm \big{(}\sum_{k=0}^{L-1}\|W_{k+1}-W_{k}\|_{\infty}^{2}\big{)}^{1/2}. The factor $\lambda$ is scaled appropriately to reflect the scale difference of the penalizations. The results are mixed: the $L_{2}$ norm of the max-norm is effective contrarily to the maximum max-norm. Further investigation of the properties of these norms is left for future work.

Conclusion

We provide a generalization bound that applies to a wide range of parameterized ODEs. As a consequence, we obtain the first generalization bounds for time-independent and time-dependent neural ODEs in supervised learning tasks. By discretizing our reasoning, we also provide a bound for a class of deep residual networks. Understanding the approximation and optimization properties of this class of neural networks is left for future work. Another intriguing extension is to relax the assumption of linearity of the dynamics at time $t$ with respect to $\theta_{i}(t)$, that is, to consider a general formulation $dH_{t}=\sum_{i=1}^{m}f_{i}(H_{t},\theta_{i}(t))$. In the future, it should also be interesting to extend our results to the more involved case of neural SDEs, which have also been found to be deep limits of a large class of residual neural networks (Cohen et al., 2021; Marion et al., 2022).

Acknowledgments and Disclosure of Funding

The author is supported by a grant from Région Île-de-France, by a Google PhD Fellowship, and by MINES Paris - PSL. The author thanks Eloïse Berthier, Gérard Biau, Adeline Fermanian, Clément Mantoux, and Jean-Philippe Vert for inspiring discussions, thorough proofreading and suggestions on this paper.

References

Organization of the Appendix

Section A contains the proofs of the results of the main paper. Section B contains the details of the numerical illustrations presented in Section 4.3.

Appendix A Proofs

is locally Lipschitz-continuous with respect to its first variable and globally Lipschitz-continuous with respect to its second variable. Therefore, the existence and uniqueness of the solution of the initial value problem (5) for $t\geqslant 0$ comes as a consequence of the Picard-Lindelöf theorem (see, e.g., Luk, 2017 for a self-contained presentation and Arnold, 1992 for a textbook).

A.2 Proof of Proposition 2

For $x\in\mathcal{X}$, let $H$ be the solution of the initial value problem (5) with parameter $\theta$ and with the initial condition $H_{0}=x$. Let us first upper-bound $\|f_{i}(H_{t})\|$ for all $i\in\{1,\dots,m\}$ and $t>0$. To this aim, for $t\geqslant 0$, we have

Next, Grönwall’s inequality yields, for $t\in$,

yielding the first result of the proposition. Furthermore, for any $i\in\{1,\dots,m\}$,

Then Grönwall’s inequality implies that, for $t\in$,

since $1\leqslant K_{f}R_{\Theta}\exp(K_{f}R_{\Theta})$ because $K_{f}\geqslant 1$, $R_{\Theta}\geqslant 1$.

A.3 Proof of Proposition 3

We first prove the result for $m=1$. Let $G_{x}$ be an $\varepsilon/2K_{\Theta}$-grid of $$and$G_{y}$an$\varepsilon/2$-grid of$[-R_{\Theta},R_{\Theta}]$. Formally, we can take

Our cover consists of all functions that start at a point of $G_{y}$, are piecewise linear with kinks in $G_{x}$, where each piece has slope $+K_{\Theta}$ or $-K_{\Theta}$. Hence our cover is of size

The bounds above show that, among those two points, at least one is at distance no more than $\varepsilon/2$ from f\big{(}\frac{(k+1)\varepsilon}{K_{\Theta}}\big{)}. This shows (15) at rank $k+1$.

To conclude, take now $x\in$. There exists $k\in\{0,\dots,\lceil\frac{2K_{\Theta}}{\varepsilon}\rceil\}$ such that $x$ is at distance at most $\varepsilon/4K_{\Theta}$ from $\frac{k\varepsilon}{2K_{\Theta}}$. Again, this is clear except perhaps at the end of the interval, where it is also true since

meaning that $1$ is located between two elements of the grid $G_{x}$, showing that it is at distance at most $\varepsilon/4K_{\Theta}$ from one element of the grid. Then, we have

A.4 Proof of Theorem 1

First note that, for any $\theta\in\Theta$, $x\in\mathcal{X}$ and $y\in\mathcal{Y}$,

Now, taking $\delta>0$, a classical computation involving McDiarmid’s inequality (see, e.g., Wainwright, 2019, proof of thm 4.10) yields that, with probability at least $1-\delta$,

according to Proposition 2. The proof for the empirical risk is very similar.

Let now $\varepsilon$ > 0 and $\mathcal{N}(\varepsilon)$ be the covering number of $\Theta$ endowed with the $(1,\infty)$-norm. By Proposition 3,

Take $\theta^{(1)},\dots,\theta^{(\mathcal{N}(\varepsilon))}$ the associated cover elements. Then, for any $\theta\in\Theta$, denoting $\theta^{(i)}$ the cover element at distance at most $\varepsilon$ from $\theta$,

The remainder of the proof consists in computations to put the result in the required format. More precisely, we have

The third step is valid if $\frac{16mR_{\Theta}}{\varepsilon}\geqslant 2$. We will shortly take $\varepsilon$ to be equal to $\frac{1}{\sqrt{n}}$, thus this condition holds true under the assumption from the Theorem that $mR_{\Theta}\sqrt{n}\geqslant 3$. Hence we obtain

since $2\leqslant 2\sqrt{\log(2)}\leqslant\sqrt{2(m+1)\log(16mR_{\Theta}\sqrt{n})}$ since $16mR_{\Theta}\sqrt{n}\geqslant 2$ by the Theorem’s assumptions, and $\sqrt{2\log(4)}\leqslant 2$. We finally obtain that

by noting that $n\geqslant 9\max(m^{-2}R_{\Theta}^{-2},1)$ implies that

A.5 Proof of Corollary 1

The corollary is an immediate consequence of Theorem 1. To obtain the result, note that $m=d^{2}$, thus in particular $\sqrt{m+1}=\sqrt{d^{2}+1}\leqslant d+1$, and besides $\log(R_{\mathcal{W}}d^{2}n)\leqslant 2\log(R_{\mathcal{W}}dn)$ since $R_{\mathcal{W}}n\leqslant R_{\mathcal{W}}^{2}n^{2}$ by assumption on $n$.

A.6 Proof of Proposition 4

where the last inequality uses that the spectral norm of a matrix is upper-bounded by its $(1,1)$-norm and that $\sigma(0)=0$. As a consequence, for any $k\in\{0,\dots,L\}$,

yielding the first claim of the Proposition.

Hence, using again that the spectral norm of a matrix is upper-bounded by its $(1,1)$-norm and that $\sigma(0)=0$,

Then, dividing by $(1+K_{\sigma}\frac{R_{\mathcal{W}}}{L})^{k+1}$ and using the method of differences, we obtain that

A.7 Proof of Proposition 5

Now, take $\mathbf{W}\in\mathcal{W}$. The second property of $\phi$ implies that $\|\phi(\mathbf{W})\|_{1,\infty}\leqslant R_{\mathcal{W}}$. Moreover, each coordinate of $\phi(\mathbf{W})$ is $K_{\mathcal{W}}$-Lipschitz, since the slope of each piece of $\phi(\mathbf{W})_{i}$ is at most $K_{\mathcal{W}}$. As a consequence, $\phi(\mathbf{W})$ belongs to

Therefore $\phi(\mathcal{W})$ is a subset of $\Theta_{\mathcal{W}}$, thus its covering number is less than the one of $\Theta_{\mathcal{W}}$. Moreover, $\phi$ is clearly injective, thus we can define $\phi^{-1}$ on its image. Consider an $\varepsilon$-cover $(\theta_{1},\dots,\theta_{N})$ of $(\phi(\mathcal{W}),\|\cdot\|_{1,\infty})$. Let us show that $(\phi^{-1}(\theta_{1}),\dots,\phi^{-1}(\theta_{N}))$ is an $\varepsilon$-cover of $(\mathcal{W},\|\cdot\|_{1,1,\infty})$: take $\mathbf{W}\in\mathcal{W}$ and consider $\theta_{i}$ a cover member at distance less than $\varepsilon$ from $\phi(\mathbf{W})$. Then

where the second equality holds by linearity of $\phi$. Therefore, the covering number of $(\mathcal{W},\|\cdot\|_{1,1,\infty})$ is upper bounded by the one of $(\phi(\mathcal{W}),\|\cdot\|_{1,\infty})$, which itself is upper bounded by the one of $(\Theta_{\mathcal{W}},\|\cdot\|_{1,\infty})$, yielding the result by Proposition 3.

A.8 Proof of Theorem 2

The proof structure is the same as the one of Theorem 1, but some constants change. Similarly to (16), we obtain that, if $\frac{16d^{2}R_{\mathcal{W}}}{\varepsilon}\geqslant 2$ (which holds true for $\varepsilon=\nicefrac{{1}}{{\sqrt{n}}}$ and under the assumption of the Theorem),

for $n\geqslant 9R_{\mathcal{W}}^{-1}\max(d^{-4}R_{\mathcal{W}}^{-1},1)$. Thus

A.9 Proof of Corollary 2

where, as in the corollary, $\widehat{\mathscr{R}}_{n}(\mathbf{W})\leqslant n^{-1}\sum_{i=1}^{n}\mathbf{1}_{F_{\mathbf{W}}(x_{i})_{y_{i}}\leqslant\gamma+\max_{j\neq y_{i}}f(x_{i})_{j}}$ and $C$ is a universal constant. Let us upper bound $A(\mathbf{W})$ to conclude. On the one hand, we have

On the other hand, for any $k\in\{1,\dots,L\}$,

under the assumption that $L\geqslant R_{\mathcal{W}}$. All in all, we obtain that

Appendix B Experimental details

https://github.com/PierreMarion23/generalization-ode-resnets.

We use the following model, corresponding to model (10) with additional projections at the beginning and at the end:

We use the initialization scheme outlined in Section 4.1: we initialize, for $k\in\{1,\dots,L\}$ and $i,j\in\{1,\dots,d\}$,

where $f_{i,j}$ are independent Gaussian processes with the RBF kernel (with bandwidth equal to $0.1$). We refer to Marion et al. (2022) and Sander et al. (2022) for further discussion on this initialization scheme. However, $A$ and $B$ are initialized with a more usual scheme, namely with i.i.d. $\mathcal{N}(0,1/c)$ random variables, where $c$ denotes the number of columns of $A$ (resp. $B$).

In Figure 1(a), we repeat training $10$ times independently. Each time, we perform $30$ epochs, and compute after each epoch both the Lipschitz constant of the weights and the generalization gap. This gives $300$ pairs (Lipschitz constant, generalization gap), which each corresponds to one dot in the figure. Furthermore, we report results for two setups: when $A$ and $B$ are trained or when they are fixed random matrices.

In Figure 1(b), $A$ and $B$ are not trained. The reason is to assess the effect of the penalization on $\mathbf{W}$ for a fixed scale of $A$ and $B$. If we allow $A$ and $B$ to vary, then it is possible that the effect of the penalization might be neutralized by a scale increase of $A$ and $B$ during training.

For all experiments, we use the standard MNIST datasplit (60k training samples and 10k testing samples). We train using the cross entropy loss, mini-batches of size $128$, and the optimizer Adam (Kingma and Ba, 2015) with default parameters and a learning rate of $0.02$.

We use PyTorch (Paszke et al., 2019) and PyTorch Lightning for our experiments.

The code takes about 60 hours to run on a standard laptop (no GPU).