{'fraction_non_alphanumeric': 0.05100687826303141, 'fraction_numerical': 0.01628408055026104, 'mean_word_length': 4.245490110845468, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hamiltonian Deep Neural Networks Guaranteeing Non-vanishing Gradients by Design

Deep Neural Networks (DNNs) training can be difficult due to vanishing and exploding gradients during weight optimization through backpropagation. To address this problem, we propose a general class of Hamiltonian DNNs (H-DNNs) that stem from the discretization of continuous-time Hamiltonian systems and include several existing architectures based on ordinary differential equations. Our main result is that a broad set of H-DNNs ensures non-vanishing gradients by design for an arbitrary network depth. This is obtained by proving that, using a semi-implicit Euler discretization scheme, the backward sensitivity matrices involved in gradient computations are symplectic. We also provide an upper bound to the magnitude of sensitivity matrices, and show that exploding gradients can be either controlled through regularization or avoided for special architectures. Finally, we enable distributed implementations of backward and forward propagation algorithms in H-DNNs by characterizing appropriate sparsity constraints on the weight matrices. The good performance of H-DNNs is demonstrated on benchmark classification problems, including image classification with the MNIST dataset.


I. INTRODUCTION
Deep learning has achieved remarkable success in various fields like computer vision, speech recognition and natural language processing [1], [2]. Within the control community, there is also a growing interest in using DNNs to approximate complex controllers [3], [4]. In spite of recent progress, the training of DNNs still presents several challenges such as the occurrence of vanishing or exploding gradients during training based on gradient descent. These phenomena are related to the convergence to zero or the divergence, respectively, of the Backward Sensitivity Matrices (BSMs) arising in gradient computations through backpropagation [5]. Both situations are very critical as they imply that the learning process either stops prematurely or becomes unstable.
Heuristic methods for dealing with these problems leverage subtle weight initialization or gradient clipping [5]. More recent approaches, instead, focus on the study of DNN architectures and associated training algorithms for which exploding/vanishing gradients can be avoided or mitigated by design [6]- [16].
Recently, it has been argued that specific classes of DNNs stemming from the time discretization of Ordinary Differential Equations (ODEs) [12]- [16] are less affected by vanishing and exploding gradients. The arguments provided in [12] rely on the stability properties of the underlying continuoustime nonlinear systems for characterizing relevant behaviors of the corresponding DNNs obtained after discretization. Specifically, instability of the system results in unstable forward propagation for the DNN model, while convergence to zero of system states is related to the occurrence of vanishing gradients. This observation suggests using DNN architectures based on dynamical systems that are marginally stable, i.e. that produce bounded and non-vanishing state trajectories. An example is provided by first-order ODEs based on skewsymmetric maps, which have been used in [12], [14] for defining anti-symmetric DNNs. Another example is given by dynamical systems in the forṁ where p, q ∈ R n and H(·, ·) is a Hamiltonian function. This class of ODEs has motivated the development of Hamiltonianinspired DNNs in [12], whose effectiveness has been shown in several benchmark classification problems [12], [14], [17]. However, these approaches consider only restricted classes of skew-symmetric weight matrices or particular Hamiltonian functions, which, together with the specific structure of the dynamics in (1), limit the representation power of the resulting DNNs. Moreover, the behavior of BSMs arising in backpropagation has been analyzed only in [14], which however focuses on DNNs with identical weights in all layers and relies on hard-to-compute quantities such as kinematic eigenvalues [18], [19].
The architectures in [12], [14] and [17] are conceived for centralized implementations of forward and backward propagation. As such, they do not cope with the constraints of large networks of geographically distributed nodes with own sensing and computational capabilities, such as Internet of Things devices and multi-vehicle systems for surveillance and scanning tasks. In such scenarios, each node captures a very large stream of input data that approaches its memory, bandwidth and battery capabilities. However, in order to take system-wide optimal decisions, the measurements gathered by all nodes must be processed simultaneously, which cannot be done in a single location due to physical and computational limitations (see the recent surveys [20]- [22]). It is therefore important to develop large-scale DNN models for which the training can be distributed between physically separated end devices while guaranteeing satisfactory system-wide predictions. Furthermore, distributed DNN architectures enhance data privacy and fault tolerance [21], facilitate the learning from graph inputs [23] and enable the execution of distributed control tasks [24], [25].

A. Contributions
Leveraging general models of time-varying Hamiltonian systems [26], in this work we provide a unified framework for defining H-DNNs, which encompass anti-symmetric [12], [14] and Hamiltonian-inspired networks [12].
By considering the dynamical system associated to Hamiltonian dynamics, we prove that the norm of the BSM can never converge to 0 by design, irrespective of the network depth and the choice of the time-varying weights. This result hinges on showing that the BSM of time-varying Hamiltonian systems enjoys a property known as symplecticity. Exploiting developments in the fields of ODEs and geometrical numerical integration [27], [28], we show that symplecticity is preserved for H-DNNs stemming from Semi-Implicit Euler (S-IE) discretization 1 . Thanks to this fact, we conclude that vanishing gradients in backpropagation cannot occur.
We then analyze the phenomenon of exploding gradients. First, we construct an example showing that it cannot be avoided for general H-DNNs. Exploding gradients can, however, be kept under control by including suitable regularization terms in the training cost. Furthermore, by requiring the weights of to lie in an orthogonal group, H-DNNs can avoid exploding gradients in continuous-time; this is obtained by exploiting the results of [10] where Ode2Ode DNNs were introduced.
Finally, we show how to design distributed H-DNN architectures by imposing sparsity constraints on matrix weights. To this purpose, we provide sufficient conditions on the sparsity patterns for distributing forward-and backward-propagation algorithms over computational nodes linked by a communication network.
At a more conceptual level, our results show the potential of combining methods from system theory (as done in [29], [30]) and numerical analysis (as done in [12], [14] and [17]) for characterizing and analyzing relevant classes of deep networks.
We also highlight that H-DNNs are fundamentally different from the neural networks proposed in [31], which have the same name but are designed to learn the Hamiltonian functions of mechanical systems.
A preliminary version of this work will be presented in the L4DC conference [32]. Compared with [32], this paper is not restricted to H-DNNs stemming from forward Euler discretization. Moreover, for studying exploding/vanishing gradients we do not focus on H-DNNs with constant weights across layers and we do not rely on eigenvalue analysis. Finally, differently from [32], we consider the design of distributed H-DNN architectures.
The reminder of our paper is organized as follows. H-DNNs are defined in Section II. Their properties are analyzed in Section III and Section IV from a continuous-time and a discrete-time perspective, respectively. Numerical examples are discussed in Section V and concluding remarks are provided in Section VI.

B. Notation
We use I n to denote the identity matrix of size n×n, 0 m×n to denote the zero matrix of size m×n, 0 n to denote the square zero matrix of dimension n × n, 1 n to denote the vector of all ones of length n, and 1 m×n to denote the matrix of all ones of dimension m × n. For a vector x ∈ R n , diag(x) is the n × n diagonal matrix with the elements of x on the diagonal. For vectors x ∈ R n , y ∈ R m , we denote the vector z ∈ R m+n stacking them one after the other as z = (x, y).
We adopt the denominator layout for derivatives, that is, the derivative of y ∈ R m with respect to x ∈ R n is ∂y ∂x ∈ R n×m .

II. HAMILTONIAN DEEP NEURAL NETWORKS
This section, besides providing a short introduction to DNNs defined through the discretization of nonlinear systems, presents all the ingredients needed for the definition and implementation of H-DNNs. Throughout the paper we focus on classification tasks since they have been used as benchmarks for similar architectures (see Section V). However, the main results apply to regression tasks as well, which only require to modify the output layer of the network. We also illustrate how H-DNNs generalize several architectures recently appeared in the literature. Finally, we introduce the problem of vanishing/exploding gradients, which is analyzed in the rest of the paper.

A. DNN induced by ODE discretization
We consider the first-order nonlinear dynamical systeṁ where y(t) ∈ R n , y(0) = y 0 and θ(t) ∈ R n θ is a vector of parameters. For specifying a DNN architecture, we discretize (2) with sampling period h = T N , N ∈ N and utilize the resulting discrete-time equations for defining each of the N network layers [33]. For instance, using Forward Euler (FE) discretization, one obtains The above equation can be seen as the model of a residual neural network [1], where y j and y j+1 ∈ R n represent the input and output of layer j, respectively.
Clearly, it may be convenient to replace FE with more sophisticated discretization methods, depending on the desired structural properties of the DNN. For instance, the authors of [12] use Verlet discretization for a specific class of DNNs. Later in this work, we show that S-IE discretization will allow us to formally prove that the phenomenon of vanishing gradients cannot occur in H-DNNs.
A remarkable feature of ODE-based DNNs is that their properties can be studied by using nonlinear system theory for analyzing the continuous-time model (2). Further, this allows one to study the effect of discretization independently. As we see in later sections, such way of thinking proves effective for the case of H-DNNs.

B. From Hamiltonian dynamics to H-DNNs
The proposed neural network architectures are inspired by time-varying Hamiltonian systems [26], [34] defined aṡ where at all times and the continuously differentiable function H : In order to recover the DNNs proposed in [12], [14], [17], we consider the following Hamiltonian function whereσ : R → R is a differentiable map, applied elementwise when the argument is a matrix, and the derivative of σ(·) is called activation function σ(·). Specifically, as it is common for neural networks, we consider activation functions σ : R → R that are differentiable almost everywhere and such that for some S > 0, where σ ′ (x) denotes any sub-derivative. This assumption holds for all the more common activation functions such as tanh(·), ReLU(·), and the logistic function, among many others. Notice that and, therefore, system (4) can be rewritten aṡ The ODE (8) serves as the basis for defining H-DNNs. Indeed, as outlined in the previous subsection, a given discretization scheme for (8) naturally leads to a neural network architecture. In the context of this work, we focus on FE and S-IE discretizations. Both discretizations lead to new DNN architectures defined as follows.
H 1 -DNN: By discretizing (8) with FE we obtain the layer equation which can be interpreted as the Hamiltonian counterpart of (3). H 2 -DNN: Assume that the number of features n ∈ N is even 2 and split the feature vector at each layer j = 0, . . . , N as y j = (p j , q j ) where p j , q j ∈ R n 2 . Further assume that J j = J does not vary across layers. Then, S-IE discretization of (8) leads to the layer equation where j = 0, 1, . . . , N − 1. In general, computing the updates (p j+1 , q j+1 ) as per (10) involves solving a system of equations. We refer the reader to [36], [37] for computational aspects of deep learning with implicit layers. To make the updates (10) easily computable, one can further assume that which yields the layer equations It is easy to see that one can first compute p j+1 through (12), while q j+1 is obtained as a function of p j+1 through (13) 3 . H 1 -DNNs are motivated by their simplicity and will be compared in Section V with existing DNNs proposed in [12] and [17] using benchmark examples. However, even if system (8) is marginally stable 4 , FE discretization might introduce instability and lead to layer features y j that grow exponentially with the network depth [32]. Instead, H 2 -DNNs do not suffer from this problem because, as shown in [27], S-IE discretization preserves the marginal stability of (8). More importantly, as we show in Section IV, H 2 -DNN architectures completely prevent the phenomenon of vanishing gradients.
We conclude the section by highlighting that H 1 -and H 2 -DNNs are more general than the DNN architectures proposed in [12], [14], [17]. A precise comparison is provided in Appendix I.
To minimize the cost (14), it is common to utilize gradient descent, which can steer the parameters to a stationary point θ ⋆ of the cost such that When using gradient descent to minimize (14), at each iteration the gradient of the loss function L with respect to the parameters need to be calculated. This gradient, for parameter i of layer j is computed according to the chain rule as: It is clear from (16) that the gradient of the loss with respect to the parameters depends directly on the quantities which can be computed iteratively as follows.
Throughout the paper, we will refer to the matrix 5 α ℓ and α N are usually called weight decays.
as the BSM at layer N − j, for j = 0, . . . , N − 1, which, according to (16), allows one to compute the gradients ∂L ∂θi,N−j with respect to the parameter θ i,N −j at every layer N − j.
As shown in the next subsection, this quantity is the key in studying the phenomena of vanishing and exploding gradients.

D. Vanishing/exploding gradients
Gradient descent methods for solving (14) update the vector θ as where k is the iteration number, γ > 0 is the optimization step size and the elements of ∇ θ L are given in (16). The problem of vanishing/exploding gradients is related to the BSM (20). Indeed, when and ∂L ∂yN , and this can make the learning process unstable or cause overflow issues. Both problems are generally exacerbated when the number of layers N is large [5].
In Sections III and IV we analyze in details vanishing and exploding gradients for H-DNNs, with the goal of showing that the former cannot occur while the latter can be mitigated. To this purpose, it is convenient to first adopt the continuoustime perspective enabled by system (4).

III. CONTINUOUS-TIME ANALYSIS
In this section we analyze the properties of H-DNNs from a continuous-time point of view. First, by using backward sensitivity analysis, we derive a continuous-time representation of the backpropagation algorithm. Specifically, the continuoustime counterpart of the BSM (20) is characterized as the solution to a Linear-Time-Varying (LTV) ODE. Second, we prove that continuous-time BSMs are lower-bounded in norm by the value 1, independent of the network depth and of the choice for the time-varying weights K(t) and b(t). This will allow us to prove in Section IV-A that arbitrarily deep H 2 -DNNs enjoy non-vanishing gradients. Third, we observe that, contrary to what has been conjectured in previous work [14], the gradients of general Hamiltonian networks may explode with the network depth, even if the weights are the same in all layers. This phenomenon is shown by providing an explicit example. Fourth, in spite of the previous point, we show that it is possible to construct useful classes of H-DNNs that provably showcase non-exploding gradients. For the general H-DNN case, we propose a simple, but informative depthdependent upper-bound on the norm of the continuous-time BSMs. Our bound suggests utilizing a regularizer on the norms for the weights. Last, motivated by large-scale applications, we conclude the section by showcasing how all relevant properties of H-DNNs are naturally ported to a distributed learning setup, by appropriately constraining the sparsities for the network weights.

A. Continuous-time backward sensitivity analysis
As discussed in Section II-D, the phenomena of vanishing and exploding gradients are tightly linked to the behavior of the BSM (20). The continuous-time counterpart of (20) that we study in this section is given by .
The following Lemma, whose proof can be found in Appendix II-A, expresses the backward sensitivity dynamics of H-DNNs as the solution to an LTV ODE. Lemma 1: Given the ODE (8) associated with a H-DNN, the time evolution of the backward sensitivity matrix ∂y(T ) is given by where t ∈ [0, T ] and Accordingly, the continuous-time counterpart of (16) is given by initialized with δ(T ) = ∂L ∂y(T ) . As previously observed, the BSM ∂y(T ) ∂y(T −t) is the unique solution to the LTV ODE (23). As such, the phenomena of vanishing and exploding gradients are avoided if the corresponding LTV system is marginally stable, i.e. neither diverging nor asymptotically converging to the origin.
When θ(t) = θ for every t ∈ [0, T ], [14] has shown that A(T − t) has all eigenvalues on the imaginary axis. By adapting the results of [14] to our setting, our work [32] has further proved that A(T − t) is diagonalizable. While [14] suggests that these spectral properties of A(T − t) may lead to non vanishing and non exploding gradients when y(t) varies slowly enough, to the best of the authors' knowledge, there is no direct link between how fast A(T − t) varies and the stability of ∂y(T ) ∂y(T −t) . We will show this fact through an example later in this section.
Indeed, as opposed to linear time-invariant (LTI) systems, the stability of (23) cannot be determined solely based on the eigenvalues of the time-varying matrix A(T − t); see [38] for a simple example. As suggested in [14], one can consider the notion of kinematic eigenvalues [18], [19]. These are determined by finding a time-varying transformation that diagonalizes A(T − t). However, such transformation depends explicitly on the LTV ODE solution ∂y(T ) ∂y(T −t) . Unfortunately, this makes the computation of kinematic eigenvalues as hard as solving (23) analytically, in general.
Motivated as above, rather than studying the properties of A(T − t), in this paper we analyze the properties of ∂y(T ) directly.

B. Non-vanishing backward sensitivity matrix
Our first main result is to establish that, under the mild assumption that J(t) = J is constant across layers, ∂y(T ) Symplectic matrices are usually defined [27] by assuming n ∈ N to be an even integer and In this respect, Definition 1 provides a slightly generalized notion of symplecticity. Lemma 2: Consider a H-DNN as per (8) with ∂y(T −t) for brevity and Since ∂y(T ) ∂y(T ) = I n by definition, then ∂y(T ) Lemma 2 is an adaptation of Poincaré theorem [27], [39] to the case of the time-varying Hamiltonian functions (5) and the notion of symplecticity provided in Defintion 1.
Next, we exploit Lemma 2 to prove that any norm of the BSM is lower-bounded by the value 1. This holds irrespective of the network depth and of the choice of the time-varying weights.
Theorem 1: Consider a H-DNN as per (8) with J(t) = J for all t ∈ [0, T ], where J is any skew-symmetric matrix. Then for all t ∈ [0, T ], where · denotes any sub-multiplicative norm. Proof: We know by Lemma 2 that (25) holds. Hence, we have for all t ∈ [0, T ]. The above inequality implies the result.
The result of Theorem 1 is crucial to prove, in Section IV, that the BSM (20) is symplectic for H 2 -DNN network implementations. Hence, all H 2 -DNNs enjoy non-vanishing gradients.

C. Towards non-exploding gradients
The next question is whether may diverge to infinity in general as the depth T of the network increases. It was conjectured in [14] that, when all the weights θ(t) = θ are time-invariant and y(t) varies slowly enough, the gradients will not explode for arbitrarily deep networks. In the next subsection we show through a simple example that, unfortunately, gradients may explode even under such conditions. 1) An example of exploding gradients: We consider the twodimensional H-DNN whose forward equation is given bẏ where J = 0 −1 1 0 and ǫ ∈ R. Clearly, (27) is an instance (27) is an instance of antisymmetric network from [14] without the input term. As in [14], all the weights are time-invariant, and if |ǫ| is very small, then y(t) varies arbitrarily slowly. It can be shown that, when the weights are time-invariant, the level of the Hamiltonian remains constant for all t ∈ [0, T ], that is, H(y(t)) = H(y(0)) [12], [32]. Further, we prove here that the solution y(t) to the ODE (27) is periodic and that its period increases or decreases as the Hamiltonian energy H(y(0)) increases or decreases, respectively. The proof is reported in Appendix II-B.
In this subsection, we denote the solution to the ODE (27) at time t, initialized at time t 0 with initial condition y(t 0 ) = y 0 as s(t, t 0 , y 0 ). Lemma 3: Consider the ODE (27). The following statements hold.
ii) Let P γ,β denote the period of the ODE (27) initialized at y(0) = y 0 +γβ, where γ > 0 and β ∈ R 2 . The period increases or decreases as the Hamiltonian energy H(y(0)) increases or decreases, respectively. More precisely Lemma 3 implies that, for arbitrarily small ǫ and γ, the original and perturbed trajectories of (27) are not synchronized in phase. Let us now reason backwards in time. As T − t decreases from T to 0, assuming without loss of generality that β points towards a larger Hamiltonian sublevel set H(y 0 + γβ) > H(y 0 ), the value increases, until reaching, for a long enough time, a value of at least where the diam(·) for the sublevel set {y|H(y) = H(y 0 )} is defined as We can connect these observations to the norm of the BSM (22) as follows. For any continuous-time β ∈ R 2 that points towards a higher sublevel set of the Hamiltonian and for any γ > 0, there exists a time τ such that ∂y(T ) Since the reached value D is independent of how small γ is, we deduce that ∂y(T ) ∂y(T −τ ) β may diverge to infinity. In Appendix III, we confirm the above argument by numerical simulation for the considered H-DNN. Specifically, for any fixed value of γ > 0, reaches a maximum value of approximately D γ , which tends to infinity as γ approaches 0.
Last, notice that the provided example is valid for arbitrarily small ǫ > 0. Hence, even if y(t) varies arbitrarily slowly, exploding gradients can still occur. Furthermore, the growth of the gradients is independent of whether the weights are chosen in a time-invariant or time-varying fashion.
2) Special H-DNNs for non-exploding gradients: Despite the example outlined above, it is possible to appropriately select the time-varying weights and the activation function of a H-DNN to guarantee non-exploding gradients for arbitrary network depths. The work [10] has identified a class of neural network ODEs based on the notion of matrix flows on compact manifold, for which a proof of non-vanishing and non-exploding gradients is presented. Our observation is that, by appropriately selecting the weights, the architectures in [10] can be seen as a special kind of H-DNN.
Proposition 1: Consider the continuous-time H-DNN (8), with n ∈ N being an even number. Assume that the activation function is defined as σ(x) = |x|. Let M(t) ∈ R n×n be such that M ⊤ (t)M(t) = I n and assume for every t ∈ [0, T ], where G = 0 1 −1 0 and ⊗ denotes the Kronecker product. Then Proof: Since M ⊤ (t)M(t) = I n and G ⊤ G = I 2 , we verify that K ⊤ (t)K(t) = I n . Hence, under such choice of weights (8) is a network in the form of [10]. Second, we have that J + J ⊤ = 0 n because G + G ⊤ = 0 2 . Hence, (8) is also a H-DNN network. By Lemma 4.1 of [10], we conclude the proof.
It has been shown that the DNNs in [10] exhibit good performance for several tasks, while enjoying the property (30). We note that fulfilling orthogonality constraints during training as per [10] may considerably increase the computational burden. Nonetheless, the result of Proposition 1 indicates that the H-DNN framework is rich enough to enclose high-performing neural networks that avoid exploding gradients as well as vanishing gradients. It also opens up the question whether H-DNN structures beyond (30) with the same properties can be constructed. We leave this topic for future work.
3) Upper-bounds for general H-DNNs: Finally, we let go of any assumptions over the parameters θ(t) and consider a general H-DNN. As showcased, we expect exploding gradients as the network depth increases to infinity. Despite this fact, it remains important to derive upper-bounds on the BSM for a fixed network depth. Characterizing how gradients grow across layers can provide network design guidelines to keep them under control. We present the following result whose proof can be found in Appendix II-C.
Proposition 2 quantifies the phenomenon of exploding gradients that we have observed in the example (27). Specifically, for a general H-DNN in continuous-time with a fixed depth T ∈ R, (31) reveals that the term K(t) 2 2 J(t) 2 is crucial in keeping the gradients under control. The same is expected for discrete-time implementations.
This fact leads to the following observation. When implementing a H-DNN, it is beneficial to add the regularizer R ℓ (θ j ) = K j 2 + J j 2 in (14) to control, albeit indirectly, the magnitude of BSMs. We exploit this regularization technique for image classification with the MNIST dataset in Section V.

