NeuPDE: Neural Network Based Ordinary and Partial Differential Equations for Modeling Time-Dependent Data

2019·arXiv

Abstract

Abstract

We propose a neural network based approach for extracting models from dynamic data using ordinary and partial differential equations. In particular, given a time-series or spatio-temporal dataset, we seek to identify an accurate governing system which respects the intrinsic differential structure. The unknown governing model is parameterized by using both (shallow) multilayer perceptrons and nonlinear differential terms, in order to incorporate relevant correlations between spatio-temporal samples. We demonstrate the approach on several examples where the data is sampled from various dynamical systems and give a comparison to recurrent networks and other data-discovery methods. In addition, we show that for MNIST and Fashion MNIST, our approach lowers the parameter cost as compared to other deep neural networks.

Keywords— partial differential equations, data-driven models, data-discovery, multilayer perceptrons, image classification

1 Introduction

Modeling and extracting governing equations from complex time-series can provide useful information for analyzing data. An accurate governing system could be used for making data-driven predictions, extracting large-scale patterns, and uncovering hidden structures in the data. In this work, we present an approach for modeling time-dependent data using differential equations which are parameterized by shallow neural networks, but retain their intrinsic (continuous) differential structure.

For time-series data, recurrent neural networks (RNN) is often employed for encoding temporal data and forecasting future states. Part of the success of RNN are due to the internal memory architecture which allows these networks to better incorporate state information over the length of a given sequence. Although widely successful for language modeling, translation, and speech recognition, their use in high-fidelity scientific computing applications is limited. One can observe that a sequence generated by an RNN may not preserve temporal regularity of the underlying signals (see, for example [5] or Figure 2.3) and thus may not represent the true continuous dynamics.

For imaging tasks, deep neural networks (DNN) such as ResNet [11, 12], FractalNet [19], and DenseNet [14] have been successful in extracting complex hierarchical spatial information. These networks utilize intra-layer connectivity to preserve feature information over the network depth. For example, the ResNet architecture uses convolutional layers and skip connections. The hidden layers take the form represents the features at layer n and F is a convolutional neural network (or more generally, any universal approximator) with trainable parameters . The evolution of the features over the network depth is equivalent to applying the forward Euler method to the ordinary differential equation (ODE): ˙). The connection between ResNet’s architecture, numerical integrators for differential equations, and optimal control has been presented in [7, 10, 23, 30, 40].

Recently, DNN-based approaches related to differential equations have been proposed for data mining, forecasting, and approximation. Examples of algorithms which use DNN for learning ODE and PDE include: learning from data using a PDE-based network [21, 22], deep learning for advection equations [6], approximating dynamics using ResNet with recurrent layers [25], and learning and modeling solutions of PDE using networks [28]. Other approaches for learning governing systems and dynamics involve sparse regularizers (or hard-thresholding approaches in [4, 29, 33, 34] and problems in [32, 35, 39]) or models based on Gaussian processes [26, 27].

Note that in [21, 22] it was shown that adding more blocks of the PDE-based network improved (experimentally) the model’s predictive capabilities. Using ODEs to represent the network connectivity, [5] proposed a ‘continuous-depth’ neural network called ODE-Net. Their approach essentially replaces the layers in ResNet-like architectures with a trainable ODE. In [5], the authors state that their approach has several advantages, including the ability to better connect ‘layers’ due to the continuity of the model and a lower memory cost when training the parameters using the adjoint method. The adjoint method proposed in [5] may not be stable for a general problem. In [8], a memory efficient and stable approach for training a neural ODE was given.

1.1 Contributions of this Work.

We present a machine learning approach for constructing approximations to governing equations of time-dependent systems that blends physics-informed candidate functions with neural networks. In particular, we construct a network approximation to an ODE which takes into account the connectivity between components (using a dictionary of monomials) and the differential structure of spatial terms (using finite difference kernels). If the user has prior knowledge on the structure or source of the data, i.e. fluids, mechanics, etc., one can incorporate commonly used physical models into the dictionary. We show that our approach can be used to extract ODE or PDE models from time-dependent data, improve the spatial accuracy of reduced order models, and reduce the parameter cost for image classification (for the MNIST and Fashion MNIST datasets).

