Accelerating Neural ODEs Using Model Order Reduction

Embedding nonlinear dynamical systems into artificial neural networks is a powerful new formalism for machine learning. By parameterizing ordinary differential equations (ODEs) as neural network layers, these Neural ODEs are memory-efficient to train, process time series naturally, and incorporate knowledge of physical systems into deep learning (DL) models. However, the practical applications of Neural ODEs are limited due to long inference times because the outputs of the embedded ODE layers are computed numerically with differential equation solvers that can be computationally demanding. Here, we show that mathematical model order reduction (MOR) methods can be used for compressing and accelerating Neural ODEs by accurately simulating the continuous nonlinear dynamics in low-dimensional subspaces. We implement our novel compression method by developing Neural ODEs that integrate the necessary subspace-projection and interpolation operations as layers of the neural network. We validate our approach by comparing it to neuron pruning and singular value decomposition (SVD)-based weight truncation methods from the literature in image and time-series classification tasks. The methods are evaluated by acceleration versus accuracy when adjusting the level of compression. On this spectrum, we achieve a favorable balance over existing methods by using MOR when compressing a convolutional Neural ODE. In compressing a recurrent Neural ODE, SVD-based weight truncation yields good performance. Based on our results, our integration of MOR with Neural ODEs can facilitate efficient, dynamical system-driven DL in resource-constrained applications.


I. INTRODUCTION
Deep Learning (DL) is reaching and surpassing human performance in domain-specific applications [1].Accordingly, there is increased demand for including DL-based algorithms into consumer devices that may contain only limited computational capacity and be battery-powered, making resource efficiency of DL-based algorithms important.An interesting new development in DL research is Artificial Neural Networks (ANNs) that employ dynamical systems, replacing traditional discrete layers with a continuous-time layer in the form of ordinary differential equations (ODEs) [2], [3], [4], [5].In these Neural ODEs, the continuous formalism allows flexible processing of time-series and irregularly sampled data, such as medical and physiological signals [6].Moreover, Neural ODEs are useful for resource-constrained and embedded applications Mikko Lehtimäki and Marja-Leena Linne are with the Faculty of Medicine and Health Technology, Tampere University, Finland.Lassi Paunonen is with the Faculty of Information Technology and Communication Sciences, Tampere University, Finland.
Correspondence: mikko.lehtimaki@tuni.fibecause they are very memory and parameter efficient [5].The ODE layer also facilitates engineering physical details, such as energy conservation laws or spectral properties, into neural networks [4].However, often a big computational bottleneck in Neural ODEs is the dynamical system layer, since propagating data through the system requires many evaluations using numerical ODE solvers.Reducing the computational cost of the ODE layer is the main motivation of our work.
In this work we show that Neural ODEs can be accelerated by compressing them using model order reduction (MOR) methods.The MOR approach is based on projecting dynamical systems onto low-dimensional subspaces.Here, we develop MOR methods that are integrated into Neural ODEs.In this manner, we lower the required storage size, memory consumption, multiply-adds and nonlinear activation count needed to compute predictions from input data.The resulting compressed Neural ODEs can be deployed for real-time computing and devices where energy efficiency is paramount.In order to validate the performance of our MOR method, we compare it to two established compression methods from the literature.Our results demonstrate that MOR is a theoretically grounded and effective method for accelerating Neural ODEs, because it achieves a favourable, adjustable and extensible balance between speedup and accuracy of compressed models.We show this result in two different classification tasks that use different Neural ODE architectures, a convolutional and a recurrent neural network (RNN).
Compressing ANNs is one of the principal ways of converting high-performing trained networks into more efficient networks, since good accuracy can often be recovered without long training times, while the size of the compressed network can be chosen in a flexible manner.Several neural network compression methods have been proposed for accelerating ANNs [7], [8], [9].These include singular value decomposition (SVD) based weight truncation [10], [11] and neuron pruning with an importance score based on zero activations [9], which we implement for comparison to our MOR method.Many prominent existing approaches rely on assigning importance scores to neurons based on their connection weights or output values.However, we argue that such methods are non-optimal for compressing Neural ODEs, as it is difficult to quantify the importance of neurons in the ODE layer based on criteria such as the output of the layer alone.This is because the final state of the nonlinear ODE system gives no information about the dynamics of the actual computations.In contrast, MOR methods are designed precisely for the approximation of state trajectories and inputoutput relationships in dynamical systems, hence they can overcome the shortcomings of existing compression methods when aiming to accelerate Neural ODEs.
Our approach to model acceleration is inspired by computational neuroscience.Due to high computational burden, modeling studies of the brain in realistic scales are limited to simulating small fractions of the brain using supercomputers [12].To overcome computational resource challenges, MOR methods have been adapted successfully for reducing single neuron [13], [14], synapse [15] and neuronal population [16], [17] models.Moreover, the similarities in models of the brain as well as artificial neural networks, which include connectivity patterns and nonlinear computation units, make us hypothesize that MOR methods may be a principled approach for accelerating Neural ODEs.
Our MOR approach for compressing Neural ODEs is based on the Proper Orthogonal Decomposition (POD) [18] with the Discrete Empirical Interpolation Method (DEIM) [19], a variant of the Empirical Interpolation Method [20].Using the POD-DEIM method, we derive reduced order models (ROMs) that can be simulated efficiently in low-dimensional subspaces, even in the presence of nonlinear activation functions.In the context of Neural ODEs, our method compresses the ODE block by projecting it onto a low-dimensional subspace using POD, which reduces the number of state variables in the ODE system.Linear operations are compressed by transformation into this subspace.When nonlinear activation functions are present, using DEIM we determine the most informative output neurons based on their time dynamics, and prune the other neurons.This reduces the dimensionality of the nonlinear operation and removes rows from the weight matrix, since connection weights of pruned neurons are discarded.We interpolate an approximate response for the pruned neurons directly in the low-dimensional subspace.In convolutional layers, the reduced model computes convolutions only at the selected interpolation points so that every kernel has its own set of evaluation coordinates.
We integrate the subspace projection and interpolation steps of POD-DEIM into the Neural ODEs as layers, and this introduces additional operations into the neural network (see Figure 2).In some network architectures these steps actually further reduce the overall number of multiply-adds in the neural network.The POD-DEIM reduction is applied after training the model and the reduced model can be fine-tuned for increased accuracy.In summary, Neural ODEs allow us to bridge a gap between ANN research and control theory research, so that a substantial amount of previously unemployed knowledge in model reduction can be utilized in ANN compression and acceleration.For example, analytical optimality results and error bounds exist for our MOR algorithms [18], [21].
In Section II we review previous work in compressing and accelerating neural networks in general and Neural ODEs specifically.In Section III we present our proposed model order reduction approach and show how to formulate it in the context of continuous neural networks.In the same section we introduce two established acceleration methods that we compare our method to in benchmark problems.In Section IV we provide theoretical compression ratios for the chosen methods and show actual accuracy and wall-time metrics of compressed Neural ODEs in two different classification tasks, one using a convolutional and the other using a recurrent ODE architecture.We discuss the success of the model reduction approach and the significance of our work in Section V with future suggestions and conclude in Section VI.