D. Distributed learning through H-DNNs
In this subsection, we consider utilizing sparse weight matrices to enable distributed implementations of H-DNNs. Sparsity structures in neural networks can also be used to encode prior information on relations among elements when learning graph data [23] or to perform distributed control tasks [24], [25].
First, we introduce the necessary notation for sparsity structures and binary matrices. Second, we characterize how we can design the sparsity patterns of the weight matrices K(t) and J(t) at each layer to enforce that nodes perform forward and backward propagation while complying with a fixed communication graph. Last, we observe that several sparsity choices lead to the same communication graph. Here, we deal with the continuous-time case for consistency. We extend the result to the DNNs obtained through S-IE discretization in Section IV-B. We start by introducing the necessary notation for defining and manipulating sparsity structures.
where g [i] denotes the forward update of node i and θ [i] (t) indicates the set of weights of node i at time t. Similarly, the backward propagation update of node i for every node i = 1, . . . , M must be expressed, for all t ∈ [0, T ], aṡ where h [i] denotes the backward update of node i. Next, we establish necessary and sufficient structural conditions on the weights θ(t) of a H-DNN to achieve localized forward and backward propagations. The proof is reported in Appendix II-D.
Theorem 2 (Distributed H-DNNs -continuous time): Consider a distributed H-DNN whose feature sub-vector for node i has dimension n i , and M i=1 n i = n. Assume that for each t ∈ [0, T ] we have