2 Modeling Via Ordinary Diﬀerential Equations

Given discrete-time measurements generated from an unknown dynamic process, we model the time-series using a (first-order) ordinary differential equation, ˙1. The problem is to construct an approximation to the unknown generating function f, i.e. we will learn networks net(t, x) such that ˙Essentially, we are learning a neural network approximation to the velocity field. Following the approach in [5], the system is trained by a ‘continuous’ model and the function f is parameterized by multilayer perceptrons (MLP). Since a two-layer MLP may require a large width to approximate a generic (nonlinear) function f, we purpose a different parameterization. Specifically, to better capture higher-order correlations between components of the data and to lower the number of parameters needed in the MLP (see for example, Figure 2.2), a dictionary of candidate inputs is added. Let D(t, x; p) be the collection (as a matrix) of the pth order monomial terms depending on t and x, i.e. each element in D can be written as:

One has freedom to determine the particular dictionary elements; however, the choice of monomial terms provides a model for the interactions between each of the components of the time-series and is used for model identification of dynamical systems in the general setting [4, 35, 39]. For simplicity, we will suppress the (user-defined) parameter p.

In [4, 29, 35, 39], regularized optimization with polynomial dictionaries is used to approximate the generating function of some unknown dynamic process. When the dictionary is large enough so that the ‘true’ function is in the span of the candidate space, the solutions produced by sparse optimization are guaranteed to be accurate. To avoid defining a sufficiently rich dictionary, we propose using an MLP (with a non-polynomial activation function) in combination with the monomial dictionary, so that general functions may be well-approximated by the network. Note that the idea of using products of the inputs appears in other network architectures for example, the high-order neural networks [9, 38].

In many DNN architectures, batch normalization [15] or weight normalization [31] are used to improve the performance and stability during training. For the training of NeuPDE, a simple (uniform) normalization layer, N(x), is added between the input and dictionary layers, which maps x to a vector in [(using the range over the all components). Specifically, let M and m be the maximum (and minimum) value of the data, over all components and samples and define the vector N(x) as:

This normalization is applied to each component uniformly and enforces that each component of the dictionary is bounded by one (in magnitude). We found that this normalization was sufficient for stabilizing training and speeding up optimization in the regression examples. Without this step, divergence in the training phase was observed.

To train the network: let be the vector of learnable parameters in the MLP layer, then the optimization problem is:

where 0 are regularization parameters set by the user and F is an MLP. Specifically, let be a smooth activation function, for example, the exponential linear unit (ELU)

or the hyperbolic tangent, tanh, which will be sufficiently smooth for integration using Runge-Kutta schemes. The right-hand side of the ODE is parameterized by a fully connected layer - activation layer - fully connected layer, i.e.

where ), i.e. the vectorization of all components of the matrices and biases . Therefore, the first layer of the MLP in the form linear combination of candidate functions (applied to normalized data). Note that the dictionary does not include the constant term since we include a bias in the first fully connected layer. The function r is a regularizer on the parameters (for example, the norm) and the time-derivative is penalized by the norm. When used, the parameters are set to tuning is performed).

The constraints in Eqn. (2.1) are written in continuous-time, i.e. the value of x(t) is defined by the ODE and thus can be evaluated at any time ]. For a given set of parameters the values ) are obtained by numerical integration (for example, using a Runge-Kutta scheme). To optimize Eqn. (2.1) using a gradient-based method, the back-propagation algorithm or the adjoint method (following [5]) can be used. The adjoint method requires solving the ODE (and its adjoint) backward-in-time, which can lead to numerical instabilities. Following the approach in [8], checkpointing can be used to mitigate this issue.

For all experiments, we take the ‘discretize-then-optimize’ approach. The objective function, Eqn. (2.1), is discretized as follows:

where Φis an ODE solver (i.e. a Runge-Kutta scheme) applied rescaled by the time-step, and the time-derivative is discretized on the time-interval with the integral approximated by piece-wise constant quadrature. The constraint that the ODE ˙) is satisfied at every time-stamp has been transformed to the constraint that the sequence is generated by the forward evolution of an ODE solver. The ODE solver takes (as its inputs) the initial data ) and the function F that defines the RHS of the ODE. Note that the ODE solver can be ‘sub-grid’ in the sense that, over a time interval [], we can set the solver to take multiple (smaller) time-steps. This will increase storage cost needed for back-propagation; however, taking multiple time-steps can better resolve complex dynamics embedded by examples below). Additionally, the time-derivative regularizer helps to control the growth of the generative model, which yields control over the solution x(t) and its regularity.

Eqn. (2.2) is solved using the Adam algorithm [16] with temporal mini-batches. Each mini-batch is constructed by taking a random sampling of the time-stamps then collecting all data-points from 0. This can be easily extended to multiple time-series, by sampling over each trajectory as well. Note that, in experiments, taking non-overlapping sub-intervals did not change the results. The back-propagation algorithm [24] applied to each of the subintervals [] can be done in parallel. For all of the regression

Figure 2.1: Time-series data generated by the 3d Lorenz system and the corresponding learned processes using our approach. The original data (dashed) and the learned series (solid) are plotted, where red, green, and blue curves correspond to the components, respectively.

examples, we set is the given (sequential) data. For the image classification examples, we used the standard cross-entropy for L.

Remark 2.1. Layers: The right-hand side of the ODE is parameterized by one set of parameters Therefore, in terms of DNN layers, we are technically only training one “layer”. However, changes in the structure between time-steps are taken into account by the time-dependence, i.e. the dictionary terms D that contain t. Thus, we are embedding multiple-layers into the time-dependent dictionary.

’: Eqn. (2.2) is fully discrete when using a Runge-Kutta solver for Φ and its gradient can be calculated using the back-propagation algorithm. If we used an adaptive ODE solver, such as Runge-Kutta 45, the forward propagation would generate a new set of time-stamps (which always contain the time-stamps ) in order to approximate the forward evolution ˙), given an error tolerance and a set of parameters tested the continuous-depth versions, using back-propagation and the adjoint method (see Appendix A). Although the backward evolution of the ODE constraint may not be well-posed (i.e. numerical unstable or unbounded), our experiments lead to similar results. This could be due to the temporal mini-batches which enforce short-integration windows. Following [8], a checkpointing approach should be used, which would also control the memory footprint. It should be noted that, the adjoint approach was tested for the ODE examples but not for PDE examples, since the PDE examples are not time-reversible and will lead to backward-integration issues.

In addition, when the network is discrete, one may still consider it as a ‘continuous-depth’ approximation, since the underlying approximation can be refined by adding more time-stamps, without the need to retrain the MLP.

2.1 Autonomous ODE.

When the time-series data is known to be autonomous, i.e. the ODE takes the form ˙x = f(x), one can drop the t-dependency in the dictionary. In this case, the monomials take the form We train this model by minimizing:

where ˜is the given data (corrupted by noise) over the time-stamps

which emits chaotic trajectories.

In Figure 2.1(a), we train the model with 20 hidden nodes per layer using a quadratic dictionary, i.e. there are 9 terms in the dictionary, with 20 bias parameters, with 3 bias parameters, for a total of 263 trainable parameters. The solid curves are the time-series generated by a forward pass of the trained model. The learned system generates a high-fidelity trajectory for the first part of the time-interval. In Figure 2.1(b-c), we investigate the effect of the degree in the dictionary. In Figure 2.1(b), using a degree 1 monomial dictionary with 38 hidden nodes per layers, i.e. 3 terms in the dictionary, with 38 bias parameters, bias parameters (for a total of 269 trainable parameters), the generated curves trace a similar part of phase space, but are point-wise inaccurate. By increasing the hidden nodes to 100 per layer (3 terms in the dictionary, with 100 bias parameters, with 3 bias parameter, for a total of 703 trainable parameters), we see in Figure 2.1(c) that the method (using a degree 1 dictionary) is able to capture the correct point-wise information (on the same order of accuracy as Figure 2.1(a)) but requires more than double the number of parameters.