II. RELATED WORK
Several studies have reported that ANNs contain structures and redundancies that can be exploited for reducing memory and power usage [7].Here we focus on structural acceleration approaches that modify trained networks to achieve a more efficient architecture, and leave efforts such as low precision algebra, quantization of weights [22], binarization [23], hashing [24], vectorization [25], frequency space evaluation [26] and adjusting ODE solver tolerances or step size [5] out of the present study since those are complementary to the approaches presented here.Moreover, several hardware accelerators have been proposed (for example [27]), and those will not be addressed here.Furthermore, computational bottlenecks in ANNs have also been addressed by first reducing data dimension and then training simpler models.These techniques range from feature engineering to data dimensionality reduction.However, in this work our focus is on compressing trained networks, and hence data compression is considered complementary to out approach.Structural compression methods can be further grouped into several categories.a) Pruning weights: Prior work has addressed weight pruning [28], [29] and enforcing weight sparsity [30], [31] during training to obtain weight matrices that require less storage space and memory than dense weight matrices.Several methods to evaluate weight importance have been proposed, see for example a recent review addressing 81 pruning studies [7] and compares the achieved compression rates and speedups for several ANN models.It is common to include pruning in the training loop, because altering weights after training leads to accuracy loss.In order to achieve significant acceleration with weight-pruning methods the use of special software, masking strategies, sparse algebra or hardware accelerators is recommended [32], [33].
b) Decomposition: Decompositions and low-rank matrix factorizations have been used for compress network weights so that linear operations in a layer are computed efficiently [34], [10], [35], [36], [11].Decomposition is based on the observation that weight matrices, especially as their size increases, are seldom full rank.The idea is that a fully connected weight matrix or a tensorized convolutional kernel can be decomposed into compact low-rank matrices that require less storage space, memory and multiply-adds during training and testing.An approximation of the original operation is then obtained as a product of the low-rank tensors.In addition to compressing linear operations, prior work has also addressed decomposition by taking the nonlinear response into account [37].The cost of decomposition approaches is that a single layer is replaced with multiple smaller ones, which trades parallelism for forward operations, canceling out some of the obtained acceleration.
c) Pruning neurons: Weight pruning and decomposition approaches maintain the structure of the compressed network in that the number of nonlinear activation functions is not changed.Eliminating activation functions leads to acceleration as an entire row of weights from a fully connected layer can be removed [8], [38].As neuron pruning changes the number of neurons in a layer, the next layer must take this into account, for example by deleting columns (inputs) from the weight matrix in the case of fully connected layers.Alternatively, interpolation can be used to approximate the original response, possibly introducing additional computation.In fully connected architectures large compression rates in storage space and memory are obtained when rows or columns are deleted from fully connected layers.Importance scores, such as percent of zero activations have been developed for identifying prunable neurons [9].Pruning and quantization have been combined with Huffman coding into a compression framework that delivers memory and energy efficiency [39].An optimization approach has been used to enforce sparse columns in a fully connected layer so that the corresponding input neurons can be pruned [38].
d) Pruning filters: Convolutional Neural Networks (CNNs) have attracted a lot of attention in the compression literature.This is not surprising given their good rate of success in real-life tasks and their high computational cost.Convolutional operations make the bulk of modern image and video processing networks and are used in many other applications, such as sequence processing.It is possible to approach CNN compression by analysing either the convolutional filters or feature maps obtained by applying the filters on input data [40], [41], [42].The filters can be compressed with decomposition methods, similarly to fully connected layers.Alternatively, entire kernels can be pruned from filters.Pruning kernels is very effective as intermediate feature maps are eliminated.The number of methods and criteria available for identifying pruning targets in CNNs is very high [42].In convolutional networks, neuron pruning corresponds to skipping a kernel evaluation at a single spatial location.However, such an operation is rarely supported by deep learning frameworks, and pruning in convolutional networks has focused on eliminating parts of or entire filters or feature maps.
In the context of Neural ODEs, a few studies have addressed the acceleration of inference and training times.One approach focuses on learning simple dynamics that do not burden ODE solvers as much as stiff systems [43], [44], [45].Additionally, training times have been reduced by further improving the adjoint method, for example by relaxing error criteria [46] while inference times have been accelerated by introducing hypersolvers -neural networks that solve Neural ODEs [47].These methods are complementary to ours, since we aim at compressing the learned architecture with MOR methods.Our approach is comparable to pruning neurons based on an importance score [8], [9], so that a number of rows from the weight matrix and nonlinear activation functions can be removed altogether.The POD-DEIM method improves on the existing methods of determining neuron importance, as it accounts for the complete dynamics of the ODE block.Additionally, POD-DEIM can be combined with many of the Fig. 1.A plain discrete feedforward network on the left and a residual network on the right.The defining feature of residual networks is that the output from an earlier layer skips layers and is added directly to a later layer.
prior methods such as quantization of weights.