K(t) ∈ BlkSprs(R(t)), J(t) ∈ BlkSprs(T(t)) ,
where T(t), R(t) ∈ {0, 1} M×M are such that T(t), R(t) ≥ I M and T(t) is symmetric, as required by the skewsymmetricity of J(t). Then, the forward update (8) and the backward update (24) are computed solely based on exchanging information with the neighbors encoded in S ∈ {0, 1} M×M if, for every t ∈ [0, T ], For instance, one can describe the case of completely isolated nodes by selecting S = I M ; in this case T(t) = I M and R(t) = I M are the only admissible choices according to (34). Naturally, good global classification performance cannot be achieved unless nodes are allowed to exchange information at least with a small subset of others. As we will also see in numerical examples, it might be convenient to allow nodes to communicate according to a graph that is connected. We leave the selection of high-performing communication topologies for future work.
Observe that, for a desired S, there may be multiple choices for T(t) and R(t) that satisfy (34). One trivial choice is T(t) = S and R(t) = I M at all times. Another simple choice is T(t) = I M and R(t) such that R ⊤ (t)R(t) ≤ S. In general there are many other choices that all lead to distributed computations based on the desired communication graph, as we illustrate in the following example. In Section IV, we see how the result of Theorem 2 allows for distributed H 1 -DNN and H 2 -DNN implementations, and in Section V, we validate the effectiveness of distributed H-DNNs through numerical experiments.

