François-Xavier Vialard
email: francois-xavier.vialard@u-pem.fr

Roland Kwitt
email: roland.kwitt@sbg.ac.at

Susan Wei
email: susan.wei@unimelb.edu.au

Marc Niethammer

A Shooting Formulation of Deep Learning

A residual network may be regarded as a discretization of an ordinary differential equation (ODE) which, in the limit of time discretization, defines a continuous-depth network. Although important steps have been taken to realize the advantages of such continuous formulations, most current techniques assume identical layers. Indeed, existing works throw into relief the myriad difficulties of learning an infinite-dimensional parameter in a continuous-depth neural network. To this end, we introduce a shooting formulation which shifts the perspective from parameterizing a network layer-by-layer to parameterizing over optimal networks described only by a set of initial conditions. For scalability, we propose a novel particle-ensemble parameterization which fully specifies the optimal weight trajectory of the continuous-depth neural network. Our experiments show that our particle-ensemble shooting formulation can achieve competitive performance. Finally, though the current work is inspired by continuous-depth neural networks, the particleensemble shooting formulation also applies to discrete-time networks and may lead to a new fertile area of research in deep learning parameterization.

Introduction

Deep neural networks (DNNs) are closely related to optimal control (OC) where the sought-for control variable corresponds to the parameters of the DNN [START_REF] Liu | Deep learning theory review: An optimal control and dynamical systems perspective[END_REF][START_REF] Li | Maximum principle based algorithms for deep learning[END_REF][START_REF] Haber | Stable architectures for deep neural networks[END_REF]. To be able to talk about an optimal control requires the definition of a control cost, i.e., a norm on the control variable. We explore the ramifications of such a control cost in the context of DNN parameterization. For simplicity, we focus on continuous formulations in the spirit of neural ODEs [START_REF] Chen | Neural ordinary differential equations[END_REF]. However, both discrete and continuous OC formulations exist [START_REF] Bryson | Applied optimal control: optimization, estimation, and control[END_REF][START_REF] Athans | Optimal control: an introduction to the theory and its applications[END_REF][START_REF] Troutman | Variational calculus and optimal control: optimization with elementary convexity[END_REF]; our approach could be developed for both.

Initial work on continuous DNN formulations was motivated by the realization that a ResNet [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] He | Identity mappings in deep residual networks[END_REF] resembles Euler forward time-integration [START_REF] Haber | Stable architectures for deep neural networks[END_REF][START_REF] Li | Maximum principle based algorithms for deep learning[END_REF]. Specifically, the forward pass of some input vector x ∈ R d through a network with L layers, specified as x(0) = x and x(j + 1) = x(j) + f (x(j), θ(j)), j = 0, 1, . . . , L, closely relates to an explicit Euler [START_REF] Trefethen | Finite difference and spectral methods for ordinary and partial differential equations[END_REF] discretization of the ODE ẋ(t) = f (t, x(t), θ(t)), x(0) = x, 0 ≤ t ≤ T .

(1.1)

In the continuous DNN formulation, we seek an optimal θ such that the terminal prediction given by x(T ), i.e., the solution to Eq. (1.1) at time T , minimizes (x(T )) for a task-specific loss function .

Standard NODE

x i (0)

x i (T )

p i (0) t = 0 t = T ṗi = -∂ x f (x i (t), θ) (p i ) ẋi = f (x i (t), θ) Shooting NODE x i (0) p i (0) t = 0 t = T ṗi = -∂ x f (x i (t), θ(t)) (p i ) ẋi = f (x i (t), θ(t))
Particle shooting

x i (0)

x i (T )

t = 0 t = T ẋi = f (x i (t), θ(t)) p j (0) p i (T )
x i (T ) p i (T ) q j (0) q j (T ) p j (T ) θ = θ({(q j , p j )} K j=1 ) ṗj = -∂ q f (q j (t), θ(t)) (p j ) qj = f (q j (t), θ(t))

Figure 1: Optimization in the neural ODE (NODE) framework [START_REF] Chen | Neural ordinary differential equations[END_REF] (left) amounts to a forward pass with the gradient computed via backpropagation ( ). Optimization under the shooting principle (middle) turns the forward-backward system into a forward second-order system, where we essentially run the backpropagation equation forward. We use a Hamiltonian particle ensemble (right) consisting of K (position, momentum) pairs (q j , p j ) to make shooting efficient. Note that we write θ = θ({(q j , p j )} K j=1 ) since θ satisfies a compatibility equation which involves all K particles. In shooting θ is time-dependent, in standard NODE θ(t) = θ ∀t.

Although Eq. (1.1) with time-varying parameter θ(t) can be considered as a neural network with an infinite number of layers, current implementations of ODE-inspired networks largely assume the parameters θ are fixed in time, i.e., ∀t : θ(t) = θ [START_REF] Chen | Neural ordinary differential equations[END_REF][START_REF] Dupont | Augmented neural ODEs[END_REF], or follow some prescribed dynamics [START_REF] Zhang | ANODEV2: A coupled neural ODE framework[END_REF]. Instead, we explore time-varying θ(t) by employing regularization (i.e., a control cost) to render the estimation well-posed and to assure regularity of the resulting flow. Specifically, (for a single data point) we propose minimizing over θ the regularized loss

E(θ) = T 0 R(θ(t)) dt + γ (x(T )), γ ∈ R + , subject to Eq. (1.1) , (1.2)
where R is a real-valued complexity measure of θ corresponding to the control cost. We will mostly work with the Frobenius norm but R(θ(t)) can be more general (see Appendix B). Instead of directly optimizing over the set of time-dependent θ(t) as in standard ResNets, we restrict the optimization set to those θ which are critical points of E(θ), thereby dramatically reducing the number of parameters. In doing so, one can describe the optimization task as an initial value problem. Namely, we show that we can rewrite the loss in Eq. (1.2) solely in terms of the input x(0) and a corresponding finite-dimensional momentum variable, p(0). Such an approach, just like optimizing the initial speed of a mass particle to reach a given point, is called a shooting method in numerical analysis [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] and control [START_REF] Bonnans | The shooting approach to optimal control problems[END_REF], giving its name to our new formulation.

The first two panels of Fig. 1 illustrate the difference between the optimization of a neural ODE (NODE) via [START_REF] Chen | Neural ordinary differential equations[END_REF] and our shooting formulation. Since in practice, we have multiple inputs xi , i = 1, . . . , n, there is an initial momentum vector p i corresponding to each of them. If the shooting formulation is to scale up to a large sample size n, we must take care that the parameterization does not grow linearly with n. To this end, we propose what we call the Hamiltonian particle-ensemble parameterization. It is a finite set of particles, where each particle is a (position, momentum) pair. The initial conditions of these particle pairs {(q j , p j )} K j=1 (where K n) completely determine θ(t). This is illustrated in the rightmost panel of Fig. 1. Once the optimized set of particles has been computed, the computational efficiency of the forward model, similarly to NODE [START_REF] Chen | Neural ordinary differential equations[END_REF], is retained for vector fields f that are linear in their parameters θ(t).

Our contributions are as follows: 1) We introduce a shooting formulation for DNNs, amounting to an initial-value formulation for neural network parameterization. This allows for optimization over the original network parameter space via optimizing over the initial conditions of critical networks only; 2) We propose an efficient implementation of the shooting approach based on a novel particle-ensemble parameterization in which a set of initial particles (the (position, momentum) pairs) describe the space of putative optimal network parameters; 3) We propose the UpDown model which gives rise to explicit shooting equations; 4) We prove universality for the flows of the UpDown vector field and demonstrate in experiments its good performance on several prediction tasks.

Related work

We draw inspiration from two separate branches of research: 1) continuous formulations of neural networks [START_REF] Chen | Neural ordinary differential equations[END_REF] and 2) shooting approaches for deformable image registration [START_REF]Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation[END_REF][START_REF] Miller | Geodesic shooting for computational anatomy[END_REF][START_REF] Niethammer | Geodesic regression for image time-series[END_REF].

Continuous-depth neural networks. Continuous equivalents of ResNets [START_REF] He | Deep residual learning for image recognition[END_REF][START_REF] He | Identity mappings in deep residual networks[END_REF] have been developed in [START_REF] Ruthotto | Deep neural networks motivated by partial differential equations[END_REF][START_REF] Haber | Stable architectures for deep neural networks[END_REF], but naïve implementations are memory-demanding since backpropagation requires differentiating through the numerical integrator. Two approaches can address this unfavorable memory footprint. NODE [START_REF] Chen | Neural ordinary differential equations[END_REF] does not store intermediate values in the forward pass, but recomputes them by integrating the forward model backward. This is easily possible only if the forward model is numerically invertible and the formulation is time-continuous [START_REF] Gholami | ANODE: Unconditionally accurate memoryefficient gradients for neural odes[END_REF] 1 . Instead, checkpointing [START_REF] Gholami | ANODE: Unconditionally accurate memoryefficient gradients for neural odes[END_REF] is a general approach to reduce memory requirements by selectively recomputing parts of the forward solution [START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF]. Our work can easily be combined with these numerical approaches.