III. METHODS
In this section, we present the theory of the continuous interpretation of ANNs.We then describe our model order reduction method and show how to formulate it in the ANN and Neural ODE context.Finally, we present two other ANN acceleration methods from the literature, which we use as benchmarks to our method, and their implementation in the Neural ODE setting.

A. Continuous networks
Since the success of Residual Neural Networks (ResNets) [48], a continuous interpretation of neural networks has gained traction.The ResNet architecture utilizes skip connections to deal with the problem of vanishing gradients, allowing the training of extremely deep neural networks.In a ResNet block, the output is the sum of the usual feedforward operations on the data and the unprocessed data entering the block.The ResNet architecture is illustrated in Fig. 1 on the right, with a plain feedforward network on the left.If the hidden layers of the plain network implement a nonlinear function x k+1 = f (x k ), the skip connection of the ResNet implement x k+1 = x k + f (x k ).By introducing a constant h = 1 to obtain x k+1 = x k + hf (x k ), the resemblance of the skip connection to Euler's formula for solving Ordinary Differential Equations (ODEs) is seen [2], [4], [3].This has lead to the continuous-time dynamical system interpretation of ANNs.
In the continuous interpretation of ANNs, a set of layers is replaced by one layer that parameterizes a dynamical system as a group of ODEs with state variable x(t).We assume this ODE system has the form in general and in our models specifically where In Neural ODEs, the activation function f a (χ i ) is commonly the hyperbolic tangent, although any differentiable activation function, or a combination of different activation functions, can be used.Here, θ = (A, b, Z) are learned parameters of the ODE.In feed-forward Neural ODEs there are typically no time-dependent inputs and hence Z = 0. On the other hand, in RNNs where x (t) is the hidden state, the input data enters the system through Z = 0.The output of the layer is x(t end ), which is the state of the dynamical system at the user specified final time t end .Initial values of the ODE system can be obtained in two ways, either as the output of the layer preceding the ODE or set explicitly.The former is more common in feed-forward architectures and the latter in RNNs that receive time-dependent input data.The output x(t end ) of the ODE layer is computed using numerical methods, taking discrete or adaptive steps with an ODE solver to solve an initial value problem.Neural ODEs can use the adjoint method of calculating gradients [5] and have enabled parameterizing ODEs as several different ANN operations, or chains of them.The primary restriction is that the output size of the ODE layer must match the input-size of the layer.Overall, it is possible to train a variety of ANN architectures for different tasks as continuous networks [5].
Training Neural ODEs with the adjoint method is memory efficient compared to deep discrete networks [5], as the backpropagation algorithm does not need to store intermediate ODE states to calculate gradients on the loss function.However, the use of the adjoint method is not required for training ODE networks.Parameterizing an ODE in place of several discrete layers may also lead to parameter efficiency and correspondingly require less storage space and memory, since the continuous layer can replace several individually parameterized discrete layers.Other benefits of Neural ODEs include using ODE solvers for speed versus accuracy tuning, and enabling ANNs to flexibly process continuous and irregularly sampled time-series data [5], [6].Neural ODEs have also improved on existing ANN-based density estimation models [5].However, the dynamics learned by Neural ODEs may be unnecessarily complex [45], [44] and result in stiff systems [49].Often the ODE block is the most computationally demanding part of the network, since many numerical evaluations of the state of the ODE are needed to obtain the output.Hence, Neural ODEs require more time to evaluate data, both in training and testing, compared to traditional discrete networks [44].
It is possible to parameterize the ODE block so that the model in training has favourable properties that facilitate learning.Antisymmetric networks are a step toward this direction, since they guarantee that the state x(t) of the ODE system does not diverge far from or decay to the origin as t → ∞ [4].This prevents the gradient of the loss function from vanishing or exploding and makes training the network well-posed.Antisymmetric networks use an antisymmetric weight matrix A = W −W T that by definition has eigenvalues λ i so that for all i, Re(λ i (A)) = 0. Furthermore, a small shift of the eigenvalues by γ may be applied so that Re(λ i (A−γI)) = −γ < 0, which improves behavior of the ODE system in the presence of noisy data [4].Notice that in practical applications t is finite, hence a small γ will not make the gradients of the loss function vanish.In [50] it is demonstrated that the shifted antisymmetric weight matrix A − γI gives the hidden state of recurrent neural networks a favourable property of longterm dependence on the inputs to the system, which helps classifying data with temporal relationships.In this work, we implement an ODE-RNN as where the weight matrix A gives the network the desired properties that facilitate learning and during training the parameters W are learned.Other parameters are similar to Eq. ( 2).In Section IV, we will demonstrate compression and acceleration of this ODE-RNN.Such networks make an interesting model compression target since their architecture gives the system temporal memory capacity that the compressed model must retain.