IV. DISCRETE-TIME ANALYSIS
Having established useful properties of H-DNNs from a continuous-time perspective, our next goal is to preserve such properties after discretization. In this section we achieve this goal by using the S-IE discretization method.
Specifically, we show that, under mild conditions, S-IE discretization conserves the symplectic property (25) for the BSM; in turns, this allows us to show that gradients cannot vanish for any H-DNN architecture based on S-IE discretization. Firstly, we analyze the symplecticity of ∂y k+1 ∂y k . Secondly, we formally prove that the BSMs in discrete-time needed for backpropagation are lower-bounded in norm by the value 1, independent of the network depth and of the choice for the time-varying weights. Finally, we extend the analysis to the distributed learning setup in discrete time.

A. Non-vanishing gradients for H 2 -DNNs
In Section III-B, we analyzed the symplecticity of the continuous-time sensitivity matrix ∂y(T ) ∂y(T −t) . We now analyze this property directly in the discrete-time case. We prove through the following Lemma that, for (10), the matrix ∂yj+1 ∂yj is symplectic. This property will allow us to show that the discrete-time backward sensitivity matrices (20) are also symplectic and thus lower bounded in norm by the value 1.
Lemma 4: Consider the time-varying system (10) and assume that J has the block structure in (11). Then, one has The proof of Lemma 4 can be found in Appendix II-E and is built upon the result of Theorem 3.3 of Section VI in [27] and the definition of extended Hamiltonian systems [28]. Theorem 3.3 in [27] proves that the numerical flow of a time-invariant Hamiltonian system with J = 0 I −I 0 is symplectic. Moreover, by defining an extended Hamiltonian system, we can embed the study of a time-dependent Hamiltonian function into the time-independent case by defining an extended phase space of dimension n + 2 instead of n.
The property (4) allows to prove that the BSMs of H 2 -DNNs are lower bounded by 1 irrespective of the parameters and the depth of the network.