Solving implicit equations. A recent line of works, including deep equilibrium models [START_REF] Bai | Deep equilibrium models[END_REF] and implicit residual networks [START_REF] Reshniak | Robust learning with implicit residual networks[END_REF], has shown that it may not always be necessary to freely parameterize all the layers in the network. Specifically, in [START_REF] Bai | Deep equilibrium models[END_REF] and [START_REF] Reshniak | Robust learning with implicit residual networks[END_REF], the parameters of each layer are defined via an implicit equation motivated by weight tying thus improving expressiveness and reducing the number of parameters while decreasing the memory footprint via implicit differentiation. Instead, our work increases expressiveness and reduces the number of parameters via particle-based shooting.

Invertibility and expressiveness. Based on similarity with continuous time integration, constraining the norm of a layer in a ResNet will result in an invertible network such as in [START_REF] Behrmann | Invertible residual networks[END_REF][START_REF] Jacobsen | i-RevNet: Deep invertible networks[END_REF]. Invertibility is also explored in [START_REF] Younes | Diffeomorphic learning[END_REF], where it is enforced (as in our setting) via a penalty of the norm. These works show that standard learning tasks can be performed on top of a one-to-one transformation. Recent theoretical developments [START_REF] Zhang | Approximation Capabilities of Neural ODEs and Invertible Residual Networks[END_REF] show that indeed capping a NODE or i-ResNet [START_REF] Behrmann | Invertible residual networks[END_REF] with a single linear layer gives universal approximation for non-invertible continuous functions. Further, expressiveness can be increased by moving to more complex models, e.g., by introducing additional dimensions as explored in augmented NODE [START_REF] Dupont | Augmented neural ODEs[END_REF]. In [START_REF] Zhang | ANODEV2: A coupled neural ODE framework[END_REF] (AnodeV2), Zhang et al. treat time-dependent θ(t). Weights are evolved jointly with the state of the continuous DNN. While this weight evolution could, in principle, also be captured by a learned weight network, the authors argue that this would result in a large increase in parameters and therefore opt for explicitly parameterizing these evolutions (e.g., via a reaction diffusion equation). In contrast, our method does not rely on learning a separate weight-network or on explicitly specifying a weight evolution. Instead, our evolving weights are a direct consequence of the shooting equations which, in turn, are a direct consequence of penalizing network parameters (the control cost) over time; a large increase in parameters does not occur.

Hamiltonian approaches. Toth et al. [START_REF] Toth | Hamiltonian generative networks[END_REF] proposed Hamiltonian generative networks to learn the Hamiltonian governing the evolution of a physical system. Specifically, they learn Hamiltonian vector fields in the latent space of an image encoder-decoder architecture. Saemundsson et al. [START_REF] Saemundsson | Variational integrator networks for physically meaningful embeddings[END_REF] also learn the underlying dynamics of a system from time-dependent data, starting from a discrete Lagrangian combined with a variational integrator. This motivates particular network structures; e.g., Newtonian networks where the potential energy is learned via a neural network. Although sharing common tools, our work completely differs from this line of research in the sense that we exploit Hamiltonian mechanics to parameterize general continuous neural networks. In principle, our work applies to most network architectures and is not specific to physical data.

Finally, we mention that shooting approaches have been applied successfully in other areas such as diffeomorphic image matching [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF][START_REF]Diffeomorphic 3D image registration via geodesic shooting using an efficient adjoint calculation[END_REF][START_REF] Niethammer | Geodesic regression for image time-series[END_REF]. However, the decisive difference here is in the dimensionality of the underlying space: in diffeomorphic image registration, the data are points in a 3D volume i.e., d = 3; for DNNs applications, data points usually lie in a much higher-dimensional space, i.e., d is very large.

Shooting formulation of ODE-inspired neural networks

We consider, for simplicity, a supervised learning task where the input and target spaces are X ⊂ R d and Y, resp., and sampled data are denoted by {(x i , ỹi )} n i=1 ⊂ X × Y. The goal is to learn the weight θ(t) in the following flow equation

ẋi (t) = f (x i (t), θ(t)), x i (0) = xi , 0 ≤ t ≤ T, i = 1, . . . , n (3.1)
such that it minimizes the loss n i=1 (x i (T ), ỹi ) for some loss function . In existing works, the weight is chosen independent of time, i.e., θ(t) = θ [START_REF] Chen | Neural ordinary differential equations[END_REF], or specific evolution equations are postulated for it [START_REF] Massaroli | Dissecting neural ODEs[END_REF][START_REF] Zhang | ANODEV2: A coupled neural ODE framework[END_REF]. Such strategies show the difficulty of addressing infinite dimensional parameterizations of time-dependent θ and the need for regularization for wellposedness [START_REF] Finlay | How to train your neural ODE: the world of Jacobian and kinetic regularization[END_REF][START_REF] Massaroli | Dissecting neural ODEs[END_REF][START_REF] Haber | Stable architectures for deep neural networks[END_REF]. Instead of parameterizing θ(t) directly, we aim at penalizing θ(t) according to the regularity of f (•, θ(t)) to arrive at a well-posed problem. Specifically, we consider a regularization term R(θ(t)) (discussed in §3.1) and propose to minimize over θ

E n (θ) = T 0 R(θ(t)) dt + γ n i=1 (x i (T ), ỹi ), γ ∈ R + , subject to Eq. (3.1) . (3.2)
Note that upon discretizing the time t (i.e., having a number of parameters proportional to the number of timesteps) this is similar to a ResNet with weight decay. For a ResNet or a NODE, optimization is based on computing the parameter gradient via a forward pass followed by backpropagation (see left panel of Fig. 1). Optimality equations. The optimality conditions for Eq. (3.2) in continuous time are:

     ẋi (t) -f (x i (t), θ(t)) = 0, x i (0) = xi , Data evolution ṗi (t) + ∂ x f (x i (t), θ(t)) (p i ) = 0, p i (T ) = -γ∇ (x i (T ), ỹi ), Adjoint evolution ∂ θ R(θ(t)) -n i=1 ∂ θ f (x i (t), θ(t)) (p i (t)) = 0 . Compatibility (3.3)
The first equation describes evolution of the input data and the second equation is the adjoint equation solved backward in time in order to compute the gradient with respect to the parameters. At convergence, the third equation is also satisfied. This last equation encodes the optimality of the layer at timestep t, as it is the case for an argmin layer or weight tying [START_REF] Amos | OptNet: Differentiable optimization as a layer in neural networks[END_REF]. Its left hand side corresponds to the gradient with respect to the parameter θ, but as we shall see it will allow us to compute θ directly via our (position, momentum) pairs in our particle shooting formulation. The shooting approach simply replaces the optimization set by the set of critical points of Eq. (3.2) expressed in these optimality conditions. That is, we only optimize over solutions fulfilling Eq. (3.3).

Shooting principle. The shooting method is standard in optimal control [START_REF] Bonnans | The shooting approach to optimal control problems[END_REF] and can be formulated as follows: since, at optimality, the system in Eq. (3.3) is satisfied, one can turn this system into a forward model defined only by its initial conditions {(x i (0), p i (0))} n i=1 which specify the entire trajectory of optimal parameters. We evolve both the data and adjoint evolution equations forward in time and compute at each time, t, θ(t) from the compatibility Eq. (3.3) via the current values of {(x i (t), p i (t))} n i=1 . We refer to the forward model defined by Eq. (3.3) as the shooting equations. Unfortunately, this initial-condition parameterization still requires all initial conditions x i (0) and their corresponding momenta p i (0) for i = 1, . . . , n. Since this does not scale to very large datasets, we propose an approximation using a collection of particles, as described next.

Hamiltonian particle ensemble. In the limit and ideal case where the data distribution is known, the optimality equations can be approximated using a collection of particles which follow the Hamiltonian system (see Appendix A). We thus consider a collection of particles {(q j , p j )} K j=1 ∈ R d × R d that drive the evolution of the entire population {x i } n i=1 ⊂ R d through the following forward model

               ẋi (t) -f (x i (t), θ(t)) = 0, x i (0) = xi Data evolution qj (t) -f (q j (t), θ(t)) = 0, ṗj (t) + ∂ q f (q j (t), θ(t)) (p j (t)) = 0, ∂ θ R(θ(t)) - K j=1 ∂ θ f (q j (t), θ(t)) (p j (t)) = 0 ,          Hamiltonian equations (3.4)
with initial conditions {(q j (0), p j (0))} K j=1 , where the gradient with respect to this new parameterization is computed via backpropagation, and typically K n. This set of (position, momentum) pairs is termed the Hamiltonian particle ensemble. As the number of particles is reduced, so are the number of free parameters, see Appendix C. Indeed, varying the Hamiltonian particle ensemble allows for controlling the tradeoff between reconstruction and network complexity. Note that the main difference to the shooting formulation of Eq. (3.3) is that the parameterization, θ(t), is now retrieved from the shooting equations as specified by the particle collection. The original data samples, xi , are simply propagated via these parameters.

Choices of regularization, parameterization and conserved quantities

The main computational bottleneck in the forward model of Eq. (3.4) is the implicit parameterization of θ by the last equation. Making it explicit is key to render shooting computationally tractable.

Linear in parameter 2 -quadratic penalty. In the simplest case, the space of functions f is a linear space parameterized by θ(t). In this case, a quadratic penalty amounts to a kinetic penalty. Specifically, as a motivating example, consider the forward model