Figure 2.2: Extracting and modeling of a nonlinear time-dependent spiral. Noisy data is given in green and learned model is given in blue.

2.2 Non-Autonomous ODE and Noise.

corrupted by noise. The third coordinate of Eqn. (2.5) is time-dependent, which can be challenging for many recovery algorithms. This is partly due to the redundancy introduced into the dictionary by the time-dependent terms. To generate Figure 2.2, we set:

(a) the dictionary degree to 1, with 4 terms, with 46 bias parameters, with 3 bias parameters (371 trainable parameters in total);

(b) the degree to 4, with 69 terms, with 4 bias parameter, bias parameters (295 trainable parameters in total), and the regularization parameter to ˜

(c) the degree to 4, with 69 terms, with 5 bias parameter, with 3 bias parameters (368 trainable parameters in total).

For cases (b-c), we set the degree of the dictionary to be larger than the known degree of the governing ODE in order to verify that we do not overfit using a higher-order dictionary and that we are not tailoring the dictionary to the problem. In Figure 2.2(a), the dictionary of linear monomials with a moderately sized MLP seems to be insufficient for capturing the true nonlinear dynamics. This can be observed by the over-smoothing caused by the linear-like dynamics. In Figure 2.2(c), a nonlinear dictionary can fit the data and extract the correct pattern (the ‘squared’-off corners). Figure 2.2(b) shows that we are able to decrease the total number of parameters and fit the trajectory within the same tolerance as (c) by penalizing the derivative. Both (b) and (c) have achieved a mean-squared loss under 0.015.

2.3 Comparison for Extracting Governing Models.

Comparison with SINDy. We compare the results of Figure 2.2 with an approximation using the SINDy algorithm from [4] (theoretical results of convergence and relationship to the appear in [41]). These approaches differ, since the SINDy algorithm seeks to recover a sparse approximation to the governing system given one tuning parameter and is restricted to the span of the dictionary elements. To make the dictionary sufficiently rich, the degree is set to 4 as was done for Figure 2.2 (b-c). Since the sparsity of the first two components is equal to one, we search over all parameter-space (up to 6 decimals) that yields the smallest non-zero sparsity. The smallest

which is independent of x and y and does not lead to an accurate approximation to the nonlinear spiral. This is likely due to the level of noise present in the data and the incorporation of the time-component.

Comparison with LASSO-based methods. We compare the results of Figure 2.2 with LASSO-based approximations for learning governing equations [32]. The LASSO parameter is chosen so that the sparsity of the solution matches the sparsity of the true dynamics (with respect to a dictionary of degree 4). In addition, the coefficients are ‘debiased’ following the approach in [32]. The learned system is:

which matches the profile of the data in the (x, y)-plane; however, it does not predict the correct dynamics for z (emitting seemingly periodic orbits). While the LASSO-based approach better resolves the state-space dependence, it does not correctly identify the time-component.

Comparison with RNN. In Figure 2.3 the Lorenz system (see Figure 2.1) is approximated by our proposed approach and a standard LSTM (RNN), with the same number of parameters. Although the RNN learns internal hidden states, the RNN does not learn the correct regularity of the trajectories thus leading to sharp corners. It is worth noting that, in experiments, as the number of parameters increases, both the RNN and our network will produce sequences that approach the true time-series.

Figure 2.3: Comparing dynamics between the true (green) time-series, the solution generated by our network (red), and a solution generated by an RNN (blue) with the same number of parameters.

2.4 ODE from Low-Rank Approximations

For certain spatio-temporal systems, reduced-order models can be used to transform complex dynamics into low-dimensional time-series (with stationary spatial modes). One of the popular methods for extracting the spatial modes and identifying the corresponding temporal-dynamics is the dynamic mode decomposition (DMD) introduced in [36]. The projected DMD method [37] makes use of the SVD approximation to construct the modes and the linear dynamical system. Another reduced-order model, known as the proper orthogonal decomposition (POD) [13], can be used to construct spatial modes which best represent a given spatio-temporal dataset. The projected DMD and the POD methods leverage low-rank approximations to reduce the dimension of the system and to construct a linear approximation to the dynamics (related to the spectral analysis of the Koopman operator), see [17] and the citations within.