B. Distributed H 2 -DNN implementation
In the following proposition we show how to port the distributed setup developed in Section III-D to H 2 -DNNs. Moreover we analyze how the choice of the sparsity structures impacts the weights θ ⋆ verifying (15).

Proposition 3 (Distributed H 2 -DNN):
Consider the S-IE discretization of the continuous-time system (8) and split the feature vector as y j = (p j , q j ). Let the weights be in the form yielding the forward propagation updates (12)-(13) and the backward propagation updates (19). Consider nodes i = 1, . . . , M whose feature sub-vector is denoted as y j ) at every layer j. Assume that at each layer we impose where T, R j ∈ {0, 1} M×M are such that T, R j ≥ I M and T is symmetric. Then, the following two facts hold.
i) The forward update (12)- (13) and the backward update (19) are computed solely based on exchanging information with the neighbors encoded in S ∈ {0, 1} M×M if, for every j = 0, . . . , N − 1, ii) The training of this H-DNN asymptotically steers the weights θ j towards a θ ⋆ j such that for all j = 0, . . . , N −1 , where ⊙ denotes the Hadamard product andL = L + R. Proof: The proof of point i) is analogous to that of Theorem 2 by splitting the weight matrix K into submatrices K p and K q . For point ii), the regularized lossL can be rewritten as a function of only those weight entries that are not constrained to be equal to 0. The neural network training will thus steer these non-zero entries to a point where the gradient value is zero along the corresponding components. However, nothing can be said about the gradient value along the components such that R j (i, l) = 0 or T(k, v) = 0, as the corresponding weights cannot be selected to further decrease the cost.
The second point of Proposition 3 shows that, in a distributed H-DNN complying with a given structure, the achieved stationary points may be very different depending on the specific choice of weight sparsities complying with (38). Hence, the achieved performance will also strongly depend on this choice. Simple algorithms to efficiently determine optimized choices for T and R j were proposed in the context of structured feedback control; we refer the interested reader to the concept of sparsity invariance [40], [41].
To conclude, we note that a distributed H 1 -DNN based on FE discretization can also be obtained in a straightforward way. Indeed, the FE update equations (9)-(18) preserve the same structure as their continuous-time counterparts (8)- (24). We will test the performance of distributed H-DNNs in the following Section V.