B. Model Order Reduction
A key contribution of our work is the formalization of MOR in the context of Neural ODEs and the realization that the necessary MOR operations can be expressed as layers of ANNs.A powerful method for MOR of general nonlinear systems is Proper Orthogonal Decomposition (POD) [18] coupled with the Discrete Empirical Interpolation Method (DEIM) [19], [20], developed in the fields of systems and control theory.In order to compress Neural ODEs, we construct reduced order models (ROMs) of the ODE block in trained Neural ODEs with the POD-DEIM method.
Both POD and DEIM utilize the method of snapshots [51].In POD, snapshots are values of the state x t of the ODE system at discrete times t.A snapshot matrix X = [x 1 , x 2 , • • • , x s ] is collected using numerical simulation of the ODE block with different initial values and time-dependent input data.ANNs provide a very natural setting for gathering snapshots, since we have access to training data that can be propagated through the trained network and the ODE block of Neural ODEs.However, it is important that the snapshots are collected after the model is trained, so that the snapshots reflect true learned dynamics of the ODE layer.Moreover, the snapshots are not used for optimization or model training.The following explains the purpose of the snapshots.
POD approximates the original system of dimension n via projection using a low-dimensional subspace.A k-dimensional POD basis with orthonormal column vectors V k ∈ R n×k , where k < n, is computed using Singular Value Decomposition (SVD) of snapshots V ΣΨ T = SVD(X) [51].V k is then the first k left singular vectors of the snapshot matrix, equaling the first k columns of V .A reduced state vector V T k x(t) = x(t) ∈ R k is obtained by a linear transformation and at any time t an approximation of the original, fulldimensional state vector can be computed with x(t) ≈ V k x(t).Projecting the system in Eq. ( 2) onto V k by Galerkin projection results in a reduced system where we can precompute the transformation of the weight and input matrices into matrices solve (P T U )c = P T u l for c 5: U ← [U u l ], P ← [P e p l ], p ← [ p p l ] 7: end for used in place of the original weight matrices in the reduced models.The reduced POD model approximates the original system optimally in the sense that the POD subspace has minimum snapshot reconstruction error [18].However, the nonlinear form of the equation prevents computational savings as the number of neurons has not been reduced.The size of A is still n×k and f (•) is computed in the original dimension n.This is known as the lifting bottleneck in reducing nonlinear models.
Efficient evaluation of the nonlinear term can be achieved with DEIM [20], [19].DEIM extends the subspace projection approach with an interpolation step for general nonlinear functions.To construct an interpolated approximation of the nonlinear term, we use where the DEIM basis vectors < n are the first m left singular vectors of the snapshot matrix of nonlinear vector-valued function outputs ) is a nonlinear function with m components chosen from f according to DEIM determined interpolation points p = p 1 , • • • , p m and P m = [e p1 , e p2 , • • • , e pm ] with e pi being the standard basis vector i of R n .Together POD and DEIM form a ROM where N ∈ R k×m can be precomputed.Now in the prediction phase only m nonlinear functions are evaluated.Due to the structure of f m : R m → R m we can then further compress Ã so that only m rows at indices p remain.Correspondingly, we only select m elements from the bias vector at indices p.
Thus we have obtained the reduced order model where Ãm ∈ R m×k and b m ∈ R m .In Eq. ( 6), the dimension k of the POD projection subspace does not need to equal the number m of the DEIM interpolation points, as N can be computed even if k = m, although in practice it is common to choose m = k.This means that m, k can be chosen either smaller or larger with respect to each other to reflect the complexity of linear and nonlinear functions of the model.Algorithm 1 shows the process of determining the DEIM interpolation points and matrix P .The ordered linearly independent basis vector set given as input is the basis U m that is obtained from snapshots of the nonlinear function.The argmax-function returns the index of the largest value in a given vector and e i is the i:th standard basis vector.For more details see the original reference [19].Therein an error bound for the DEIM approximation in Euclidean space is also given.For an error estimate of the reduced state in POD-DEIM models specifically, see [21].