f (x(t), θ(t)) = A(t)σ(x(t)) + b(t), (3.5)
where σ is a component-wise activation function,

A ∈ L 2 ([0, 1], R d 2 ), b ∈ L 2 ([0, 1], R d ) and θ(t) = [A(t), b(t)]. With the quadratic regularizer R(θ(t)) = 1 2 Tr A(t) M A A(t) + 1 2 b(t) M b b(t)
, where M A , M b are positive definite matrices, the particle shooting equations are qj (t) = A(t)σ(q j (t)) + b(t), ṗj (t) = -dσ(q j (t)) A(t) p j (t),

A(t) = M A -1 (-K j=1 p j (t)σ(q j (t)) ) b(t) = M b -1 (-K j=1 p j (t)) , (3.6)
2 Obviously, an affine function of the parameters also works similarly.

with given initial conditions (p j (0), q j (0)). We emphasize that θ(t) is explicitly defined by {(p j (t), q j (t))} K j=1 and the computational cost is reduced to matrix multiplications. As is well-known [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF], the Hamiltonian flow preserves the Hamiltonian function. In the "linear in parameter -quadratic penalty" case, this preserved quantity, denoted H(p(t), q(t)) = R(θ(t)), corresponds to a (kinetic) energy of the system of particles. As a first consequence, the objective functional can be rewritten as

H(p(0), q(0))) + γ n i=1 (x i (T ), ỹi ) .
This clearly allows for direct optimization on (p(0), q(0)), i.e., shooting. As a second consequence, since the vector field has constant norm (its squared norm is the Hamiltonian), it gives a quantitative bound on the regularity of the flow map at time t = T explicit in terms of H(p(0), q(0)). In addition (Appendix A), the Rademacher complexity of the generated flows with bounded H(p(0), q(0))) can also be controlled.

Nonlinear in parameter and non-quadratic penalty. A standard ResNet structure uses vector fields of the type (in convolutional form or not)

f (x(t), θ(t)) = θ 1 (t)σ(θ 2 (t)x(t) + b 2 (t)) + b 1 (t) , (3.7)
where

θ 1 (t) ∈ L(R d , R d ) and θ 2 (t) ∈ L(R d , R d ).
We will refer to Eq. (3.7) as the singlehidden-layer vector field. This model can also be handled in our shooting approach since the shooting equations in Eq. (3.3) are completely specified by the Hamiltonian

H(p, q, θ) = R(θ) -p f (q, θ).
Automatic differentiation can be used (see Appendix D) to implement the forward model

q(t) = ∂H ∂p (p(t), q(t), θ(t)), ṗ(t) = - ∂H ∂q (p(t)
, q(t), θ(t)), θ(t) ∈ arg min H(p(t), q(t), θ(t)).

(3.8) Note that a necessary condition for solving the third equation above is in fact the compatibility equation in Eq. (3.4). Important bottlenecks appear since the third equation is nonlinear and potentially associated with a non-convex optimization problem. This could be addressed by unrolling the optimization corresponding to the last equation, resulting in increased computational cost. In addition, in this nonlinear case, the Hamiltonian function is no longer (in general) equal to R(θ(t)) even in the quadratic regularization setting. Therefore, results on the smoothness or Rademacher complexity would no longer be guaranteed as for the linear -quadratic penalty case. Last, quadratic regularization has no known theoretical results for the Rademacher complexity of functions generated by Eq. (3.7) with bounded norm. Norms for which the Rademacher complexity of this class of functions is known [START_REF] Ma | Barron spaces and the compositional function spaces for neural network models[END_REF] to be bounded are called Barron norms, which are non-smooth and non-convex, and which would add to the difficulty. To circumvent these issues while retaining expressiveness and theoretical guarantees in the linear parameterization setting, we next introduce the UpDown model.

The UpDown model

The key idea is to transform the vector field of Eq. (3.7) into a model which is linear in parameters on which the quadratic regularization can be applied. To this end, we introduce the additional state

v(t) = θ 2 (t)x(t) + b 2 (t)
which we differentiate with respect to time to obtain

v(t) = θ2 (t)x(t) + ḃ2 (t) + θ 2 (t) ẋ(t) .
Replacing ẋ(t) by its formula, we get

v(t) = θ2 (t)x(t) + ḃ2 (t) + θ 2 (t)(θ 1 (t)σ(v(t)) + b 1 (t)) .
Now overloading on notation slightly, we use the additional state variable v(t) to propose the following ODE system, denoted the UpDown model:

ẋ(t) = θ 1 (t)σ(v(t)) + b 1 (t), v(t) = θ 2 (t)x(t) + b 2 (t) + θ 3 (t)σ(v(t)) , (3.9)
with x(t) ∈ R d , v(t) ∈ R αd and introducing the (integer-valued) inflation factor α ≥ 1. For the data evolution, x i (0) are given by the data {x i }. We parameterize the v i (0) using an affine map g Θ , i.e.,

v i (0) = g Θ (x i (0)) = Θ 12 (x i (0)) + b 12 ,
where

Θ 12 ∈ L(R d , R αd ) and b 12 ∈ L(R αd ).
In Appendix E, we prove the following theorem:

Theorem 1. Given a time-dependent vector field defined on a compact domain C of R d , which is time continuous and Lipschitz, we denote by ϕ(T, x(0)) its flow at time T from starting value x(0). Then, there exists a parameterization of the UpDown model for which its solution is ε-close to the flow,

sup x(0)∈C ϕ(T, x(0)) -x(T ) ≤ ε.
Notably, in the proof, the dimension of the hidden state v is used twice: first, for having a sufficient number of neurons in Eq. (3.7) to approximate a stationary vector field (standard universality property of multilayer perceptron) and, second, for approximating time-dependent vector fields. Therefore, at the cost of introducing a possibly large number of dimensions, the UpDown model is universal in the class of time-dependent NODEs. As shown in Appendix E, this universality result transfers to our shooting formulation. Due to its additional dimensions, it is also likely to be universal in the space of functions (i.e., not necessarily injective). We focus on the UpDown model in our experiments. Note also that while we derived our theory for vector-valued evolutions for simplicity, similar linear in parameter evolution equations can for example be derived for convolutional neural networks.

Experiments

Our goal is to demonstrate that it is possible to learn DNNs by optimizing only over the initial conditions of critical networks. This is made possible via shooting and efficient via our particle parameterization. A key difference to prior work is that our approach allows to capture time-dependent (i.e., layer-dependent in the discrete setting) parameters without discretizing these parameters at every time-point. Comparisons to other NODE like methods are not straightforward due to hyper-parameters and different implementations. For consistency, we therefore provide four different formulations (based on the UpDown model of §3.2).

• The static direct model forgoes the Hamiltonian particle ensemble, and instead directly optimizes over time-constant parameters: θ(t) = θ for all t. Everything else, including the UpDown model, stays unchanged. This model is most closely related to NODE [START_REF] Chen | Neural ordinary differential equations[END_REF] and augmented NODE [START_REF] Dupont | Augmented neural ODEs[END_REF].

• We call our proposed shooting model dynamic with particles. It is parameterized via a set of initial conditions of (position, momentum) pairs, which evolve over time and fully specify θ(t). • The static with particles model is similar to the static direct model. However, instead of directly optimizing over a time-constant θ, it uses a set of (position, momentum) pairs (i.e., particles, as in our dynamic with particles model above) to parameterize θ indirectly.

• Finally, we consider the dynamic direct model which uses a piece-wise time-constant θ(t). It essentially chains together multiple static direct models and is closely related to a discrete ResNet in the sense that multiple blocks (we use five) are used in succession. However, each block involves time-integrating the UpDown model. While the dynamic with particles model captures θ(t) indirectly via particles and shooting, the dynamic direct model requires many more parameters as it represents θ(t) directly. We show results for the dynamic direct model for a subset of the experiments.

All experiments use the UpDown model with quadratic penalty function R. Detailed experimental settings, including weights for the quadratic penalty function, can be found in Appendix F.

Simple 1D function regression. We approximate a simple quadratic-like function y = x 2 + 3/(1 + x 2 ) which is non-invertible. We use 15 particles for our experiments. Fig. 2 shows the test loss and the network complexity, as measured by the log Frobenius norm integrated over time [START_REF] Neyshabur | Exploring generalization in deep learning[END_REF], for the different models as a function of the inflation factor α (cf. §3.2). On average, the dynamic with particles model shows the best fits with the lowest complexity measures, indicating the simplest network parameterization. Note that the static with particles approach results in the lowest complexity measures only because it cannot properly fit the function as indicated by the high test loss. Additional results for a cubic function y = x 3 are in Appendix G.