We apply our approach to construct a neural network approximation to the time-series generated by a low-rank approximation of the von K´arm´an vortex sheet. We explain the construction for this example here but for more details, see [17]. Given a collection of measurements where (are the spatial coordinates and are the time-stamps, define X as the matrix whose columns are the vectorization of each where m is the number of grid points used to discretize Ω. The SVD of the data is given by are unitary matrices and Σ is a diagonal matrix. The best r-rank approximation of X is given by restriction of Σ to the top r singular values and are the corresponding singular vectors. The columns of the matrix represent the r spatial modes that can be used as a low-dimensional representation of the data. In particular, we define the vector projection of the data (i.e. the columns of X) onto the span of

Thus, we can construct the time-stamps ˜) from the measurements X and can train the system using a version Eqn. (2.1) with the constraint that the ODE is of the form:

The additional matrix resembles the standard linear structure from the DMD approximation and the function f can be seen as a nonlinear closure for the linear dynamics. The function f is approximated, as before, by ). To train the model, we minimize:

where also includes the trainable parameters from Note that, to recover an approximation to the original measurements ), the vector is mapped back to the correct spatial ordering (inverting the vectorization process).

In Figure 2.4, our approach with an 8 mode decomposition is compared to an 8 mode DMD approximation. The DMD approximation in Figure 2.4(a) introduces two erroneous vortices near the bottom boundary. Our approach matches the test data with higher accuracy, specifically, the relative error between our generated solution at the terminal time is 0.049 compared to DMD’s relative error of 0.060. It is worth noting that this example shows the benefit of the additional term ) in the low-mode limit; however, using more modes, the DMD method becomes a very accurate approximation. Unlike the standard DMD method, our model does not require the data to be uniformly spaced in time.

Figure 2.4: Learning reduced-order dynamics from fluid simulations.

3 Partial Diﬀerential Equations

A general form for a first-order in time, a-th order in space, nonlinear PDE is:

where denotes the collection of all i-th order spatial partial derivatives of We form a dictionary ]) as done in Sec. 2, where the monomial terms now apply to . The spatial derivatives as well as calculated numerically from data using finite differences. We then use an MLP, F, to parametrize the governing equation:

see also [29, 32]. In particular, the function F can be written as:

where are collections of 11 convolutions, are biases, are all the parameters from is ELU activation function. The input channels are the monomials determined by is extended to a constant 2d array. The first linear layer maps the dictionary terms to multiple hidden channels, each defined by their own 1 1 convolution. Thus, each hidden channel is a linear combination of input layers. Then we apply the ELU activation, followed by a 1 1 convolution, which is equivalent to taking linear combinations of the activated hidden channels. Note that this differs from [29, 32] in several ways. In the first linear layer, our network uses multiple linear combinations rather than the single combination as in [29, 32]. Additionally, by using a (nonlinear) MLP we can approximate a general function on the coordinates and derivative; however, previous work defined approximations that model functions within the span of the dictionary elements.

To illustrate this approach, we apply the method to two examples: a regression problem using data from a 2d Burgers’ simulation (with noise) and the image classification problem using the MNIST and MNIST-Fashion datasets.

3.1 Burgers’ Equation

We consider the 2d Burgers’ equation,