ODE ODE
An essential part of our work is building the POD-DEIM subspace projection operations V T k , V k , reduced matrices Ãm , Z and low-dimensional interpolation N of Eq. ( 7) into the Neural ODE.Consider a Neural ODE with an input layer, an ODE block, and a readout layer, as on the left of Fig. 2. We add two new layers, one for projecting the ODE block into a low-dimensional subspace and one for transforming the result of the ODE block back to the original space, as seen on the right of Fig. 2.These layers implement the computations V T k x(0) and V k x(t end ), respectively.The cost of these operations is offset by the cheaper evaluation of the ODE layer, which typically makes the bulk of the compute time of Neural ODEs.In the ODE layer, the reduced matrices Ãm , Z replace the corresponding original large weight matrices and a new linear layer is added to implement the low-dimensional interpolation N .Furthermore, if the operations around the ODE block are linear, V T k , V k can be computed into those existing layers, resulting in even cheaper online evaluation.Finally, we emphasize that this POD-DEIM process is widely applicable to different model architectures, as long as Eq. ( 1) defines the ODE block.
The DEIM method is known to exhibit lower accuracy in the presence of noise or perturbations in the snapshots.The input function Zu(t) in ODE-RNNs can cause such non-smooth behavior and a large number of timesteps taken is sensitive to instabilities.Hence for compressing ODE-RNNs, we use an oversampling strategy that stabilizes the method as shown in [52].In this oversampled DEIM (ODEIM), the number of nonlinear sampling points in p becomes m+o and is decoupled from the number of basis vectors m of the DEIM subspace.The matrix P then has size n×(m+o) and the matrix inverse of Eq. ( 6) is replaced with the Moore-Penrose pseudoinverse (P T U m ) † .In ODEIM, interpolation becomes approximation via regression.During model evaluation, m + o activation functions are evaluated.In our results, we evaluate m + 1 nonlinear activation functions, which compared to vanilla DEIM has a slightly negative effect on model acceleration but a positive effect on accuracy of the compressed model.
The POD-DEIM method has adjustable parameters that affect the performance of the reduced model.The dimension k of the linear projection matrix and the number of interpolation points m are typically set according to the decay of the singular values of the snapshot matrices X, F , but both can be adjusted independently.The number of snapshots can be controlled in two ways, the number of samples of training data to use and the number of snapshots of model dynamics to save per sample.Computing the SVD of the snapshot matrix is the most computationally demanding step of the MOR pipeline, with memory being the typical bottleneck.

C. Benchmark compression methods
We compare the performance of the POD-DEIM method to two other ANN compression methods presented in the literature.To the best of our knowledge, no other studies using these methods to compress Neural ODEs have been published before.
The first method compresses the weight matrix of the ODE block into two low-rank matrices [10], [11] using SVD based truncation.This reduces the overall number of multiply-adds needed to evaluate the layer but, unlike POD-DEIM, does not change the number of activation functions.Given a fully connected layer with activation f (Ax + b) and A ∈ R n×l , we apply SVD to the weight matrix to obtain ΦΣΨ T = A.Then, we only keep the first k < n singular values by truncation so that Σ k = diag(σ 1 , . . .σ k ) ∈ R k×k and the first k left and right singular vectors so that In ANNs we can implement the approximation in two consequent fully connected layers that only use activation and bias after the last layer.The connection weights of the first layer are Σ k Ψ T k ∈ R k×l and the second weights are Φ k , which changes the parameter count from nl to k(n + l).Notice that in a single-layer ODE block we have l = n.This method is theoretically sound for compressing the ODE block.It retains the most parameters out of the compression methods we test, meaning that it is expected to yield low to moderate acceleration.Since the weight truncation only uses the learned weights as inputs, the method does not depend on any other data, unlike POD-DEIM which requires gathering snapshots of the ODE dynamics.
The second method we used for comparison is based on computing an importance score for neurons, called average percent of zeros (APoZ) method, following [9].There it is suggested to prune neurons that have the largest number of zero activation values over the training or testing data set.The corresponding rows of the weight, bias and input matrices are then removed altogether.Moreover, with Neural ODEs specifically, for each removed row of the weight matrix we also remove the corresponding column.The importance score of the c:th neuron in a layer is defined as where N is the number of examples used to calculate the score and M is the dimension of the output feature map O.In Neural ODEs it is common to use hyperbolic tangent activations, hence we adapt the scoring so that neurons with lowest absolute activation magnitudes are eliminated.We compute the score after the last timestep of the ODE block, as this is the value that gets propagated to the following layers in the network.We only use training data for computing the score.This method does not account for the dynamics of the ODE block, as scoring is computed only at the last timestep of the ODE block, hence we evaluate the performance in detail here.
Overall, this method applies the largest amount of compression out of the methods we test here.
IV. RESULTS We implemented two Neural ODE models and trained them on different classification tasks; one Neural ODE with a convolutional ODE block to classify digits of the MNIST dataset and one ODE-RNN to classify digits in the pixel-MNIST task, where the network is given one pixel at a time as input.We trained both models using data-augmentation with random rotations and affine translations in order to prevent overfitting, using a decaying learning rate starting from 0.04, stochastic gradient descent optimizer and cross-entropy loss function on image labels.After training, the POD-DEIM snapshots were collected by feeding each image in the training data to the network, and saving the ODE state and activation function values at every second timestep.POD and DEIM bases were computed using a memory efficient partitioning strategy [53].Next, we compressed the models using POD-DEIM, APoZ trimming and weight truncation and evaluated their performance directly after compression, after 3 epochs of fine-tuning and after excessive fine-tuning on training data.The performance of each compressed model was compared to the respective original model using Top-1 accuracy, Top-3 accuracy and wall-time as metrics.We restrict fine-tuning to the layers following the ODE block in order to get a better view of how well the compressed blocks maintain their computational capacity.Our results below have been implemented using the PyTorch machine learning framework [54].To implement Neural ODEs, we additionally used the TorchDiffEq [5] and Torchdyn [55] python packages.Execution times are true walltimes realized on Intel Xeon E5-2680 v3 2.5GHz processors and measured as the time it takes to classify every item in the test dataset.For the convolutional network, we use a single core, and for the recurrent network four processor cores.
In Table I, we show the number of parameters in a theoretical Neural ODE and after it has been compressed with each compression method.We consider a nonlinear ODE block of the form x (t) = f (Ax(t)) where the weight matrix has dimensions A ∈ R n×n , hence there are n activation functions.The ODE block is preceded by a layer with weight matrix W 1 ∈ R n×i and followed by a layer with weight matrix operations compared to APoZ.Moreover, for POD-DEIM the best case is realized when the preceding and following operations are linear, whereas the APoZ method reaches the best case even if they are nonlinear.With the APoZ method, some nonlinear operations, such as max pooling, require that we maintain a lookup table of pooled indices or insert the k compressed values to their respective indices in an n size vector, which is the worst case result for the method.SVD truncation is not affected by the surrounding layers.