V. NUMERICAL EXPERIMENTS
In this section, we demonstrate the potential of H-DNNs on various classification benchmarks, including the MNIST dataset 6 . Our first goal is to show that H-DNNs are expressive enough to achieve state-of-art performance on these tasks, even when distributed architectures are used. Then, we validate our main theoretical results by showing that gradients do not vanish despite considering very deep H 2 -DNNs. Instead, when using the same data, standard multilayer perceptron networks do suffer from this problem, which causes an early termination of the learning process.

A. Binary classification examples
We consider two benchmark classification problems from [12] with two categories and two features. First, we show that H 1 -and H 2 -DNNs perform as well as the networks MS i -DNNs, i = 1, 2, 3 given in Appendix I and introduced in [12], [14] and [17]. Second, we test a distributed H 2 -DNN architecture, showing that it can achieve excellent accuracy despite local computations. In all cases, we optimize over the weights K j and b j with for j = 0, . . . , N − 1.
As outlined in Section II-C, we complement the DNNs with a soft-max output layer y N +1 = f N (y N , θ N ). In (14), we use a standard cross-entropy loss function L [5]. The optimization problem is solved by using the Adam algorithm with minibatches; we refer to Appendix IV-A for details on the implementation. 1) Comparison with existing networks: We consider the "Swiss roll" and the "Double moons" datasets shown in Figure  1. By performing feature augmentation [35], we use input feature vectors given by (y k 0 , 0 2×1 ) ∈ R 4 where y k 0 ∈ R 2 , k = 1, . . . , s are the input datapoints.
In Table I, we present the classification accuracies over test sets when using MS-and H-DNNs with different number of layers. We also give the number of parameters per layer of each network. It can be seen that, for a fixed number of layers, the performances of H 1 -DNN and H 2 -DNN are similar or better than the other networks. This can be motivated by the fact that the new architectures are more expressive than MS i -DNNs, since they have, in general, more parameters per layer. Although not shown, H-and MS-DNNs with same number of parameters have similar performance.
2) Distributed vs. centralized training: We test the effectiveness of H 2 -DNNs for a distributed learning task under constraints (11). We assume that 8 nodes have access to their own local features y [i] ∈ R 2 , i = 1, . . . , 8. Moreover, they communicate with a small set of neighbors y [j] according to the connected graph in Figure 2, whose adjacency matrix is denoted with S. Therefore, each node computes local forward and backward updates only based on the features known to itself and to its first-and second-order neighbors. for each layer j. The first three columns represent existing architectures (MS i -DNN) while the two last columns provide the results for the new H i -DNNs. The first two best accuracies in each row are in bold. The last row provides the number of parameters per layer of each network. (a) "Swiss roll" -64 layers (b) "Double moons" -4 layers In order to comply with the communication constraints, we exploit (38) in Proposition 3, which is verified, for j = 0, . . . , N , by T = I 8 and R j equal to the adjacency matrix of the subgraph including only the blue edges in Figure 2. Since J j in the form of (11) is chosen equal to (39), then X ∈ BlkSprs(T). Accordingly, the weights K p,j , K q,j are constrained to lie in BlkSprs(R j ).
Then, we initialize the augmented input feature vectorỹ 0 = (p 0 , q 0 ) ∈ R 16 with zeros, except for the two entries p [1] 0 = y k 0,0 and q [5] 0 = y k 0,1 . Table II reports the prediction accuracies for the distributed and centralized learning settings when using the "Swiss roll" and "Double circles" datasets shown in Figure 3. Moreover, the colored regions in Figure 3 provide an illustration of the   predictive power of the distributed H 2 -DNNs with four layers. From Table II, even if the centralized architectures perform better for shallow networks (e.g. 2 layers), it is worth noticing that the distributed architectures use less than half of the parameters. Moreover, if the networks are deep enough, the classification task is successfully solved in both settings.