Spiral. Next, we revisit the spiral ODE example of [START_REF] Chen | Neural ordinary differential equations[END_REF] following the nonlinear dynamics ẋ = Ax 3 , x ∈ R 2 (where the power is component-wise). We fix x(0) = [2, 0] T , use A = [-0.1, 2.0; -2, -0.1] and evolve the dynamics for time T = 10. The training data consists of snippets from this trajectory, all of the same length. We use an L 2 norm loss (calculated on all intermediate time-points) and 25 particles. Our goal is to show that we can obtain the best fit to the training data due to our dynamic model. Fig. 3 (top) shows that we can indeed obtain similar or better fits (lower losses) for a similar number of parameters while achieving the lowest network complexity measures. Fig. 3 (bottom) shows the corresponding results for the validation data consisting of the original long trajectory starting from initial value x(0). Interestingly, by pasting together short-range solutions we are successful in predicting the long-range trajectory despite training on short-range trajectory snippets.

Concentric circles. To study the impact of the inflation factor α in a classification regime, we replicate the concentric circles setting of [START_REF] Dupont | Augmented neural ODEs[END_REF].

The task is learning to separate points, sampled from two disjoint annuli in R 2 . While we are less interested in the learned flow (as in [START_REF] Dupont | Augmented neural ODEs[END_REF]), we study how often the proposed of α and 20 particles. Notably, the effect of α is only visible if the classification loss is down-weighted so that the regularization, R, dominates. Otherwise, for the tested α, the model always fits the data. The experiment is consistent with [START_REF] Dupont | Augmented neural ODEs[END_REF], where it is shown that increasing the space on which an ODE is solved allows for easy separation of the data and leads to less complex flows. The latter is also observed for our model. Rotating MNIST. Here, we are given sequences of a rotating MNIST digit (along 16 angles, linearly spaced in [0, 2π]). The task is learning to synthesize the digit at any rotation angle, given only the first image of a sequence. We replicate the setup of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF] and consider rotated versions of the digit "3". We identify each rotation angle as a time point t i and randomly drop four time points of each sequence during training. One fixed time point is consistently left-out and later evaluated during testing. We use the same convolutional autoencoder of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF] with the UpDown model operating in the internal representation space after the encoder.

α =128 α =64 α =32 α =16 α =128 α =64 α =32 α =16 log 2 (

complexity) Loss

Short-range results

Long

α =128 α =64 α =32 α =16 α =128 α =64 α =32 α =16

Number of parameters

Static direct

Static with particles

Dynamic with particles

During training, the encoder receives the first image of a sequence (always at angle 0 • ), the UpDown model integrates forward to the desired time points, and the decoder decodes these representations. As loss, we measure the mean-squared-error (MSE) of the decoder outputs. Fig. 4 lists the MSE (at the left-out angle), averaged over all testing sequences and shows two example sequences with predictions for all time points (100 particles, α = 10).

While all UpDown variants substantially lower the MSE previously reported in the literature, they exhibit comparable performance. To better understand the differences, we visualize the internal representation space of the autoencoder by projecting all 16 internal representations (i.e., the output of the UpDown models after receiving the output of the encoder) of each testing image onto the two largest principal components, shown to the right (different colors indicate the different rotation angles). This qualitative result shows that allowing for a time-dependent parameterization leads to a more structured latent space of the autoencoder.

Bouncing balls. Finally, we replicate the "bouncing balls" experiment of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF]. This is similar to the rotating MNIST experiment, but the underlying dynamics are more complex.

In particular, we are given 10,000 (training) image sequences of bouncing balls at 20 different time points [START_REF] Sutskever | The recurrent temporal restricted Boltzmann machine[END_REF]. The task is learning to predict, after seeing the first three images of a sequence, future time points. We use the same convolutional autoencoder of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF] and minimize image (per-pixel) MSE (using all 20 time points for training). Our UpDown model operates in the internal representation space of the encoder (50-dimensional in our experiments3 ). In test mode, the network receives the first three image of a sequence and predicts 10 time points ahead. We measure the image (per-pixel) MSE and average the results (per time point) over all 500 testing sequences. For model selection, we rely on the provided validation set. Our UpDown (dynamic with particles) model uses 100 particles. Fig. 5 (left) lists the averaged MSE per time point, plotted against the approaches listed in [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF]. Fig. 5 (right) shows two testing sequences with predictions (the three input time points are not shown). Results for the UpDown static and static with particles model are 0.0154 and 0.0150, respectively. Computational cost. The computational cost of the UpDown model consists in storing the particles and running forward the model for the collection of particles and the data. Hence, computational cost scales linearly in the number of particles. To get rid of this linear relationship (in case only a forward pass is needed), the ODE can be discretized in time and the ResNet with its weights is obtained.

Discussion and Conclusions

We demonstrated that it is possible to parameterize DNNs via initial conditions of (position, momentum) pairs. While our experiments are admittedly still simple, results are encouraging as they show that 1) the particle-based approach can achieve competitive performance over direct parameterizations and that 2) time-dependent parameterizations are useful for obtaining simpler networks and can be realized with significantly fewer parameters using particle-based shooting.

Our work opens up many different follow-up questions and formulations. For example, we presented our approach for a model with continuous dynamics, but the particle and the shooting formalism can also be applied to discrete-time models. Further, we focused, for simplicity, on continuous variants of multi-layer perceptrons, but similar linear-in-parameter models can be formulated for convolutional neural networks. Models that are nonlinear in their parameters hold the promise for connections with optimal mass transport theory and to theoretical complexity results, which we touched upon for our UpDown model. Indeed, this change of paradigm in the parameterization may result in new quantitative results on network generalization properties. Lastly, how well the approach generalizes to more complex problems, how many particles are needed to switch from a standard deep network to its shooting formulation, and how optimizing over critical points of the original optimization problem via shooting relates to network generalization will be fascinating to explore.

Source code is available at: https://github.com/uncbiag/neuro_shooting

Broader Impact

One goal of this work is to enrich the understanding of continuous depth neural networks and to open a different (or alternative) perspective on its parameterization. Specifically, we shift the parameterization of deep neural networks from a layer-by-layer perspective to an initial-value perspective and Hamiltonian dynamics. At this point, our work is conceptual and theoretical in nature; broader impact emerges most likely as a consequence of better understanding the role of neural network parameterizations.

where γ is a positive regularization parameter. We now consider the particular case of a ResNet model where each layer is given by an UpDown model Eq. (3.9) or even a single hidden layer Eq. (3.7).

Without loss of generality, set the terminal time to T = 1. Letting ρ 0 denote the probability density of X, minimizing Eq. (A.2) is equivalent to minimizing

inf ϕ Reg(ϕ) + γ R d (ϕ(x))ρ 0 (x) dx .
This can be rewritten as

inf ϕ Reg(ϕ) + γ R d (x )ρ 1 (x ) dx ,
where ρ 1 (x) := ρ(1, x) is the flow of the continuity equation

∂ t ρ(t, x) + div(ρ(t, x)f (x, θ)) = 0 , ρ(0, x) = ρ 0 (x) ,
where div is the divergence operator on vector fields. Note that ρ 1 can be regarded as the density representing the data at time 1. In the following, we deal with a general regularization term Reg(ϕ) = 1 0 R(θ(t), ρ(t)) dt, where the R term can depend on the density of data at time t. A particular though important case is when the regularization R does not depend on ρ,

1 0 R(θ(t)) dt.

A.2 Optimality equations and Hamiltonian ensemble approximation

We detail the optimality equations when data points are represented by a probability measure. As mentioned above, the regularity of the map is enforced via a penalty on the weights at each timepoint and is the integral 1 0 R(θ(t)) dt or even more generally 1 0 R(θ(t), ρ(t)) dt. Using Lagrange multipliers, this constraint can be enforced and minimizers of the energy should be saddlepoints of the energy

L(ρ, θ, p) := γ R d (x)ρ 1 (x) dx + 1 0 R(θ(t), ρ(t)) dt + 1 0 R d p(t, x)(∂ t ρ(t, x) + div(ρ(t, x)f (x, θ(t)))) dx dt ,
where p(t, x) is a time and space dependent function. The optimality equations are then

     ∂ t ρ(t, x) + div(ρ(t, x)f (t, x, θ(t))) = 0 , ∂ t p(t, x) + ∇p(t, x) • f (t, x, θ(t)) = δR δρ (θ(t), ρ(t)) , ∂ θ R(θ(t), ρ(t)) - R d ∂ θ f (x, θ(t)) (∇p(t, x)ρ(t, x)) = 0 , (A.3)
where ∇p is the gradient w.r.t. x of p(t, x) and δ denotes differentiation w.r.t. the indicated parameter. The notation δR δρ means the Fréchet derivative of the penalty w.r.t. the density ρ. Note that in our current work, R is independent of ρ. However, this more general setup encompasses optimal transport models, see Section B.2.

In practice, one does not have access to the full distribution and the variational setup needs to be approximated. As proposed in the main text, we approximate it using a collection of particles that follow the optimality equations which are Hamiltonian evolution equations for this collection of particles. The collection of particles {(q j , p j )} are defined by their state and costate. We estimate ρ using the empirical measure 1 K K j=1 δ q j (t) (•). Writing the optimality equation for this particular empirical measure leads to the equation Eq. (3.4). When the number of particles tends to infinity, we can hope to recover the optimal trajectory. However, we do not explore this question formally here. We simply remark that this question is directly connected to expressiveness and generalization properties of the constructed neural network and is also probably data dependent.

B Choice of regularization

The simplest regularization on the flow ϕ is given by