The training and test data are generated on (, with time-step ∆t = 132 uniform grid. To make the problem challenging, the training data is generated using a sine function in x as the initial condition, while the test data uses a sine function in y as the initial condition. We generate 5 training trajectories by adding noise to the initial condition. Our training set is of size [5, 100, 32, 32] and our test data is of size [1, 100, 32, 32]. To train the parameters we minimize:

where Ωis a discretization of Ω.

Training, Mini-batching, and Results. The mini-batches used during training are constructed with mini-batches in time and the full-batch in space. For our experiment, we set a (temporal) batch size of 16 with a length of 3, i.e. each batch is of size [16, 3, 32, 32] containing 16 short trajectories. The points are chosen at random, without overlapping. The initial points of each mini-batch are treated as the initial conditions for the batch, and our predictions are performed over the length of the trajectory. This is done at each iteration of the Adam optimizer with a learning rate of 0.1.

In Figure 3.1, we take 2000 iterations of training, and evaluate our results on both the training and test sets. Each of the 1 1 convolutional layers have 50 hidden units, for a total of 2301 learnable parameters. For visualization, we plot the learned solution at the terminal time on both the training and test set. The mean-squared error on the full training set is 0.005 and on the test set is 3.6 (for reference, the mean-squared value of the test data is over 1000).

3.2 Image Classification: MNIST Data.

Another application of our approach is in reducing the number of parameters in convolutional networks for image classification. We consider a linear (spatially-invariant) dictionary for Eqn. (3.1). In particular, the right-hand side of the PDE is in the form of normalization, ReLU activation, two convolutions, and then a final normalization step. Each convolutional layer uses a 3 of the form with 6 trainable parameters, where 3 kernels that represent the identity and the five finite difference approximations to the partial derivatives and . In CNN [11, 12, 20], the early features (relative to the network depth) typically appear

Figure 3.1: Burgers’ Equation Example. The surfaces at the terminal time simulated by the NeuPDE, with (a) training data and (b) the learned surface on the test data.

to be images that have been filtered by edge detectors [30, 40]. The early/mid-level trained kernels often represent edge and shape filters, which are connected to second-order spatial derivatives. This motivates us to replace the 3 3 convolutions in ODE-Net [5] by finite difference kernels.

Result shows that even though the trainable set of parameters are decreased by a third (each kernel has 6 trainable parameters, rather than 9), the overall accuracy is preserved (see Table 3.1). We follow the same experimental setup as in [5], except that the convolutional layers are replaced by finite differences. We first downsample the input data using a downsampling block with 3 convolutional layers. Specifically, we take a 3 3 convolutional layer with 64 output channels and then apply two 3 3 convolutional layers with 64 output channels and a stride of 2. Between each convolutional layer, we apply batch-normalization and a ReLU activation function. The output of the downsampling block is of size 8 8 with 64 channels. We then construct our PDE block using 6 ‘PDE’ layers, taking the form:

We call each subinterval (indexed by i) a PDE layer since it is the evolution of a semi-discete approximation of a coupled system of PDE (the particular form of the convolutions act as approximations to differential operators). The function G takes the form:

where BN(x) is batch-normalization, is a collection of 3 3 kernels of the form contains all the learnable parameters, and ) is the ReLU activation function. The PDE block is followed by batch-normalization, the ReLU activation function, and a 2d pooling layer. Lastly, a 64 10 fully connected layer is used to transform the terminal state (activated and averaged) of the PDE blocks to a 10 component vector.

For the optimization, the cross-entropy loss is used to compare the predicted outputs and the true label. We use the SGD optimizer with momentum set to 0.9. There are 160 total training epochs; we set the learning rate to 0.1 and decrease it by 1/10 after epoch 60, 100 and 140. The training is stopped after 160 epochs. All of the convolutions performed after the downsampling block are linear combinations of the 6 finite difference operators rather than the traditional 3 convolution.

Table 3.1: Comparison Between Networks on MNIST

3.3 Image Classification: Fashion MNIST.

For a second test, we apply our network to the Fashion MNIST dataset. The Fashion MNIST dataset has 50000 training figures and 10000 test figures that can be classified into 10 different classes. Each data is of size 32 32. For our experiment, we did not use any data augmentation, except for normalization before training. The structure of NeuPDE in this example differs from the MNIST experiment as follows: (1) we downsample the data once instead of twice and (2) after the 1st PDE block, which has 64 hidden units, we use a 3 3 kernel with 64 input channels and 128 output channel, and then use one more PDE block with 128 hidden units followed by a 128 fully connected layer. There are numerous reported test results on the Fashion MNIST dataset [1]; we compare our result to ResNet18 [2] and a simple MLP [3]. The test results are presented in Table 3.2 and all algorithms are tested without the use of data augmentation. In both the MNIST and Fashion MNIST experiments, we show that the NeuPDE approach allows for less parameters by enforcing (continuous) structure through the network itself.