B. Experiments with the MNIST dataset
We evaluate our methods on a more complex example: the image classification benchmark MNIST 7 . The dataset consists of 28 × 28 digital images in gray scale of hand-written digits from 0 to 9 with their corresponding labels. It contains 60,000 training examples and 10,000 test examples.
Following [12], we use a network architecture consisting of a convolutional layer followed by a Hamiltonian DNN and an output classification layer. The convolutional layer is a linear transformation that expands the data from 1 to 8 channels, and the network output is a vector in R 10 representing the probabilities of an image to belong to each of the 10 classes.
We compare the performance of MS 1 -DNNs and H 1 -DNNs 8 with 2, 4, 6 and 8 layers. We set h = 0.4 for MS 1 -DNNs and h = 0.5 for H 1 -DNNs. The implementation details can be found in the Appendix IV-B.
In Table III, we summarize the train and test accuracies of the considered networks. The first row of the table provides, as a baseline, the results obtained when omitting the Hamiltonian DNN block, i.e., when using only a convolutional layer followed by the output layer. We observe that both MS 1 -DNN and H 1 -DNN achieve similar performance. Note that, while the training errors are almost zero, the test errors are reduced when increasing the number of layers, hence showing the benefit of using deeper networks. Moreover, these results are in line with test accuracies obtained when using standard convolutional layers instead of H-DNNs [17].

C. Gradient analysis
Our next aim is to provide a numerical validation of the main property of H-DNN networks: the absence of vanishing gradients. We consider the "Double moons" dataset shown in Figure 1b and analyze the norm of BSMs during the training of a H 2 -DNN and a fully connected multilayer perceptron (MLP) network 9 . Figure 4a displays the BSM norms during the 960 iterations of the training phase of a 32-layer H 2 -DNN and Figure 4b presents the same quantities for a MLP network with 8 and 32 layers. While the H 2 -DNN and the 8-layer MLP network achieve good performance at the end of the training (100% and 99.7% accuracy over the test set, respectively), the 32-layer MLP network fails to classify the input data (50% accuracy). Figure 4a validates the results of Theorem 3 since no BSM norm is smaller than 1 at any time 10 . However, this is not the case for MLP networks where gradient norms can be very small and may prevent the network to succeed in the training. Note that the main inconvenient for the 32-layer MLP network to achieve good performance is that after 400 iterations, only a few layers still show a gradient norm different from zero i.e. ∂yN ∂y N −ℓ ≈ 0 for ℓ = 7, . . . , 32.
In addition, it can be seen from Figure 4a that the BSM norms do not explode and remain bounded as This is in line with Proposition 2, where we show that, in each iteration of the training phase, BSM norms are upper bounded by a quantity depending on the network depth and the parameters K j and b j of that specific iteration. 9 The layer equation of the MLP network is given by y j+1 = σ(K j y j + b j ), with activation function σ(·) = tanh(·) and trainable parameters K j and b j . We use the same implementation as described in Appendix IV-A. 10 We highlight that L(θ) achieves a stationary point θ ⋆ since ∇ θ L(θ ⋆ ) = 0 in approximately 350 iterations, obtaining a 100% accuracy over the test set.