Reg(ϕ) = 1 0 R(θ(t)) dt , (B.1)
where R does not depend on ρ(t, •). The first possibility is a quadratic penalty for the single-hidden-layer vector field of Eq. (3.7), where R(θ(t)) = 1 2 θ(t) 2 2 is the Frobenius norm of the parameter θ. Since the space of vector fields is a finite dimensional linear space, it can be endowed with a scalar product, which turns this space into a Reproducing Kernel Hilbert Space (RKHS). Therefore, the linear in parameter -quadratic penalty setting of §3.1 is a particular case of vector fields encoded by f (•, θ(t)) ∈ H, with H a RKHS embedded in W 1,∞ vector fields. This setting leverages strong analytical and geometrical foundations [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Bruveris | On completeness of groups of diffeomorphisms[END_REF]:

1) When the activation function is smooth, the resulting vector field is smooth 4 , and consequently the associated flow map ϕ is guaranteed to be a one-to-one smooth map (i.e., a diffeomorphism). For instance, with the UpDown model, it is a homeomorphism in (x, v). Moreover, the quadratic penalty induces a right-invariant distance on the set of flows generated by Eq. (A.1) and the distance to identity of the resulting flow can be bounded by Reg(ϕ) (see [START_REF] Younes | Shapes and Diffeomorphisms[END_REF][START_REF] Bruveris | On completeness of groups of diffeomorphisms[END_REF] for more details in a Sobolev setting). 2) When the activation function is of ReLU type, the resulting map is still a W 1,∞ one-to-one map (i.e., a homeomorphism) and has Lipschitz regularity.

Another type of regularization for the single-hidden-layer vector field of Eq. (3.7) we discuss is based on the Barron norm [START_REF] Ma | Barron spaces and the compositional function spaces for neural network models[END_REF]:

θ 2 B := 1 d d j=1 θ j 1 2 2 ( [θ 2 ] j 1 + b j 2 1 ) 2 ,
where θ j 1 denotes the j th column of θ 1 and [θ 2 ] j denotes the j th row of θ 2 . As discussed in the main text, the reason we might consider a Barron norm penalty for the single-hidden-layer vector field in Eq. (3.7) rather than the quadratic penalty is because of its theoretical results. Indeed, the Rademacher complexity is bounded for the combination of a single-hidden-layer vector field with a Barron norm penalty, but not when combined with a quadratic penalty.

B.1 Linear in parameters -quadratic energy

Now let us examine in detail models that are linear in parameters and have quadratic energy on parameters: this case is the simplest to be studied, and computationally not as demanding as the nonlinear case. As mentioned above, the set of possible vector fields f (•, θ(t)) is a finite dimensional linear space, which is a reproducing kernel Hilbert space when endowed with an L 2 norm. Since all Hilbert norms in finite dimensions are equivalent, this choice of regularization is universal in this class of quadratic penalties.

1. The vector field is f (•, θ(t)) = θ • σ, where σ is a vector of maps. In this case, the optimality equation reads

∂ θ f (x, θ(t)) (∇p(t, x)ρ(t, x)) = R d σ(x) (∇p(t, x)ρ(t, x)) dx .
2. If the penalty R only depends on θ and is quadratic: R(θ(t)) = 1 2 1 0 θ(t) 2 dt, then one has δR δθ (θ(t), ρ(t)) = θ(t). Thus, under these two conditions, the parameters are explicit in terms of p, ρ and σ:

θ(t) = R d σ(x) (∇p(t, x)ρ(t, x)) dx . (B.2)
Two observations are warranted. First, if, instead of quadratic regularization on the parameters, we were to choose a RKHS norm (in the infinite dimensional case) as penalty, it would result in the introduction of the kernel applied to the R.H.S. of Eq. (B.2). Second, from Eq. (B.2), one could be tempted to derive an evolution equation for θ. This equation is known as the EPDiff equation [START_REF] Younes | Shapes and Diffeomorphisms[END_REF] and is unfortunately not a closed equation on the set of parameters θ(t) themselves. Therefore, our approach is a possible way to approximate it.

An important property of this simple setting is that the norm of the vector field is preserved by the forward model defined by the collection of Hamiltonian particles and it also holds in the continuous setting. As stated in Section 3.1, the Hamiltonian is given by

R(θ(t)) = 1 2 Tr A(t) M A A(t) + 1 2 b(t) M b b(t)
where A, b are the optimal parameters given by

A(t) = M A -1 (-K j=1 p j (t)σ(q j (t)) ) b(t) = M b -1 (-K j=1 p j (t)) . (B.3)
The Hamiltonian R(θ(t)) being constant gives a constant norm vector field.

B.2 Nonlinear in parameters -energy which depends on the distribution

For exposition purposes, we present two cases of interest which we have not well explored numerically.

Example of the Barron norm. Obviously, the single-hidden-layer vector field in Eq. (3.7) is not linear in parameters. We have already discussed that it is proper in this case to endow the space with norms such as the Barron norm [START_REF] Ma | Barron spaces and the compositional function spaces for neural network models[END_REF]. For simplicity, consider the single-hidden-layer vector field in Eq. (3.7) without b 1 , i.e., f (x(t), θ(t)) = θ 1 σ(θ 2 (x) + b 2 ). A simple upper bound for the Barron norm 5 is

f (•, θ) 2 B := 1 d d j=1 θ j 1 2 2 ( [θ 2 ] j 1 + [b 2 ] j 1 ) 2 . (B.4)
Again, θ j 1 denotes the j th column of θ 1 and [θ 2 ] j denotes the j th row of θ . Let us consider the case of R(θ

(t)) = 1 2 f (•, θ) 2 B .
In this case, one has the following optimality equations to solve

θ j 1 ( [θ 2 ] j 1 + |[b 2 ] j 1 ) 2 = R d σ([θ 2 ] j x + [b 2 ] j ) (∇p(t, x)ρ(t, x)) dx , θ j 1 2 2 ( [θ 2 ] j 1 + [b 2 ] j 1 )∂ [θ 2 ] k j 1 = R d [ dσ([θ 2 ] j x + [b 2 ] j )(x k )] (∇p(t, x)ρ(t, x)) dx , θ j 1 2 2 ( [θ 2 ] j 1 + [b 2 ] j 1 )∂ [b 2 ] j k 1 = R d [ dσ([θ 2 ] j x + [b 2 ] j )(x k )] (∇p(t, x)ρ(t, x)) dx .
These equations involve the subdifferential of the L 1 norm, and optimization of this type of functions, which involves sparsity, is a well-explored field [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF]. We leave experiments with this norm for future work. Note that in this case the norm of the vector field is not equal to the Hamiltonian and it is not a constant of the flow.

B.3 L 2 regularization, optimal transport

Last, we briefly mention a model that is part of our framework which has the advantage of not specifying the penalty on the space of parameters encoding the vector field. In case there is no obvious norm to be used on the space of vector fields, it is possible to use an L 2 type of penalty on the vector fields themselves instead of on the parameters.

Indeed, one way to be rather independent of the choice of the parameterization of the map consists in introducing a cost that represents the L 2 norm of the map. However, L 2 depends on the choice of a measure and this measure can be chosen as the density of the data, ρ(t, x). More precisely, one can use

R(f, ρ(t)) = 1 2 R d f (x, θ) 2 ρ(t, x) dx . (B.5)
In such a case, this formulation resembles finding an optimal transport (OT) map between ρ 0 and ρ 1 . Specifically, optimal transport is an optimization problem which can be solved via a fluid dynamic formulation [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] introducing the kinetic penalty above. However, the two models (OT and the one defined by the regularization Eq. (B.5)) differ since the optimization set for optimal transport is the set of L 2 vector fields with respect to measure ρ and the above formulation is a parameterized approximation of this set. This parameterized approximation needs to retain generalization properties of the optimized map. Note however, that in the limit where the number of neurons goes to infinity, optimal transport will be well-approximated since the optimization is performed on a dense subset of all vector fields. Obviously, fixing the choice to a single-hidden-layer design implies a choice for [START_REF] Bai | Deep equilibrium models[END_REF], which thus gives a regularization of the computed approximation of the optimal transport map.