Table 3.2: Comparison on Fashion MNIST

4 Discussion

We propose a method for learning approximations to nonlinear dynamical systems (ODE and PDE) using DNN. The network we use has an architecture similar to ResNet and ODENet, in the sense that it approximates the forward integration of a first-order in time differential equation. However, we replace the forcing function (i.e. the layers) by an MLP with higher-order correlations between the spatio-temporal coordinates, the states, and derivatives of the states. In terms of convolution neural networks, this is equivalent to enforcing that the kernels approximate differential operators (up to some degree). This was shown to produce more accurate approximations to complex time-series and spatio-temporal dynamics. As an additional application, we showed that when applied to image classification problems, our approach reduced the number of parameters needed while maintaining the accuracy. In scientific applications, there is more emphasis on accuracy and models that can incorporate physical structures. We plan to continue to investigate this approach for physical systems.

In imaging, one should consider the computational cost for training the networks versus the number of parameters used. While we argue that our architecture and structural conditions could lead to models with fewer parameter, it could be potentially slower in terms of training (due to the trainable nonlinear layer defined by the dictionary). Additionally, we leave the scalability of our approach for larger imaging data set, such as ImageNet, to future work. For larger classification problems, we suspect that higher-order derivatives (beyond second-order) may be needed. Also, while higher-order integration methods (Runge-Kutta 4 or 45) may be better at capturing features in the ODE/PDE examples, more testing is required to investigate their benefits for imaging classification.

5 Acknowledgements

The authors would like to acknowledge the support of AFOSR, FA9550-17-1-0125 and the support of NSF CAREER grant #1752116. We would like to thank Scott McCalla for providing feedback on this manuscript.

References

[1] Fashion MNIST Dataset. https://github.com/zalandoresearch/fashion-mnist.

[2] Fashion MNIST ResNet18. https://github.com/kefth/fashion-mnist.

[3] MLP for Fashion MNIST. https://github.com/heitorrapela/fashion-mnist-mlp.

[4] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.

[5] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems, pages 6571– 6583, 2018.

[6] Emmanuel de Bezenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical pro- cesses: Incorporating prior scientific knowledge. arXiv preprint arXiv:1711.07970, 2017.

[7] Weinan E. A proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics, 5(1):1–11, 2017.

[8] Amir Gholami, Kurt Keutzer, and George Biros. ANODE: Unconditionally accurate memory- efficient gradients for neural ODEs. arXiv preprint arXiv:1902.10298, 2019.

[9] C Lee Giles and Tom Maxwell. Learning, invariance, and generalization in high-order neural networks. Applied optics, 26(23):4972–4978, 1987.

[10] Eldad Haber and Lars Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34(1):014004, January 2017.

[11] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. ArXiv e-prints, December 2015.

[12] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Identity mappings in deep residual networks. ArXiv e-prints, March 2016.

[13] Philip Holmes, John L Lumley, Gahl Berkooz, and Clarence W Rowley. Turbulence, coherent structures, dynamical systems and symmetry. Cambridge university press, 2012.

[14] Gao Huang, Zhuang Liu, Laurens van der Maaten, and Kilian Q. Weinberger. Densely con- nected convolutional networks. ArXiv e-prints, August 2016.

[15] Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. arXiv preprint arXiv:1502.03167, 2015.

[16] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.

[17] J. Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua L. Proctor. Dynamic Mode Decomposition: Data-driven modeling, Equation-free modeling of Complex systems. SIAM, 2016.

[18] J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor. Dynamic mode decomposition: data-driven modeling of complex systems. SIAM, 2016.

[19] Gustav Larsson, Michael Maire, and Gregory Shakhnarovich. FractalNet: Ultra-deep neural networks without residuals. ArXiv e-prints, May 2016.

[20] Yann LeCun, L´eon Bottou, Yoshua Bengio, Patrick Haffner, et al. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 1998.

[21] Zichao Long, Yiping Lu, and Bin Dong. Pde-net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network. arXiv preprint arXiv:1812.04426, 2018.