VI. CONCLUSION
In this paper, we proposed a class of H-DNNs obtained from the time-discretization of Hamiltonian dynamics. We proved that H-DNNs stemming from S-IE discretization do not suffer from vanishing gradients and also provided methods to control the growth of associated BSMs. We further derived sufficient structural conditions on weight matrices to enable distributed learning of H-DNNs. Although we limited our analysis to S-IE discretization, one can leverage the rich literature on symplectic integration [27] for defining even broader classes of H-DNNs. This avenue will be explored in future research. Moreover, we will study the use of distributed H-DNNs architectures for parametrizing distributed controllers for multi-agent systems. tion method to discretize where y, z ∈ R n 2 , obtaining the layer equations Note that (40) is an instance of H 2 -DNN when assuming K j to be invertible and setting J j K ⊤ j = I n , and K j = 0 K0,j −K ⊤ 0,j 0 for all j = 0, . . . , N − 1. MS 2 -DNN: In [12], [14], the authors propose to use FE to discretizeẏ where K(t) is skew-symmetric ∀t ∈ [0, T ], obtaining the following layer equation In this case, (41) is an instance of H 1 -DNN when assuming K j to be invertible and setting J j K ⊤ j = I n and K j = −K ⊤ j for all j = 0, . . . , N − 1. MS 3 -DNN: In [17], the authors propose to use Verlet integration method to discretize where y, z ∈ R n 2 , obtaining the layer equations Note that (42) is an instance of H 2 -DNN when setting In [12] and [17], the DNNs MS 1 and MS 3 have been called Hamiltonian-inspired in view of their similarities with Hamiltonian models, although a precise Hamiltonian function for the corresponding ODE has not been provided. Moreover, note that the Verlet discretization used coincides with S-IE.
We highlight that a necessary condition for the skewsymmetric n × n matrix K j to be invertible is that the size n of input features is even 11 . If n is odd, however, one can perform input-feature augmentation by adding an extra state initialized at zero to satisfy the previous condition [35].
It can be verified that the value sech 2 (x)x − tanh(x) is negative for every x > 0. We conclude that the angular time derivative decreases as the perturbation γ increases. Since the above holds for any orientation of the Cartesian axes, we conclude that the period increases as γ increases. Analogous reasoning holds if β points towards a lower sublevel set; in this case, the period decreases as γ increases.

C. Proof of Proposition 2
We will first state the following Lemma needed for proving Theorem 2.
Proof: Consider x ∈ R n where the i-th entry of x is denoted as x i for i = 1, . . . , n. Since Ax = n i=1 x i a i , then, It is easy to prove using the Cauchy-Schwarz inequality that 1 √ n x 1 ≤ x 2 . Then, we conclude: We are now ready to prove Theorem 2.
Proof: Consider the backward gradient dynamics ODE (23) and subdivide ∂y(T ) ∂y(T −t) into columns as per ∂y(T ) ∂y(T −t) = z 1 (t) z 2 (t) . . . z n (t) . Then, (23) is equivalent tȯ for i = 1, 2 . . . , n and z i (0) = e i , where e i is the unit vector with a single nonzero entry with value 1 (i.e. e 1 = 1 0 . . . 0 ⊤ ). The solution of (48) is given by Assuming that A(T − t) 2 ≤ Q for all t ∈ [0, T ], and applying the triangular inequality in (49), we obtain where the last equality comes from the fact that z i (0) 2 = e i 2 = 1 for all i = 1, 2, . . . , n. Then, applying Gronwall inequality, we have To characterize Q, note that The last inequality is obtained by applying Lemma 5 to D and noticing that each column of D is expressed as d i = e i σ ′ (Ky + b) i , where |σ ′ (x)| ≤ S for every x ∈ R by assumption. Hence, we can characterize Q as Last, having bounded the column vectors z i , (50), we apply Lemma 5 to conclude the proof.

D. Proof of Theorem 2
Define It is easy to verify thatẏ . , M such that T(i, r) = R ⊤ (r, s) = R(s, h) = 1, or equivalently (TR ⊤ R)(i, h) = 1. We conclude that the forward propagation can be implemented in a localized way according to the graph S if T(t)R(t) T R(t) ≤ S holds. Similar reasoning holds for the backward propagation. Notice that the sparsity pattern of J ⊤ (t) is the same as that of J(t) by skew-symmetricity. Define

E. Proof of Lemma 4
We study the Hamiltonian system (4) in the extended phase space [28], i.e., we define an extended state vector y = (p, q, ε, t) 12 , an extended interconnection matrixJ = where the equality ΛJΛ = ΓJΓ is verified by direct inspection. We conclude that ∂yj+1 ∂yj is a symplectic matrix. (27) We simulate the system (27) for ǫ = 1. In Figure 5 we compute the value s(T, T − τ, y 0 + γβ) − s(T, T − τ, y 0 ) γ , for γ = 0.005, β = 1 0 , 0 1 . This is obtained by selecting y(0) at random, letting y 0 = s(T − t, 0, y(0)) and then computing both s(T, T − t, y 0 + γβ) and s(T, T − t, y 0 ) by appropriate numerical integration of (27). The numerical experiment confirms that the entries of ∂y(T ) ∂y(T −t) diverge as t increases. We note that for a large enough t the values reach a maximum value, as expected because γ is not infinitesimal. One can further inspect that, when γ is chosen closer to 0, the maximum values achieved by the trajectory diverges.

APPENDIX IV IMPLEMENTATION DETAILS
DNN architectures and training algorithms are implemented using the PyTorch library 14 .

A. Binary classification datasets
For two-class classification problems we use 8000 datapoints and a mini-batch size of 125, for both training and test data. Training is performed using coordinate gradient descent, i.e. a modified version of stochastic gradient descent (SGD) with Adam (β 1 = 0.9, β 2 = 0.999) [12]. Following [12], in every iteration of the algorithm, first the optimal weights of the output layer are computed given the last updated parameters of the hidden layers, and then, a step update of the hidden layers' parameters is performed by keeping fixed the output parameters. The training consist in 50 epochs and each of them has maximum 10 iterations to compute the output layer 14 https://pytorch.org/ weights. The learning rate, or optimization step size as per γ in (21), is set to 2.5 × 10 −2 . For the regularization, α ℓ = 0, the weight decay for the output layer (α N ) is constant and set to 1 × 10 −4 and α is set to 5 × 10 −4 .