d in θ 1 (t) ∈ L(R d , R d ) and θ 2 (t) ∈ L(R d , R d ) of Eq. (3.
Computational burden. In either case of the Barron norm or the optimal transport type of penalty, the implicit equation corresponding to the third equation in Eq. (A.3) has to be solved at each layer of the discretization. We experimented with a simple strategy of unrolling the related minimization scheme. An efficient approach to solve such implicit equations will be necessary for practical implementations.

B.4 Rademacher complexity of bounded energy flows.

In this section, given a set of vector fields with bounded Rademacher complexity, we show that the resulting flows also have bounded Rademacher complexity. The flow of a vector field f (•, θ(t)) is a vector valued map denoted by ϕ. Let us first treat the case of the Rademacher complexity of a component of the flow map ϕ k . Theorem 2. Let F be a space of vector fields defined on a compact space C ⊂ R d . Assume that the Rademacher complexity on n points in C of each component of the vector fields f k (t, •) for k = 1, . . . , d is controlled by M (n, t) which depends on n, then the Rademacher complexity of each component of the flows at time 1 is bounded by

1 0 M (n, t) dt.
Proof. Recall that Rademacher complexity, see [START_REF] Wainwright | High-Dimensional Statistics: A Non-asymptotic Viewpoint[END_REF], of a class of functions F is defined as, for

Z = (z 1 , . . . , z n ) ∈ C, Rad Z (F) def. = E sup g∈F n i=1 ε i g(z i ) ,
where the {ε i } n i=1 are i.i.d. Rademacher random variables. Our hypothesis ensures Rad Z (F) ≤ M (n). Apply the definition of the flow to get

ϕ(1, x) = x + 1 0 f (ϕ(t, x), θ(t)) dt .
Therefore,

E sup ϕ∈F n i=1 ε i ϕ k (z i ) ≤ Rad Z ({Id}) + 1 0 E sup f (•,θ(t)) n i=1 ε i f k (ϕ(t, z i ), θ(t)) dt , ≤ 0 + 1 0 M (n, t) dt .
In the previous formula, we used the fact that the Rademacher complexity of a set comprised of a single map is zero. Corollary 3. Let H be a RKHS of vector fields whose kernel k is bounded on the diagonal k(x, x) k ∞ < ∞ , then, the set of flows denoted by F at time 1 of time-dependent vector fields in B(0, R), the ball of radius R centered at the origin satisfies Rad Z (F)

≤ 2R √ k(x,x) k ∞ √ n
, where Rad Z (F) is the Rademacher complexity for n points.

Proof. The Rademacher complexity of the ball of radius R in the RKHS H [8, Lemma 22] is

upper bounded: Rad Z (B(0, R)) ≤ 2R √ k(x,x) k ∞ √ n
. We then directly apply Theorem 2.

A similar result also holds for vector fields generated by the single-hidden-layer vector field in Eq. (3.7), see [START_REF] Ma | Barron spaces and the compositional function spaces for neural network models[END_REF]. Last, we note that the result and its proof also hold if one uses the following Rademacher complexity for vector valued functions [START_REF] Maurer | A vector-contraction inequality for rademacher complexities[END_REF],

Rad Z (F) def. = E   sup g∈F n i=1 d j=1 ε j ε i g j (z i )   , for g = (g j ) j=1,...,d ∈ F.

B.5 Consequences for the UpDown model

We put together the previous results on the UpDown model. First, the space of vector fields endowed with the quadratic penalty on the parameters forms a RKHS. The variational formulation implies that the norm of the velocity field generated by a given collection of Hamiltonian particles {(q j (0), p j (0))} is preserved. In addition, this norm can be explicitly computed since the parameters at time 0 can be computed in terms of {(q j (0), p j (0))}. Last, the generated space of maps has a Rademacher complexity which is linearly bounded by this norm. In order to be fully explicit on the constant for the Rademacher complexity, we need to compute sup x∈C k(x, x) where k is the kernel associated with the RKHS. Without making this quantity explicit here, we simply mention that the bound degrades (i.e. increases) with increasing inflation factor α, as it can be expected.

C Analysis of the number of free parameters

It is instructive to understand the number of parameters for a shooting approach in comparison to the typical approach of optimizing a neural network (where the parameter-dependency at optimality is only considered implicitly at convergence of the numerical solution rather than explicitly during the shooting). We focus on the cases of affine and convolutional layers for illustration.

Consider a DNN with a depth of L layers, where each hidden layer has P parameters. The number of free parameters is then LP , compared to 2KS where K is the number of active particles, each of them of size 2S 6 . Hence, solutions with less than LP/(2S) particles provide benefits in the number of free parameters. Therefore, as the number of particles is reduced, we may parameterize the DNN with a smaller number of parameters. Most remarkably, the number of free parameters is always 2KS regardless of the number of parameters of a particular layer as the layer parameters are obtained via the shooting equations based on the particle states. This is a consequence of regularizing the parameters in our loss which couples them across time at optimality. We make this clearer in what follows. 6 For example, S for our UpDown model simply corresponds to the dimension of its state space: S = (α + 1)d, where α is the inflation factor and d the data dimension. Note that in our experiments with the UpDown model we also learned an affine map from the initial conditions x(0) to the initial conditions v(0). Such a map has αd(d + 1) parameters. These parameters are included in the table of Fig. 3 and in Tables 1/2 summarizing the number of mode parameters. However, we will not consider parameters in our discussion here, as they would equally apply to both a shooting and a direct optimization approach and could also be avoided by simply initializing v(0) to zero. A similar initialization to zero approach is, for example, commonly taken in ResNets when increasing the number of feature channels [START_REF] He | Deep residual learning for image recognition[END_REF].

Affine layers. Recall that in our simple example of §3.1 the parameters θ(t) = [A(t), b(t)] of our affine7 model are given as

A(t) = M A -1   - K j=1 p j (t)σ(q j (t))   , b(t) = M b -1   - K j=1 p j (t)   . (C.1)
Here, A(t) and b(t) have d 2 and d parameters, respectively; these parameters are indirectly given by the set of particles {(q j (t), p j (t))} at any given time. Hence, for this model S = d and P = d(d+1). If we assume we have K particles and compare to a discrete layer implementation of this model then the particle-based approach will have less free parameters if

2Kd < Ld(d + 1) .
Importantly, the state-space dimension, d, only enters the number of free parameters linearly for the particle approach (2Kd), while there is a quadratic dependence for direct optimization (Ld(d + 1)). This is a direct consequence of the optimality condition which couples the parameters θ(t) across time. One can see this phenomenon in action in Eq. (C.1), where the matrix A is expressed as the sum of matrices p j (t)σ(q j (t)) with rank ≤ 1. Concretely, a particle-based shooting approach uses less parameters if the number of particles K < L(d+1)/2. Another interesting observation based on this example is that even if we would have only considered a linear model (i.e., without the bias term, b(t)) the number of parameters for the particles would have still remained at 2KS. This is again a consequence of optimality and of our parameter regularization. Note that this also means that even though our UpDown model

ẋ(t) = θ 1 (t)σ(v(t)) + b 1 (t), v(t) = θ 2 (t)x(t) + b 2 (t) + θ 3 (t)σ(v(t)) , has significantly more parameters θ(t) = [θ 1 (t), b 1 (t), θ 2 (t), b 2 (t), θ 3 (t)
] when directly optimized, this has no direct impact on the number of free parameters of its particle-based parameterization. Only the state-space dimension matters. Concretely, if we were to instead consider a model of the form

ẋ(t) = θ 1 (t)σ(v(t)), v(t) = θ 2 (t)x(t) ,
the particle-based parameterization would stay unchanged! Only the way how one infers θ(t) from the particles changes. Convolutional layers. Shooting approaches for convolutional models can also be derived. We did not experiment with such models in this work. However, we show here that the number of free parameters may also be decreased with a particle-based approach. This will be interesting to explore in future work. Specifically, for convolutional layers a particle-based parameterization could be particularly effective as one typically has quadratic complexity in the number of filters between convolutional layers (i.e., if a layer with N feature channels is followed by a layer with M feature channels, this will induce the estimation of N × M convolutional filters and hence will drastically influence the number of parameters for large N or M ). In contrast, a particle-based shooting approach does not increase the number of parameters as it ties them together via the optimality conditions expressed by the shooting equations. As a rough estimate for a standard convolutional ResNet for L = 50, P = 100 2 ×16, LP ≈ 8.10 6 . Thus, if particles have size 40, we end up with at most 10 5 active particles.

General remarks. Nevertheless, all model parameters (e.g., [A(t), b(t)] or all convolutional filters for a convolutional layer) are still instantiated during computation. It is important to note that regardless of the chosen number of particles, a shooting neural network solution is a possible optimal solution (for a given data set) at any given time, not only at convergence. One optimizes over the family of possible neural network models with the goal of finding the element within this family that best matches the observations. The system can be rewritten as

q(t) = θ 1 (t)σ(v(t)) + b 1 (t) , v(t) = θ 3 (t)q(t) + θ 4 (t)σ(v(t)) + b 3 (t) . (E.1)
Therefore, with the initial condition v(0) = θ 2 (0)q(0) and q(0) = q 0 , the two systems of ordinary differential equations are equivalent.

Note that the key point in Lemma 4 is the loss of regularity in the evolution of θ 2 since we differentiated once in time. For that reason, we now show that adding more dimensions using the inflation factor α alleviate this issue. It is likely possible that one could prove a universality result using only α = 1 but we shall leave this question for future work 8 . However, experimentally, the inflation factor has a crucial effect on the performance of the optimization, as discussed in §4. Lemma 4 helps us establish the next result.

Lemma 5. Consider the single-hidden-layer vector field in Eq. (3.7) with θ(t) being piecewise continuous. Then, there exists a parameterization of the UpDown model that gives the same flow.

Proof. Without loss of generality, we only treat the case of one discontinuity in time of the parameterization; We thus assume that θ(t) is continuous on [0, t 1 [ and [t 1 , 1]. We consider q, v 1 , v 2 ∈ R d such that q, v 1 are defined as in Lemma 4. We now define, up to time t 1 , v 2 (t) = θ(t 1 )v 1 (t) + θ 2 (t 1 ) which implies (differentiating w.r.t. time) that v 2 follows an evolution equation similar to v 1 and thus can be encoded in the general form of Eq. (E.1). Now, q(t), v 2 (t) are defined on [t 1 , 1] by the evolution Eq. (E.1) in order to coincide with the flow of single-hidden-layer vector field on [t 1 , 1], q(t) = θ 1 (t)σ(v 2 (t)) + b 1 (t) and v2 (t) = θ 3 (t)q(t) + θ 4 (t)σ(v 2 (t)) + b 3 (t) for well chosen parameters as in Lemma 4. Since the value of v 1 (t) is not used in the evolution equation of q(t), we can simply extend it by v 1 (t) = v 1 (t 1 ) which is a valid evolution equation for Eq. (E.1).

In the general case, we decompose the time interval [0, 1] into k intervals [t i , t i+1 [ on which θ(t) is continuous and the proposed method can be directly extended using an inflation factor

α = k, introducing v k ∈ R d .
Note that the result of this lemma gives an equality between the two flows defined on the whole space R d . The next result is an approximation result which holds on a compact domain C ⊂ R d . For a function f : R d → R, we denote f C,∞ = sup x∈C |f (x)|. Proposition 6. The UpDown model is universal in the class of time-dependent vector fields. Let C ⊂ R d be a compact domain. For every time-dependent vector field (such that it is time continuous and is Lipschitz in space) w : [0, 1] × R d → R d and its associated flow ϕ(t, x) there exist time dependent parameters of the UpDown model such that

q(t) = θ 1 (t)σ(v(t)) + b 1 (t) , v(t) = θ 2 (t)(q(t)) + b 2 + θ 3 (t)σ(v(t)) , is ε-close to the solution ϕ(1, x), e.g. ϕ(1, x) -q(1, x) C,∞ ≤ ε.
Proof. The proof is standard and we include it here for self-containedness. It is the consequence of [START_REF]Approximation by superpositions of a sigmoidal function[END_REF] and Lemma 5. Let B(0, r) a ball of radius r in R d which contains ϕ(t, x) for all time t ∈ [0, 1]. The flow associated with a given time-dependent vector field v(t, •) can be approximated by a vector field which is piecewise constant in time; i.e. let ε > 0 be a positive real, (by continuity in time of v(t, •)) there exists a decomposition of [0, 1] into k intervals [t i , t i+1 ] and Lipschitz vector fields v i (x) = f (x, θ i ) where f is the single hidden layer in Eq. (3.7) such that v i (x) -v(t, x) B(0,r),∞ ≤ ε for t ∈ [t i , t i+1 ]. Denote by w(t, •) the time-dependent vector field defined by w(t, •) = v i (•) for all t ∈ [t i , t i+1 ]. Thus, denoting the flow of v(t, •) by ϕ v and the flow of w(t, •) by ϕ w , we get

ϕ v (1, x) -ϕ w (1, x) ≤ 1 0 v(t, ϕ v (t, x)) -v(t, ϕ w (t, x)) + v(t, ϕ w (t, x)) -w(t, ϕ w (t, x)) dt ϕ v (1, x) -ϕ w (1, x) C,∞ ≤ 1 0 Lip(v) ϕ v (t, x) -ϕ w (t, x) C,∞ + v(t, •) -w(•) B(0,r),∞ dt ≤ 1 0 Lip(v) ϕ v (t, x) -ϕ w (t, x) C,∞ dt + ε ,
where Lip(v) denotes a bound on the Lipschitz constant of v(t, x) w.r.t. x ∈ B(0, r) for all t ∈ [0, 1]. Then, the Grönwall lemma [START_REF] Younes | Shapes and Diffeomorphisms[END_REF] gives

ϕ v (1, x) -ϕ w (1, x) C,∞ ≤ εe Lip(v) . (E.2)
By Lemma 5, ϕ w (1, x) can be approximated by the flow of the UpDown and the result is obtained via the triangle inequality.

In this section, we focused on a universality result in the space of time-dependent vector fields. Interestingly, due to the additional dimensions, it is likely that the model is universal in the space of functions as well. This conjecture is supported by the quadratic 1D function regression example which shows that the UpDown model is able to capture some maps which are not homeomorphic. We leave this question for future work.

F Experimental settings

This section describes our experimental settings. We use our UpDown model for all experiments and simply use a weighted Frobenius norm penalty for all parameters. Specifically, we weigh this penalty for all parameters with 1 except for, θ 3 which we penalize by 10. In our experiments, we have observed better convergence properties for higher penalties on θ 3 . This might be due to the special role that θ 3 plays in the model as it subsumes a quadratic term in the original derivation of the UpDown model (see §3.2). In all experiments, we also optimize over the affine map from x(0) to v(0) for the data evolution.

Simple function regression. We use 500 epochs for all experiments. For all particlebased experiments we freeze the positions of the particles for the first 50 epochs. We use a ReLU activation function and the MSE loss. We weigh the MSE loss by 100 and the parameter norm loss by 1. We use 500 training samples, 1,000 testing samples and 1,000 validation samples and a batch size of 50. Note that for these simple examples there is, in practice, no real difference between the training, testing, and validation data, as the number of samples is large and the domain is [-1.5, 1.5]. We initialize the particle positions uniformly at random in [-1.5, 1.5] and draw the momenta from a Gaussian distribution with zero mean and standard deviation 0.1. All time-integrations are done via a fourth-order Runge-Kutta integrator with time-step 0.1. For optimization, we use Adam with a learning rate of 0.01 and the ReduceLROnPlateau learning rate scheduler of PyTorch with a learning rate reduction factor of 0.5.

Spiral. The spiral data is generated between time t = 0 and t = 10 with 200 uniformly spaced timepoints. Training is only on small time snippets with an approximate length of 0.25 time-units. Evaluation is on these short time snippets as well as on the entire trajectory by pasting together solutions for these short time snippets, i.e., an individual short solution starts where the previous one ends. Settings for the spiral are the same as for the simple function regression with the following exceptions. We use 1,500 epochs and the step-size for the fourth-order Runge-Kutta integrator is 0.05. The MSE loss is still weighted by 100, but the parameter norm loss only by 0.01. We randomly draw 100 new training samples during each epoch and use 100 evaluation samples and 1,000 short range samples and 1 long-range testing sample. All samples are randomly drawn from the trajectory. However, as the trajectory is traversed at highly nonuniform speed the samples are drawn from a uniform distribution across the trace of the spiral. As for the simple function regression experiment, there is little practical difference between the training, validation, and testing data as the problem is so simple. However, this is not of concern in these experiments as the prime objective is to study the fitting behavior of the different models. We use a batch size of 100.

Rotating MNIST. We use the data provided by the authors of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF] and follow the same autoencoder architecture, except that our encoder maps into a 20-dimensional representation space. The number of particles is set to 100 and the inflation factor α is set to 10. For optimization, we use Adam with a learning rate of 0.001 and the CosineAnnealingLR learning rate scheduler of PyTorch. We train for 500 epochs with a batch size of 25 and the parameter norm loss set to 0.1.

Bouncing balls. As in the rotating MNIST experiment, we rely on the data provided by the authors of [START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF], follow their autoencoder architecture and set the dimensionality of the representation space of the encoder to 50. The first three images of each sequence are provided to the encoder by concatenating the images along the channel dimension. The inflation factor α is set to 20 and we use 100 particles. We optimize over 100 epoch using Adam with the CosineAnnealingLR learning rate scheduler of PyTorch, the initial learning rate is set to 0.001 and the parameter norm loss is set to 0.0001.

G Additional results

Simple function regression. In Section 4, we considered approximating a quadratic-like function. Here we show parallel results for approximating a cubic function y = x 3 . We will also include some additional figures for the quadratic-like regression function. Note that whereas the cubic function is invertible (but not diffeomorphic), the quadratic-like one considered in Section 4 is a simple example of a non-invertible function. Tab. 1 shows the number of parameters for the four different formulations for both regression functions. Fig. 6 shows for the cubic regression the test loss and the network complexity, as measured by the Frobenius norm [START_REF] Neyshabur | Exploring generalization in deep learning[END_REF], for the four formulations. On average the particle-based approaches show the best fits with the lowest complexity measures, indicating the simplest network parameterization. Note however that while the dynamic particle approach greatly outperformed the static particle approach for the quadratic-like function (see Fig. 2) this is not the case here. In fact, the static particle approach shows slightly better fits than the dynamic one. This might be because the cubic function is significantly simpler to fit and hence may not benefit as much from the dynamic approach. To illustrate that fitting the quadratic-like function is indeed harder, Figs. 7 and8 show function fits for different numbers of particles for the cubic function and the quadratic-like function, respectively. All these fits are for the particle-based dynamic UpDown model. Clearly, very few particles can achieve reasonable fits for a simple function. As little as two particles already show a good fit for the cubic function, whereas the quadratic-like function requires with more particles. This supports our hypothesis that fitting more complex functions may require more particles.

Since our approach is based on the time-integration of the UpDown model it is interesting to see 1) how the mapping is expressed across time and 2) how the parameters, θ(t), of the UpDown model change over time. Fig. 9 shows example mappings for the cubic and the quadratic-like function, respectively. The estimated mappings are highly regular. Lastly, Figs. 10 and 11 show the time-evolutions of the model parameters for the cubic and the quadratic-like function for two different inflation factors. While different parameters show different dynamics, clear changes over time can be observed. In particular, θ 2 (t) and b 2 (t) show strong changes. These parameters mostly control the behavior of the hidden high-dimensional state, v, as θ 3 (t) is penalized significantly more in our model (see Sec.F) and consequently shows more moderate changes. Spiral. Tab. 2 shows the number of parameters in each of the four formulations for the spiral experiment. This table complements the Table in Fig. 3 which only showed the number of parameters when using 15 particles. As can be seen, the mappings are highly regular.

Figure 2 :

 2 Figure 2: Fit for quadratic-like y = x 2 + 3/(1 + x 2 ) for 10 random initializations. Left: Test loss; Right: time-integral of log 2 of the Frobenius norm complexity. Lower is better for both measures. * indicates number of removed outliers (outside the interquartile range (IQR) by ≥ 1.5× IQR); α denotes the inflation factor.

FitFigure 3 :

 3 Figure3: Fit for spiral (short-and long-range). Losses for the different models as well as the time-integral of log 2 of the Frobenius norm complexity measure. Lower is better for both measures. The * symbol indicates how many outliers were removed and α denotes the inflation factor.

Figure 4 :

 4 Figure 4: Left: Image (per-pixel) MSE (measured at the marked time point) averaged over all testing sequences of the rotated MNIST dataset. Right: Two testing sequences and predictions (marked blue)for all 16 time points when the image at t = 0 is given as input (marked red). Results marked with † are taken from[START_REF] Yildiz | ODE 2 VAE: Deep generative second order ODEs with Bayesian neural networks[END_REF].

Figure 5 :

 5 Figure 5: Left: Image (per-pixel) MSE for predicting 10 time points ahead (after receiving the first three inputs of a sequence), averaged over all testing sequences (numbers in parentheses indicate the MSE when additionally averaged over all prediction time points). Results marked with † are taken from [40]. Right: Two testing sequences with predictions (marked blue).

Figure 6 :

 6 Figure 6: Function fit (15 particles) for cubic y = x 3 for 10 random initializations. Left: Test loss; Right: time-integral of log 2 of the Frobenius norm complexity. Lower is better for both measures. * indicates number of removed outliers (outside the interquartile range (IQR) by ≥ 1.5× IQR); α denotes the inflation factor.

Figure 7 :

 7 Figure7: Fits for the cubic function with inflation factor 16 and for different numbers of particles. Vertical lines indicate particle positions after optimization. While subtle, the figures suggest that using more particles allows for better approximation of the function. This is confirmed by the test loss values in Fig.6(bottom left).

Figure 8 :

 8 Figure 8: Fits for the quadratic-like function for inflation factor 16 with different of particles.Vertical lines indicate particle positions after optimization. As this function is more complex than the cubic function 2 and 5 particles is not sufficient for a fit. But 15 and 25 particles result in a well-fitting approximation.

Figure 9 :

 9 Figure 9: Mapping of the cubic function (left) and the quadratic-like function (right). As can be seen, the mappings are highly regular.

  -range results

	α
	16	1,282	1,116	6,026
	32	4,610	2,172	22,282
	64	17,410	4,284	85,514
	128	67,586	8,508	334,858

  MSE ± σ † GPPVAE-DIS 0.0309 ± 0.00002 † GPPVAE-JOINT 0.0288 ± 0.00005 † ODE 2 VAE 0.0194 ± 0.00006 † ODE 2 VAE-KL 0.0184 ± 0.0003 Ours (stat. direct) 0.0126 ± 0.0064 Ours (stat. w particles) 0.0125 ± 0.0063 Ours (dyn. w particles) 0.0122 ± 0.0064 ↔

	Test time point

Table 1 :

 1 Number of parameters for the simple function regression cubic and quadratic experiments.

	Inflation factor

Table 2 :

 2 Number of parameters for the spiral experiment.

				Inflation factor
		#Particles	16	32	64	128
		15	1,116 2,172	4,284	8,508
	static/dynamic w/ particles	25	1,796 3,492	6,884	13,668
		50	3,496 6,792 13,384 26,568
	static direct	n/a	1,282 4,610 17,410 67,586
	dynamic direct	n/a	6,026 22,282 85,514 334,858

In a discrete setting, resolving the forward model in the backward direction generally requires costly solving of implicit equations. This can be done (it is, e.g., done for invertible ResNets[START_REF] Behrmann | Invertible residual networks[END_REF]). In general, an explicit numerical solution for forward time-integration becomes implicit in the backward direction and vice versa.

We did not further experiment with this hyperparameter, so potentially better results can be obtained.

I.e smoothness asks for Lipschitz regularity vector field, which ensures existence and uniqueness of the flow.

The actual Barron norm is defined as the infimum of the r.h.s. in Eq. (B.4) on all the possible representations of the function f (•, θ) as a single-hidden-layer.

In this section, we mean affine with respect to (σ(qj(t)))j∈1,...,K .

Note that the case α = 1 is similar in its formulation to a second-order model on q.

Acknowledgments and Disclosure of Funding

This research project was initiated during a one-month invitation of M. Niethammer by the Labex Bézout, supported by the French National Research Agency ANR-10-LABX-58. R. Kwitt is partially funded by the Austrian Science Fund (FWF): project FWF P31799-N38 and the Land Salzburg (WISS 2025) under project numbers 20102-F1901166-KZP and 20204-WISS/225/197-2019. S. Wei is the recipient of an Australian Research Council Discovery Early Career Award (project number DE200101253) funded by the Australian Government.

Supplementary material

The following sections discuss in more detail the theoretical guarantees of our approach. §A presents the optimality conditions underlying our shooting formulation and it is shown how these optimality equations can be approximated via a collection of particles. §B proposes different regularizations, whose choice is key for practical and theoretical results. We show that, under some conditions, the Rademacher complexity of the set of flows can be bounded and apply our results in §B.5 to the UpDown model. §C discusses the number of free parameters of our shooting approach in relation to the number of free parameters for direct optimization. §D explains how the shooting equations can be automatically derived via automatic differentiation.

§E shows the universality of our UpDown model. §F provides details on our experimental setup. Lastly, §G shows some additional experimental results.

A Expectation approximation of optimality equations

We first discuss a general variational setup of supervised learning including regularization.

A.1 Variational setup

Suppose the data consists of input X ∈ R d . Let f (•, θ(t)) be a vector field on R d , e.g. the single hidden layer of Eq. (3.7) or a linear (in parameter) layer. Consider the flow ϕ := ϕ(T, •) generated by f according to

We consider the general task of minimizing,

D Automatic shooting

The general shooting equations were presented in Eq. (3.3). We then proceeded to explicitly derive the shooting equations for a continuous DNN with linear-in-parameter layers and UpDown layers in §3.1 and §3.2, respectively. While this was instructive, it is somewhat cumbersome, in particular, for more complex models or when moving to convolutional networks. Fortunately, in practice these shooting equations do not need to be derived by hand. Indeed, they are completely specified by the Hamiltonian

in the sense that the shooting equations in Eq. ( 3.3) are computed via differentiation of H. Specifically, the shooting equations in Eq. ( 3.3) are equivalently given by

As discussed above, the last equation can be replaced by solving

Automatic differentiation can be used to automatically obtain the shooting equations. As fitting a shooting model requires differentiating the shooting equations, we in effect end up with differentiating twice. This can be done seamlessly using modern deep learning libraries, such as PyTorch.

E Universality of the UpDown model

In this section, we set out to demonstrate that the UpDown model is universal in the sense that its associated flow can come ε-close to the flow of any well behaving time-dependent vector field.

Recall the single-hidden-layer vector field in Eq. (3.7) with time-varying parameters θ(t) = (θ 1 (t), θ 2 (t), b 1 (t), b 2 (t)). While shooting with the single hidden layer vector field is theoretically appealing as it is universal [START_REF]Approximation by superpositions of a sigmoidal function[END_REF], it would result in implicit shooting equations. We first show that the UpDown model introduced in 3.2 can give the same flow as the single hidden layer (Lemma 5) and then leverage this relationship to show that the UpDown model inherits the universality of the single hidden layer (Proposition 6). Lemma 4. Consider the single-hidden-layer vector field in Eq. (3.7) with θ 2 (t) and b 2 (t) being piecewise C 1 and θ 1 (t), b 1 (t) continuous. Then, there exists a parameterization of the UpDown model that gives the same flow at a fixed time, T = 1.

Proof. We rewrite the differential equation

by introducing the additional state variable v(t) = θ 2 (t)q(t) + b 2 (t) which we differentiate w.r.t. time. We obtain v(t) = θ2 (t)q(t) + ḃ2 (t) + θ 2 (t) q(t) . Replacing q(t) by its formula, we get v(t) = θ2 (t)q(t) + ḃ2 (t) + θ 2 (t)θ 1 (t)(σ(v(t)) + b 1 (t)) .

Figure 10: Weight evolution across time (i.e., continuous depth) for 15 particles when fitting the cubic function using the UpDown model:

Results are for the dynamic with particles approach. Top: Inflation factor 16. Bottom: Inflation factor 64. Changes in parameter values can clearly be observed.