A. Convolutional Neural ODE
We illustrate the performance of reduced models on the MNIST dataset of handwritten digits [56].For our first task we train a convolutional Neural ODE that uses 16 convolutional kernels and tanh-activation in the ODE block, with Z = 0 (see Eq. ( 2).Before the ODE block we use a convolutional layer with 3 × 3 kernels and rectified linear unit (ReLU) activations that outputs 16 channels into 3 × 3 max pooling.After the ODE block we use 3 × 3 max pooling into a linear readout layer with ReLU activations.The ODE block propagates the data for t = [0, 1] with timestep dt = 0.1 using a fourth order fixed-step Runge-Kutta solver.The initial value of the ODE is the output from the prior layer.We use the adjoint method for training the ODE weights [5].In order to implement model order reduction and model compression of the convolutional ODE block, the convolution operator is needed in matrix form.After training, we replaced the original convolution operator with an equivalent operation implemented as a linear layer.The conversion yields a sparse Toeplitz matrix.With this operation, 16 convolution channels correspond to 1024 neurons and activation functions.All reported wall times are obtained with this conversion and compressed models are compared to these times.The reported wall times are medians of 10 consecutive evaluations of the entire validation data set.The trained network achieves 97.9% top-1 and 99.7% top-3 accuracies on held-out test dataset of 10 000 MNIST images and it takes 9.4 seconds to run the model on the test dataset on CPU.
In order to determine the suitability of our chosen compression methods, we looked at several metrics of the trained Neural ODEs, namely decay of singular values in the POD-DEIM snapshot matrices X, F and the weight matrix A. Moreover, we looked at the distribution of APoZ scores and their relation to DEIM interpolation indices.Fig. 3 shows these metrics for the convolutional Neural ODE.The left plot shows sorted singular values of the solution snapshot matrix X (POD) and F (DEIM) of the ODE block, gathered from training data, in purple and green colors, respectively.The logarithmic y-axis shows the magnitude of a given singular value divided by the sum of all singular values.Both POD and DEIM singular values collapse rapidly, indicating that a small number of singular vectors can span a space where the majority of dynamics are found.This is an important indicator for the suitability of the POD-DEIM method.In the same plot, the yellow plot shows the singular values of the weight matrix A of the ODE block.These values are useful for determining a rank for the SVD based truncation method that we compare against POD-DEIM.Singular values of A also decay exponentially, meaning that a good lowrank approximation to A can likely be found.The right plot of Fig. 3 shows APoZ scores for each neuron in the ODE block, with low-scoring neurons being the least important.Additionally, we have indicated the first 100 interpolation points p i chosen by the DEIM method in green color.All neurons seem to contribute to the model, with only a small number of relatively low-scoring neurons seen, which could make the APoZ method challenging to apply in the Neural ODE context.Some consensus between the methods is seen, as neurons with low APoZ scores are not in the top DEIM selections either.
We wanted to find out how the POD-DEIM method compares to existing model compression methods when compressing only the ODE block of a convolutional network.For POD-DEIM dimensions, we have chosen k = m and use dimension k compressed models, which corresponds in number of activation functions to dimension k APoZ trimmed models.The results on the MNIST dataset are seen in Fig. 4, where y-axis, model accuracy, is computed as reduced model top-1 accuracy divided by original model result and x-axis, achieved speedup, as original model test time divided by reduced model test time.Each dot represents a compressed model dimension, which decreases to the direction of increased speedup.While the dimensions are not comparable between methods, the graphs indicate how much accuracy is maintained at a given acceleration amount.Table A.II presents the absolute performance values of compressed models without fine-tuning, with three epochs of fine-tuning after model compression and with retraining (30 epochs of training after compression).
Fig. 4 shows that the POD-DEIM method retains the highest accuracy when comparing to similar speedups from other methods.Our POD-DEIM reaches over four fold speedup in runtime in this task.The SVD method is also accurate, but cannot reach speedups greater than two times the original runtime, explained by the method having to evaluate all the original nonlinear activation functions.The APoZ method reaches the fastest run times but loses the most amount of accuracy.After three epochs of fine-tuning, the APoZ method improves in accuracy, but not to the level of performance that the other methods show even without tuning.Overall, the POD-DEIM method achieves the best profile in terms of accuracy retained for speedup gained.A big contributor is the interpolation matrix N , which makes the method slower compared to APoZ, but the cost is justified by the better accuracy.When looking at absolute performance values, Table A.II shows that at dimension 350, the POD-DEIM method is two times faster to run than the original model, while achieving 93.1% top-1 accuracy.At dimension 50, the POD-DEIM model is over four times faster than the original model, still reaching 91.1% top-1 accuracy without fine-tuning.

B. Recurrent Neural ODE
For our second study, we designed a task in which Neural ODEs have an advantage over discrete networks.ODE-RNNs [6] and antisymmetric RNNs [4], [50] are one such example.Using a continuous hidden state, ODE-RNNs enable flexible inference on continuous-time data even with irregularly sampled or missing values and antisymmetric models have spectral properties that are favourable for RNNs, as described in Section III.We implemented a shifted antisymmetric ODE-RNN for the pixel MNIST task.In this task, the network sees a pixel of a handwritten digit at a time and after iterating over each pixel outputs the class of the image.The length of the input sequence is 784 steps, corresponding to a flattened MNIST image, and the network sees each pixel for dt = 0.1 seconds.Following [50], we initialize the ODE weights from zero mean, unit variance normal distribution as well as train and test the network with the Euler discretization method.We use 512 hidden units with a hyperbolic tangent nonlinearity and the shifted antisymmetric weight matrix as shown in Eq. 3. Initial value of the ODE is set to x(0) = 0.The ODE block is followed by a linear readout layer.Compared to the convolutional network, this model is smaller in neuron count but needs to incorporate the time-dependent input data.The trained network achieves 96.2% accuracy on held-out MNIST test data in 48.1 seconds.Fig. 5 shows compression metrics for the ODE-RNN.The left plot shows sorted relative singular value magnitude of the solution snapshot matrices X (POD) and F (DEIM) and the weight matrix A of the ODE block in purple, green and yellow, respectively.For this model, the singular values of the snapshot matrices as well as the weight matrix also decay rapidly.However, with regards to POD-DEIM, the singular values of the snapshot matrices do not decay as fast as those obtained from the convolutional model and there is a considerable amount of total energy contained in the singular values until rank 500.This indicates that the dynamics of the ODE block are more difficult to approximate in low dimensions than those learned by the convolutional model.Such a result is expected, because the ODE-RNN has a time-dependent input contributing to the hidden state and the forward propagations are computed for a longer time span than in our earlier example, allowing for richer dynamics.On the other hand, singular values of the shifted antisymmetric weight matrix decay more rapidly than in the earlier task, indicating suitability of SVD-based truncation of weights.The APoZ scores on the right of Fig. 5 show the same pattern as earlier, where it is not straightforward to detect the most important neurons conclusively.This is expected, given that in the context of ODE blocks, the output activity at last time step is not a conclusive measure of neuron importance.
We then evaluated the selected compression methods on the pixel MNIST task using our ODE-RNN.The results relative to the performance of the original model are seen in Fig. 6.Absolute results of compressed models without finetuning, with three epochs of fine-tuning and with exhaustive tuning in Table A.III.This task is more challenging for the reduction methods than the convolutional network.Possible factors are the dense structured antisymmetric weight matrix or the time-dependent input.Here, the SVD truncation method maintains the most accuracy.As predicted by the singular values of the snapshot matrices, the overall performance of the POD-DEIM method on this model is not as satisfying as with the convolutional model.On this model the POD-DEIM method is faster than the original model until dimension 425, at which point the accuracy is 94.6%, while the APoZ method remains faster at all times.With regards to the tradeoff between run time and accuracy, the SVD method maintains the most accuracy for speed gained.POD-DEIM shows good accuracy at high dimensions, although accuracy drops rapidly so that the overall performance is level with the APoZ method.With exhaustive fine-tuning the APoZ method may be fit for compressing the ODE-RNN to a low dimension, since APoZ yields the best absolute acceleration.The APoZ compressed model also seems easier to fine-tune than the POD-DEIM version, since more accuracy is recovered in three fine-tuning epochs.

V. DISCUSSION
In this work we have studied compression of artificial neural networks that contain ordinary differential equations as layers (Neural ODEs).After training and compressing the networks, we evaluated their accuracies as well as execution speeds.Although in the literature many compression methods are applied in the training loop, the training of Neural ODEs is slower than discrete networks to begin with, making it attractive to compress the network in a separate step after training.We also assessed fine-tuning of the compressed network to see how much accuracy can be recovered.We focused on the runtime and classification accuracy of compressed models, since Neural ODEs are already very memory and parameter efficient by design.
Here we showed how to formulate the Proper Orthogonal Decomposition with Discrete Empirical Interpolation Method (POD-DEIM) into compressing Neural ODEs.This method has been developed for accurately approximating dynamical systems in low dimensions and the continuous-time network formalism of Neural ODEs makes it possible to use POD-DEIM in deep learning applications.Additionally, we compared the POD-DEIM and APoZ methods to Singular Value Decomposition (SVD) based truncation.The performance was measured with two different architectures, a convolutional Neural ODE and a recurrent ODE architecture, which were trained on the MNIST digit classification task.
We hypothesized that some existing neuron pruning methods that measure neuron importance based on the final output of the layer only, such as the average percent of zeros (APoZ) approach, are not optimal for compressing Neural ODEs.The reason is that these methods do not account for the dynamics of the ODE block when identifying pruning targets.Based on our results, the proposed approach for Neural ODE compression depends on the model architecture.We found the POD-DEIM method to achieve the most favourable continuum between accuracy and speedup on compressing the convolutional architecture, as seen in Fig. 4. The SVDbased weight truncation is the simplest of the studied methods and it performed well on both tasks, although it yielded less absolute acceleration than the other methods.The ODE-RNN was challenging to accelerate, and here the weight truncation method is our favoured choice based on Fig. 6.In this task, the APoZ method performed better or equally to POD-DEIM.While the POD-DEIM method is based on approximating the time-trajectory in a low-dimensional subspace, there are also MOR methods, such as the Balanced Truncation, that are designed to approximate input-output behavior of highdimensional systems.Such a method could be a good tool for compressing RNN models with time-dependent input, once the applicability of these methods to nonlinear systems develops further [57].
In this study we only fine-tuned the layers following the ODE block, keeping results more indicative of raw performance of the methods.If the entire model were tuned, POD-DEIM and SVD truncation have more parameters available for fine-tuning than the APoZ method, allowing for greater accuracy.On the other hand, the small number of parameters  is what makes APoZ compressed models the fastest out of the three methods, hence it should be favored if there are no restrictions on fine-tuning and retraining time.Moreover, compared to POD-DEIM, the APoZ method is more data and memory efficient, since it does not need snapshots of trajectories of ODE dynamics nor does it compute large matrix decompositions.Overall, we found that POD-DEIM can achieve good acceleration and accuracy of compressed models, it has a solid theoretical foundation and several extensions, while also being the most adjustable and optimizable approach of the methods used in this study.
We expect our results with the POD-DEIM method motivate machine learning researchers to look into the possible applications of other MOR methods developed in control theory and related fields.For example, we assessed here only the original DEIM [19] and the ODEIM [52] methods, although in recent years many advanced versions have been developed.These include novel interpolation point selection methods [58] and methods with adaptive basis and interpolation point selection capabilities [59], [60], [61] that allow the model to change its dimension or projection subspace online.It is worth noting that the POD-DEIM method aims at approximating the dynamics of the original system and is not aware of classification or other performance indicators of the Neural ODE.A potential future development could explore connecting the low-dimensional dynamics to downstream network performance.
A limitation of the present work is that we could only study relatively small models, due to known challenges with training very large Neural ODEs [45].In small models the expected gains from compression methods are lesser than in large models.Moreover, we only discussed computation times on CPUs, as on GPUs the differences between original and compressed models were minor due to the small size of the models.While the computational speedup from POD-DEIM needs to be verified on GPUs and extremely large models, theoretical reductions in parameter and activation function count, comparable results from existing literature as well as our results on the CPU indicate that the method should perform equally well.
Our results were obtained with fixed-step ODE solvers.We found that compressing Neural ODEs sometimes induces stiffness to the low-rank models (data not shown).This could cause adaptive ODE solvers to be slower, mitigating some of the speedup from compression.However, a systematic study of this behavior was not sought here and fixed-step ODE solvers were chosen to quantitatively compare the methods across dimensions.A future study is needed to quantify the effects of stiffness that may rise when compressing Neural ODEs.

VI. CONCLUSION
In this study, we addressed speeding up and compressing artificial neural networks with continuous layers.Our POD-DEIM method, which has not previously been used to compress Neural ODEs, provides attractive results for accelerating convolutional Neural ODEs and equal results to neural pruning methods for accelerating recurrent Neural ODEs.POD-DEIM properly considers the trajectory of the ODE layer in an efficient manner.Additionally, POD-DEIM has adjustable parameters and advanced variations that can be used to optimize the speed versus accuracy trade-off for each model.
We expect our acceleration method to be useful when Neural ODEs are deployed on low-power hardware or batterypowered devices, for example wearable devices, medical monitoring instruments or smart phones.The dynamical system formalism for neural networks also allows other model order reduction methods to be applied in deep learning, which will provide for interesting future studies.Each column of Top-1, Top-3 and Runtime shows the results for the models directly after reduction, after three fine-tuning epochs and after retraining (30 epochs), in that order, separated by vertical lines.Each column of Top-1, Top-3 and Runtime shows the results for the models directly after reduction, after three fine-tuning epochs and after retraining (30 epochs), in that order, separated by vertical lines.

Fig. 2 .
Fig. 2. A Neural ODE on the left with the discretized differential equation block in orange color.A POD-DEIM reduced network on the right illustrates that the ODE block is evaluated in a low-dimensional subspace, with linear transformations around the ODE block.The networks have equal inputs and approximately equal outputs.

Fig. 3 .Fig. 4 .
Fig. 3. Metrics for the convolutional Neural ODE.Left: Singular values POD (purple) and DEIM (green) snapshot matrices and singular values of the weight matrix of the ODE block (yellow).Right: histogram of APoZ scores, with the first 100 DEIM indices indicated in green color.

Fig. 5 .
Fig. 5. Metrics of the ODE-RNN.Left: Singular values of POD (purple) and DEIM (green) snapshot matrices and singular values of the weight matrix of the ODE block (yellow).Right: histogram of APoZ scores, with the first 50 DEIM indices indicated in green color.

Fig. 6 .
Fig. 6.Relative performance of reduced ODE-RNNs on the pixel MNIST task.Left plot shows results without fine-tuning, and right plot with three epochs of tuning.The x-axis shows achieved acceleration, while the y-axis shows reduced model top-1 accuracy divided by full model accuracy.Reduced model dimension decreases to the direction of improving speedup.Absolute values can be seen in Table A.III.
2 ∈ R o×n .In TableI, k denotes the dimension of the compressed ODE block and n × i, o × n are the number of weights in the layers preceding and following the ODE block, assuming those layers exist.With regards to the model speed, the complexity of the ODE layer is significantly more important than the sizes of the surrounding layers.This is because the ODE function is evaluated multiple times; at least once for each intermediate time point in [t 0 , t end ].The architecture of the original ANN affects the achieved compression.The POD-DEIM and APoZ methods can reduce the ODE layer to a dimension that is independent of the original size n, although the interpolation in POD-DEIM requires an additional k 2 W

TABLE II ABSOLUTE
PERFORMANCE OF REDUCED CONVOLUTIONAL NEURAL ODES.

TABLE III ABSOLUTE
PERFORMANCE OF REDUCED ODE-RNNS.