[22] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning PDEs from data. arXiv preprint arXiv:1710.09668, 2017.

[23] Yiping Lu, Aoxiao Zhong, Quanzheng Li, and Bin Dong. Beyond finite layer neural networks: Bridging deep architectures and numerical differential equations. arXiv preprint arXiv:1710.10121, 2017.

[24] Dougal Maclaurin, David Duvenaud, and Ryan P Adams. Autograd: Reverse-mode differen- tiation of native python. In ICML workshop on Automatic Machine Learning, 2015.

[25] Tong Qin, Kailiang Wu, and Dongbin Xiu. Data driven governing equations approximation using deep neural networks. arXiv preprint arXiv:1811.05537, 2018.

[26] Maziar Raissi and George Em Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 357:125–141, 2018.

[27] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Machine learning of linear differential equations using Gaussian processes. Journal of Computational Physics, 348:683– 693, 2017.

[28] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations. arXiv preprint arXiv:1711.10566, 2017.

[29] Samuel H Rudy, Steven L Brunton, Joshua L Proctor, and J Nathan Kutz. Data-driven discovery of partial differential equations. Science Advances, 3(4):e1602614, 2017.

[30] Lars Ruthotto and Eldad Haber. Deep neural networks motivated by partial differential equa- tions. ArXiv e-prints, April 2018.

[31] Tim Salimans and Durk P Kingma. Weight normalization: A simple reparameterization to accelerate training of deep neural networks. In Advances in Neural Information Processing Systems, pages 901–909, 2016.

[32] Hayden Schaeffer. Learning partial differential equations via data discovery and sparse opti- mization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2197):20160446, 2017.

[33] Hayden Schaeffer and Scott G McCalla. Sparse model selection via integral terms. Physical Review E, 96(2):023302, 2017.

[34] Hayden Schaeffer, Giang Tran, and Rachel Ward. Learning dynamical systems and bifurcation via group sparsity. arXiv preprint arXiv:1709.01558, 2017.

[35] Hayden Schaeffer, Giang Tran, and Rachel Ward. Extracting sparse high-dimensional dynamics from limited data. SIAM Journal on Applied Mathematics, 78(6):3279–3295, 2018.

[36] Peter Schmid and Joern Sesterhenn. Dynamic Mode Decomposition of numerical and experi- mental data. Bulletin of the American Physical Society, 53, 2008.

[37] Peter J Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics, 656:5–28, 2010.

[38] Yoan Shin and Joydeep Ghosh. The pi-sigma network: An efficient higher-order neural network for pattern classification and function approximation. In IJCNN-91-Seattle International Joint Conference on Neural Networks, volume 1, pages 13–18. IEEE, 1991.

[39] Giang Tran and Rachel Ward. Exact recovery of chaotic systems from highly corrupted data. Multiscale Modeling & Simulation, 15(3):1108–1129, 2017.

[40] Linan Zhang and Hayden Schaeffer. Forward stability of resnet and its variants. arXiv preprint arXiv:1811.09885, 2018.

[41] Linan Zhang and Hayden Schaeffer. On the convergence of the sindy algorithm. Multiscale Modeling & Simulation, 17(3):948–972, 2019.

A Derivation of Adjoint Equations

Let be the vector of learnable parameters (that parameterizes the unknown function g, which embeds all network features), then the training problem is:

where 0 are regularization parameters set by the user. All subscripts with respect to a variable denote a partial derivative. We use the dot-notation for time-derivative for simplicity of exposition. The function r is a regularizer on the parameters (for example, the norm) and the time-derivative is penalized by the norm. Define the Lagrangian L by:

where ] is a time-dependent Lagrange multiplier. To apply gradient-based algorithms, the total derivative of the Lagrangian with respect to the trainable parameter must be calculated. Using integration by parts after differentiating with respect to

=

The initial condition ) is independent of any two time-stamps [

then

To determine ) and at the right-endpoints of [)). The derivative of the Lagrangian with respect to

We augment the evolution for ), starting at and integrating backwards. The code follows the structure found in [5].

Designed for Accessibility and to further Open Science