Neural DAEs: Constrained neural networks
Authors: Tue Boesen
Eldad Haber
Uri M. Ascher
Venue: Submitted to ICLR 2023 (2023)
Keywords: neural networks
differential algebraic equations
constraints

Abstract: In this article we investigate the effect of explicitly adding auxiliary trajectory information to neural networks for dynamical systems. We draw inspiration from the field of differential-algebraic equations and differential equations on manifolds and implement similar methods in residual neural networks. We discuss constraints through stabilization as well as projection methods, and show when to use which method based on experiments involving simulations of multi-body pendulums and molecular dynamics scenarios. Several of our methods are easy to implement in existing code and have limited impact on performance while giving significant boosts in terms of inference.

NEURAL DAES: CONSTRAINED NEURAL NETWORKS
Anonymous authors
Paper under double-blind review
ABSTRACT
In this article we investigate the effect of explicitly adding auxiliary trajectory information to neu- ral networks for dynamical systems. We draw inspiration from the field of differential-algebaai equations and differential equations on manifolds and implement similar methods in residual neu- ral networks. We discuss constraints through stabilization well projection methods, and show when to use which method based on experiments involving simulations of multi-body pendulums and molecular dynamics scenarios. Several of our methods are easy to implement in existing code and have limited impact on training performance while giving significant boosts in terms of infer- ence.
1 INTRODUCTION
Many scientific simulations of dynamical systems have natural invariants that can be expressed by constraints. Such constraints represent conservation of some quantities of the system under study. For example, in molecular dynamics, bond lengths between atoms are assumed fixed. Another example is incompressible fluid flow, where the divergence of the velocity field vanishes at any point in space and time. Similarly, in Maxwell's equations, the divergence of the magnetic field vanishes (no magnetic charge). Such additional information about the flow can be crucial if we are to keep the simulations faithful. As result, wealth of techniques have been proposed to conduct simulations that obey the constraints at least approximately (Weiglhofer 1994; Ascher & Petzold [998 Allen et al. 2004).
In recent years, machine learning based techniques, and in particular deep neural networks, have been taking a growing role in modelling physical phenomena. In some cases, such techniques are used as inexpensive surrogates of the true physical dynamics and in others they are used to replace it altogether (see e.g. Wang al. (2018); Degiacomi (2019); Miyanawala & Jaiman (2017)). These techniques use the wealth of data, either observed or numerically generated, in order to "learn" the parameters in neural network, so that the data are fit to some accuracy. The network then expected to perform well on new data, outside of the training set and yield simulation results that are accurate and reliable, in many cases, at significantly less computational effort.
From classical simulations, we know that it is often as important to accurately obey the additional constraint informa- tion as it to accurately satisfy the underlying ordinary/partial differential equation (ODE/PDE) system. Nonetheless, no neural network architecture known to us is designed to honour such constraints or invariants. The hope in standard training procedures is that by fitting the data, the network will "learn" the constraints and embed them in the weights implicitly. This, however, has been demonstrated to be insufficient in many cases (Wah & Qian 2001). As we show in this paper, on some very simple examples, neural networks may able to approximately learn the dynamics but they can drift off the constraints manifold. This leads to erroneous results that violate simple underlying physical properties. The question then is: How should additional constraints information be incorporated in a neural network architecture that this physical information is, at least approximately, honoured?
The idea of adding constraints information to a network is essentially a continuation of the ongoing process of con- necting mathematies with machine learning and explicitly adding known information into neural network rather than having the network implicitly learn it (Willard et al.] 202I). Equivariant networks are one example of this, where the symmetry of problem is explicitly built into the neural network (Thomas et al. 2018). Previous work adding con- straints to a neural network includes Stewart & Ermon (2017): Xu et al. (2018) 8) which add an auxiliary regularization term to the loss function; Raissi et al. (2019) physics informed networks shape the output of the neural network to fulfill particular partial differential equation; while Li & Srikumar (2019) incorporates first order logic directly into the neural network.

---

In this paper we introduce a new paradigm in neural network architectures by developing methodologies that allow us to incorporate such additional information on a dynamical system. Motivated by the field of differential algebraic equations (DAH), we study two different approaches for the solution of the problem. The first is the incorporation of the constraints into the network by Lagrange multipliers, and the second employs so called stabilization techniques that aim to penalize outputs that grossly violate the constraints. Both approaches have similar counterparts in the physical simulation world, and in particular in systems of DAEs (Ascher & Petzold (998). Our methodology is designed for residual neural networks, however it can be used and adopted for other architectures as well. We experiment with it on a number of well-known problems, focusing on molecular dynamics (MD) applications (Allen et al. (004), which enables to incorporate variety of constraints often resulting in significant performance improvement.
The rest of this paper is organized as follows: Section 2 describes various ways that constraints can be introduced in neural networks. Section [3 introduces our model problems, and describes relevant constraints for those problems, while Section 4 performs a series of experiments based on those problems using constrained neural networks. The paper is wrapped up in Section 5 with a discussion and conclusions.
2 RESIDUAL NETWORKS AND CONSTRAINTS
We consider the problem of statistical learning where we have the data pairs (x,.y)). i = 1.... n, and we assume that x 2. and y E V. Our goal is to find function f depending parameters @ that satisfies
f(x, 0) y.
(l)
We focus our attention on functions f(:,)) that are mapped by neural networks. Such functions typically contain sequence of multiplications by "learnable" matrices followed by nonlincarities. this work we particularly focus on the continuous form of residual network architectures that reads
Zo K,x,
(2a)
à = o(2, T [0,0,
(2b)
Yc = K.z(1).
(2c)
The matrix K. typically embeds input data x E 2 in a vector in larger space V. Next, the residual network uses m layers with learnable parameter O. Such parameters can be different weights, normalizations and attention parameters. Finally, the larger space is closed by a learnable closing matrix K..
For the class of problems we solve here, we assume that there is a given vector function c(-) such that, c(y.) = c(K,z(1)) 0. The ODE equation 2b with the initial condition equation 2a represent a trajectory, z(r) in the high dimensional space Ø. This trajectory can be projected into the low dimensional, physical space, 3/ setting, (r) (T) which is assumed to be constrained. Therefore, along trajectories z satisfying equation 2c
c c(K.z)= =0. 0.
(3)
Given the output ye trains the network by first discretizing the ODE (typically by the forward Euler method) and then, minimizing a loss function that measures the difference between ye and y. However, even when such parameters are learned, they rarely yield zero loss, especially on the validation set. Therefore, in general, the constraints c(K.z) are not automatically fulfilled and this can yield results that are physically infeasible.
Our goal is modify the architecture given in equation 2 such that the additional information given by equation B is addressed. We next discuss four such approaches that can be used.
2.1 AUXILIARY REGULARIZATION
The simplest method for incorporating constraint information is to add an auxiliary regularization term to the loss function. Let L be the traditional loss. We extend this by adding a regularization term
La = L 7 c(yy)) Tc(y.)
(4)
where the positive parameter 7) determines the strength of the regularization.
2

---

While it is possible to use auxiliary regularization successfully as evidenced by Stewart & Ermon (2017); Xu et al. (2018), we have found that auxiliary regularization by itself did not boost network's s learning significantly, as shown in Section C Furthermore, we have found that the network may still not satisfy the constraints for the validation set, and thus can fail to generalize. Finally, it can be difficult tune the parameter v) to have a meaningful balance between the data fit and the constraint.
2.2 PROJECTING THE FINAL STATE ONTO THE CONSTRAINT MANIFOLD
Another approach that has been used for problems in image processing where the output is bounded, is to project the final state onto the constraint manifold. To this end, the network used as is, obtaining a vector ye- We then look for perturbation dy such that
c(ye + -y) = 0.
This problem is not well-posed as there are many such perturbations. Therefore, we look for the perturbation with minimal norm, which results in the iteration
Yo++ Ye, By, Tc(y,,), j j 0, 1...,...,,AM
(5)
where Yeg = ye is the output of the unconstrained residual neural network, J is the Jacobian, that is Vyc(y), and B is an approximation to (J J)-1.
Since the network in this method observes the constraint only at the end r 1, solving for the projection may require many iterations. Nonetheless, it is important to note that there are no learnable weights in the projection. Thus, by using implicit differentiation it is possible to write an analytic backward function and therefore, there is no need to hold all the states when computing the projection in order to compute the derivative. Using the constraint at the end can be thought of as a "final layer" before the results are sent into the loss function.
Projecting the state at the end has advantage that it is simple and requires minimal changes to any existing code. However, it may have some serious drawbacks. In particular, the dynamics in equation_can lead to states that are very far from satisfying the constraint. In this case, the constraint due to very large changes in the last layer. We therefore next explore two techniques that can be applied order to follow the constraint throughout the network, at least approximately.
be difficult to fulfil and the training can be difficult
2.3 PROJECTING THE STATE ONTO THE CONSTRAINT MANIFOLD THROUGHOUT THE NETWORK
To change the architecture we recall that when constraints are added to an ODE one obtains DAE, To this end another parameter vector function (a Lagrange multiplier) is added to the system and equation 2b is replaced with
? = oz, (()) + K. J(K.z) X,
T T [
(6)
while requiring equation 3 to hold. (See Ascher & Petzold (998) for derivation.)
Consider first using explicit stepping method discretize the system with respect to z, obtaining
+ 2₁: hōq(zz),
where l refers to the k-th discrete layer in the neural network, and ōq is single discretization of the ODE with respect to z around Zg. With forward Euler it has the form
8q(z)) o (Zz, 0x).
However, since forward Euler may have poor stability when nearly imaginary eigenvalues are present, we use the classical Runge-Kutta 4th order (RK4) discretization which can be written as
1 (kl (ZR; 0x) + 2k2(zz, 0x) + 2k3 0x) + ku(ze, 0k))
ŒqRK, (zx)
6
where each k, is the Runge--Kutte stage that requires a single application resnet. Next, project the update onto the constraint by solving for > A
2 AK J (KZzi+1) TX.
(7)
3

---

To eliminate À we substitute the right hand side into the constraint, obtaining
c(Koz*+1) 1) (Ko(z) + hK? J X)) -
with J -(K, We can now approximately solve this equation for ôz - J T X in the same way we project the constraint. Linearizing we obtain
с (Ka(z + hK! X)) = c(K,z)+)) + KJK,K, J = 0, which gives for the Lagrange multiplier in the kth layer the expression
= h h (JKK JJ) J((K
(8)
Inserting this expression into equation we finally obtain
Zk+1 2 - K (JK.K. "
c(Kozi+1)
(9)
which performs one iteration of a linearized projectiono of our state onto the constraint. In practice, we perform several such projections iteratively
*+1 *+1 K2 JJ KJ,K.K -
j+1
j -.,, M,
(10)
where j is the projection index, and 2 = z x and J ((
Based on this, the residual network corresponding equation can now be written
Zo K,X;
(11a)
ÖqRK (zz) (k) 0x) + 2kg (zy: 0x) + 2kg 0x) (2 0x))
(11b)
Z++1 proje (zR ddRR((z)))). 0, ..
(11c)
Ye = Kezn:
(11d)
with the projection operator defined by proje(z)
farg min 1 2 [|02|1 s.l. c(K,(z 8z)) 0)-
(12)
Here we have introduced the projection to make the network inherently obey constraints at every layer.
In summary, have derived a method for approximately projecting the output of each layer onto the constraints. The projections depend on the learnable matrix Kg, and iteratively apply equation I0 until the desired accuracy is reached.
Next, we note that the projection shown in equation 9 uses Newton-like iterations which are very accurate, but can be rather expensive. We have found that for some problems, simple gradient descent iterations
Mx+1 x - K c(K,zit+1)
(13)
are sufficient for our purpose. For other applications we have also applied few conjugate gradient steps to approxi- mately solve the linear system equation 10
The projection methods introduced in Section 2.2 and here both allow obeying the constraints to arbitrary precision,
but they can be rather expensive due the iterative projections. We next introduce faster alternative.
2.4 STABILIZATION BY PENALTY
The last technique we explore is stabilization with respect to the additional information equation 3 Unlike the projec- tion, stabilization does aim to satisfy the constraints exactly but rather do approximately throughout the network. To this end, note that descent towards the constraint can be written as approximation of an ODE of the form
K JK,z) Hc(K.z),
where H is any Symmetric Positive Definite matrix.
Therefore, one way stabilize the system is to augment the original ODE with term that "flows" towards the constraint. The resulting limit ODE becomes
ă = =z 0) - ((Kz)) Hc(K,z).
(14)
Using again an explicit method to discretize the system we obtain a discrete analog that can be used to solve the problem.
In the following we use the RK4 discretization with H = I, which is often sufficient. But for some problems more advanced choices might be needed (Ascher 1997). One common choice for H is (JK.K. JT)-
4

---

3 MODEL PROBLEMS
In this section we discuss two model problems that are used to test the different architectures. The multi-body pendu- lum is simple model problem that has well studied numerical solution (Amold (017). The ubiquitous molecular dynamics problem has received much attention in the literature (Allen et al. 2004, Kadupitiya et al., 2022, Schütt et al. 2018). Both problems are constrained mechanical systems with known equations of motion that can be solved by standard integration techniques. We use numerical simulations these physical problems to generate datasets that allow us determine and evaluate constrained neural networks.
Our goal is to train a neural network to predict future states of mechanical systems given the present state.
3.1 THE MULTI-BONE BODY PENDULUM
Our first mechanical system is 2D multi-body pendulum, as shown in Figure 1 We chose this toy experiment since it has obvious constraints, as well as chaotic motion which makes its motion non-trivial to predict.
Figure I: A multi-body pendulum system with four pendulums.
A multi-body pendulum can be parameterized in Cartesian coordinates by position r E 2 and velocity v E R matrices, where N is the number of pendulums in the system. Since the distance between consecutive mass points is constant we have that at each time point along trajectory
C
rr r.-1l 4, i 1=
(15)
where ro = 0, and 4 is the length the ith pendulum piece. The Jacobian of the constraint is easily calculated and it is a sparse matrix. In practice, only matrix-vector products with the Jacobian (and its transpose) are needed and this can coded efficiently without storing the matrix.
3.2 WATER MOLECULES SIMULATION
For second experiment have created microcanonical ensemble (NVE) simulation of 32 water molecules at temperature of 300 K, approximated with the Lennard-Jones force-field described in Praprotnik & Janežid (2005). using cp2k (Kühne et al. 2020). The physical simulation is done employing step-size of 0.1 Is for 100, 000 steps. The simulation is done without periodic spatial boundaries in order to have system that fundamentally evolves over time. Constraints are commonly added MD simulations in order freeze out high frequency vibrational movement of molecules, which enables the simulations to use larger time-steps in physical space (Allen et al. 2004).
Water molecules are known to be bound in a triangle configuration with 95.7 pm between the hydrogen and oxy- gen atoms and an angle of 104.5 between the hydrogen atoms as illustrated in Figure Each water molecule is constrained separately and with identical constraints
-
(16)
Note that while the pendulum constraint is exact, the water molecule constraint is not exact, and small deviations exist in the data due to viggling/vibration of the water molecules around their bound state.
4 EXPERIMENTS
In Section 2, we discussed in consecutive subsections several ways for introducing auxiliary constraint information into a neural network, namely, auxiliary loss, end constraints, smooth constraints, and penalty method. Here, we next evaluate these methods on the two model problems discussed in the previous section.
5

---

95.7 pm
95.7 pm
104.5
H
H
Figure 2: An illustration of the bound state of water molecule.
4.1 BULTI-BODY PENDULUM
Our first experiment the multi-body pendulum, described Section the start of each experiment we simulate pendulum for sufficient number of timesteps to generate the amount of data samples needed for training, validation and testing the neural networks plus an additional 1000 steps to ensure chaotic motion. The i-th data sample contains the positions and velocities of the system steps and i k. We randomly split the data into training, validation and testing datasets. When training the neural network on the i-th data sample, we use the input position and velocity as initial best guess, x; (Ie, vi), and try predict the future position y; ri+k-
All multi-body simulations are done with a five-pendulum system, a step length of 0.01 8, a pendulum mass of kg, pendulum length of 1 m. and free-fall acceleration constant of 9.82 ms- We implement our neural networks using RK4 stepping between the layers. The neural network is trained with a batch-size of 10. For the constraints, we use maximum of 100 steepest descent projections obey the constraints with early stopping criterion if the max constraint violation is ever below x 10 m. In this experiment we can control how difficult the problem is, by changing Figure shows two snapshots of five-pendulum configuration.
(a) Snapshot 1.
(b) Snapshot 2.
Figure 3: Snapshots of five-body pendulums. The green and blue pendulums follow the red pendulums by t 20 and k = 50 steps, respectively. The arrows indicate velocity. Snapshot I shows the system before it fully enters chaotic motion, which takes about 1000 steps. Snapshot shows the system later, more chaotic state.
We start by determining appropriate values for y and used in the penalty method and the auxiliary regularization method (For more details see Appendix C We then compare the different ways of adding constraints to neural network. Figure shows the results of training a mimetic neural network without and with constraints (Eliasof et al. 2022).
Based these results, we suspect due to the fact that the neural network no longer predicts solutions even remotely close to obeying the constraints. Instead the solutions but leads to instabilities for large constraint violations. In order to counteract this, we want the network to still predict reasonable pendulum positions before constraints applied, and to that effect have two different methods we can apply. One is to also add penalty term to the end- and smooth-constraints which naturally leads the network towards solution that close to obcying the constraints, as described in Section 2.4 The second method add auxiliary regularization term to the loss function as described in Section 2. based on the constraint violation before any projections onto the constraints are applied the neural network.
that end constraints and smooth constraints lead to unstable training, which we
projected onto the constraints, which works well for small constraint violations,

---

10
Mor constrain voletoon (mm
Maain absouut ornrr (m)
30"
10-4
10-
Auxilary loss 12 tnd constraints Smooth constraints
No constraints Penalty y 10
10
Epochs
Epochs
Figure 4: Comparison of the different ways of adding constraints to neural network trained on 100 multi-body pen- dulum samples with = 20. Left: Mean absolute error. Right: Maximum constraint violation. Constraint violations for end/smooth constraints the number of projections/early stopping criteria. Each run has been repeated five times. The solid/dashed lines are the average based on validation//rainingg data, while the shaded regions around each line designate error bounds of one standard deviation.
not shown since these are several orders of magnitudes lower and determined by
Constraints
n4 = 100
14 = 1000
nt = 10000
No constraints Auxiliary loss 1] = 10
35.0
6.40
2.39
35.1
8.52
2.94
1:
Penalty y (ours)
32.6
5.72
2.32
End constraints 7, 1/ 10 (ours)
29.8
6.89
2.93
Smooth constraints 7, = 10 (ours)
27.9
7.35
2.85
No constraints
67.6
10.7
3.98
Auxiliary loss !) = 10
78.0
31.8
27.4
k 50
Penalty 7 10 (ours)
66.3
13.1
5.43
End constraints
10 (ours)
68.8
40.9
48.8
Smooth constraints 7, 77 (ours)
73.2
40.3
52.7
Table 1: Mean absolute егrоr (cm) over a test set of 1000 samples on the multi-body pendulum problem. ne is the number of training samples, while k is the number of steps predicted ahead. Each experiment was repeated three times and the average of those runs are shown.
We test both methods separately with various strengths, and also combination of the two, and find that the penalty method is insufficient to completely stabilize the training. Auxiliary loss on the other hand does stabilize the network if applied with sufficient strength. For both end constraints, and smooth constraints, auxiliary loss = 10 and penalty 7 = 10 gives the best result. We still refer to these combined methods as end constraints or smooth constraints, but specify a value for 7, m. With the modified smooth and end constraints, the initial comparative experiment as shown in Figure 5
find that combination of
We have investigated the effect of training sample size as well as the difficulty of the problem and show the results of this in table I Figure shows comparative example of multi-body pendulum prediction based of neural networks with different constraining schemes.
7

---

10a
No constraints
Auxiiary
Perilly y-10
Smouth corstraini =
Mor constrent voletion (mii
Maann absouute ornor (m)
30"
10-1
10-2
No constraints Penalty y 10
10
Epochs
Epochs
Figure 5: Comparison of the different ways of adding constraints neural network trained on 100 multi-body pen- dulum samples with = 20. Left: Mean absolute error. Right: Maximum constraint violation. Constraint violations for end/smooth constraints not shown since these are several orders of magnitudes lower and determined by the number of projections/early stopping criteria. Each run has been repeated five times, the solid/dashed lines are the average based on validation/training data while the shaded regions around each line designate error bounds of one standard deviation.
(a) No Constraints 21.2
(b) Auxiliary loss 28.5
(c) Penalty 21.0
(d) End constraints
Smooth constraints
11.5
18.5
Figure 6: Comparative snapshots of the different trained networks during test evaluation. The networks are trained with A 20 1000 training samples. The green pendulum is the desired output, while the blue pendulums the predicted output. For auxiliary loss, end constraints and smooth constraints 17 = 10. For penalty, end constraints and smooth constraints 10. The number given below each pendulum is the mean absolute error (cm) for that particular prediction.
4.2 MOLECULAR DYNAMICS SIMULATIONS WATER MOLECULES
Our second experiment is the water molecule simulation described Section 3.2 Similarly to the multi-body pendu- lum experiment, the ith data sample contains positions and velocities of the system i and + k steps. At the start of datasets based on the 100, 000 physical time steps already simulated. When training the neural network the ith data sample, vectors as initial best guess, x; (r,, v:), and try to predict the future position simulation we change from mimetic neural network to an SO3-equivariant neural network using the e3nn framework (Batzner et al. 2022). The neural network is trained with batch-size of 10. For the smooth constraints, we use a maximum of 100 steepest descent projections and an early stopping of 5 10 A.
each experiment we randomly generate training/validation/eesst
the input position and velocity ri-k- Due to the nature of the

---

Constraints
n4
100
T4
1000
T4 10000
No constraints
6.80
6.67
3.79
Auxiliary loss 1) = 10
6.82
6.70
3.28
h: = 50
Penalty 10 (ours)
4.59
4.05
1.52
End constraints 7r 7 = 10 (ours)
3.63
3.62
3.13
Smooth constraints 7, 7) = 10 (ours)
3.55
3.50
3.18
Table 2: Mean absolute error (pm) over a test set 1000 samples on the water molecule problem. ne is the number of training samples, while k is the number of steps predicted ahead. Each experiment was repeated three times and the average of those runs are shown.
We use the modified versions of smooth-constraints and end-constraints with similar hyper parameters as in the previ- ous experiment. Table shows a comparison of the different constraint methods with sizes of training data.
5 DISCUSSION & CONCLUSION
In this work have proposed different ways of incorporating constraint information in neural network and suc- cessfully demonstrated their effectiveness on two separate modelling problems, namely, the multi-body pendulum, where the constraints exact, as well as the more realistic water molecules simulation where the constraints only approximately represented in the data.
Based on the results in Tables I and 2 we see that including constraints generally gives significant inference benefits, Projection constraints are often the best choice when the amount of training data is low, while the penalty method seems be the best choice for large amounts of training data.
The explicit incorporation of constraints information into neural networks does not appear to be as essential a compo- nent success it is for classical MD simulations, likely due to the fact that this information already inherently present in the data that the networks are learning from. Constraints do nonetheless bring significant improvements, especially in cases where there limited training data. Furthermore, for certain fields, such robotic movements or molecular dynamic propagation our work could be valuable, not just because the constraints might lead to better fitting on train/test data, but also because any constraint violation could lead to critical problems in such fields.
Our method uses RK4 stepping in combined with gradient descent projections, which was sufficient for the problems we have investigated thus far. But future studies using other stepping methods as well as Newton projections should be undertaken. Similarly,
set H I in the penalty term defined in Section 2.4] but know that more advanced expressions are commonly used in traditional DAEs. refformance-wisee while the constraint projection methods in their current form do add significant performance booster, they could be optimized significantly. However, the stabi- lization method which only adds simple Jacobian penalty term, can be used with very little in terms of performance degradation, and can lead to significantly improved results as seen in Tables Tand
DATA AVAILABILITY
The code and data sets will be made available upon publication.
REFERENCES
Michael P Allen et al. Introduction to molecular dynamics simulation. Computational soft matter: from synthetic polymers to proteins, 23(1):1-28. 2004.
Martin Arnold. DAE aspects of multibody system dynamics. Surveys in differential- algebraic equations IV. pp. 41-106, 2017.
U.M. Ascher and L.R. Petzold. Computer Methods for Ordinary Differential Equations and Differential-Alebriie Equations. SIAM, Philadelphia, 1998.
Uri M. Ascher. Stabilization of invariants of discretized differential systems. Numerical Algorithms, pp. 1-24, 1997.

---

Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Kornbluth, Nicola Moli- nari, Tess E Smidt, and Boris Kozinsky. E -qu-quiiarr graph neural networks for data-efficient and accurate interatomic potentials. Nature communications, 131)):1-11,20 2022.
Matteo T Degiacomi. Coupling molecular dynamics and deep learning to mine protein conformational space. Struc- lure, (766):0334100 2019.
Moshe Eliasof, Tue Boesen, Eldad Haber, Chen Keasar, and Eran Treister. Mimetic neural networks: unified framework for protein design and folding. Frontiers in Bioinformatics, 2022.
Mario Geiger, Tess Smidt, Alby M., Benjamin Kurt Miller, Wouter Boomsma, Bradley Dice, Kostiantyn Lapchevskyi. Maurice Weiler, Michat Tyszkiewicz. Simon Batzner, Jes Frellsen, Nuri Jung, Sophia Sanbor,, Josh Rackers, and Michael Bailey. e3nn/e3nn: 2021-05-10, May 2021. URL https::/ddoo org/ /
44757844
J C S Kadupitiya, Geoffrey C Fox, and Vikram Jadhao. Solving newton's 's equations motion with large timesteps using recurrent neural networks based operators. Machine Learning: Science and Technology, 3():025002, apr 2022. doi: 10. 0188/2632--1153/a5555
URL http:::/ htt//ii.orgg/1/ org/10. 10l.10880/2222213//2122
Thomas D Kühne, Marcella Iannuzzi, Mauro Del Ben, Vladimir V Rybkin, Patrick Seewald, Frederick Stein, Teodoro Laino, Rustam Z Khaliullin, Ole Schütr, Florian Schiffmann, et al. Cp2k: An electronic structure and molecular dynamics software package-quickstep: Efficient and accurate electronic structure calculations. The Journal of Chemical Physics, 152(19):194103, 2020.
Tao Li and Vivek Srikumar. Augmenting neural networks with first-order logic. In Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, pp. 292-302, Florence, Italy, July 2019. Association for Computational Linguistics. doi: 10. 111111661//19111111
BRLLhttps://eet htnttps://aclanttolloo. logt
Konstantin Lipnikov, Gianmarco Manzini, and Mikhail Shashkov. Mimetic finite difference method. Journal of Computational Physics, 257:1163-12211 2014.
Tharindu Miyanawala and Rajeev Jaiman. An efficient deep learning technique for the navier-stokes equations: Application to unsteady wake flow dynamics. arXiv preprint arXiv:1710.09092 2017.
Matej Praprotnik and Dušanka Janežič. Molecular dynamics integration and molecular vibrational theory. iii. the infrared spectrum of water. The Journal of chemical physics, 122(17):174103, 2005.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Kristof Schütt,, Huziel Sauceda, P-J Kindermans, Alexandre Tkatchenko, and K-R Müller. Schnet-a deep learning architecture for molecules and materials. The Journal of Chemical Physics, 148(24):241722, 2018.
Russell Stewart and Stefano Emnon. Label-free supervision of neural networks with physics and domain knowledge. In Thirty-First AAAI Conference on Artificial Intelligence, 2017.
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field net- works: Rotation-and mannslaiion-qqivariaan neural networks for 3d point clouds. arXiv preprint arXiv:1802.002/9, 2018.
Benjamin W Wah and Minglun Qian. Violation-guided neural-network learning for constrained formulations in time- series predictions. International Journal of Computational Intelligence and Applications, ((4)):333337, 2001. 2001.
Han Wang, Linfeng Zhang, Jiequn Han, and E Weinan. Deepmd-kin -it: A deep learning package for many-body potential energy representation and molecular dynamics. Computer Physics Communications, 028:178-184 2018.
WS Weiglhofer. On a medium constraint arising directly from maxwell's 's equations. Journal of Physics Mathemat- ical and General, (882), 1994.
Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach, and Vipin Kumar. Integrating scientific knowledge with machine learning for engineering and environmental systems. ACM Computing Surveys (CSUR), 2021.
10

---

Jingyi Xu, Zilu Zhang, Tal Friedman, Yitao Liang, and Guy Broeck. A semantic loss function for deep learning with symbolic knowledge. In International conference on machine learning, pp. 5502-5511, PMLR, 2018.
A NETWORK ARCHITECTURE
In this work we used two different neural networks, a a3 rotation-translation equívariant network and non-equivariant mimetic network (Eliasof et al. 2022). Both networks are residual graph convolutional neural networks, with a learnable stepsize which starts out very small, which ensures that the initial network output is very similar to the best guess.
A.1 EQUIVARIANT NETWORK
Our equivariant network is written using the e3nn software (Geiger et al. 2021). which allows for any geometric tensor to be accurately represented, and allows meaningful interactions between any such objects (Thomas et al., 2018). Our network inspired by |Bttznrr et al. (2022). and reaches comparable results to the ones shown their work on MD17. The propagation block for equivariant network is simple equivariant convolutional filter, non-linear activation, and self-interacting tensorproduct. In this work we limit our network to scalars, pseudo-scalars, vectors, and pseudo-vectors. Our equiavariant network has R 2.8M parameters and 8 layers.
A.2 NON-EQUIVARIANT MIMETIC NETWORK
Our non-equivariant network is mimetic network, with a propagational block which computes node averages and gradients well higher order products by moving node information through their edge connections. The informa- tion transfer between nodes and edges is done in a mimetic fashion (Lipnikov et al. 2014). Our mimetic network has 5.8M parameters in 8 layers.
B DETERMINING PENALTY STRENGTH
The penalty stabilization method mentioned in Section 2.4 contains the hyperparameter y, which determines the strength of the penalty. In Figure we set H = I and test various values of in order to find an appropriate value for it. Note that for this method the constraint violation does not simply diminish when 7 is increased.
C TRAINING WITH AUXILIARY REGULARIZATION
The stabilization method mentioned in Section 2.T employs auxiliary regularization, which contains the hyper param- eter that determines the strength of the regularization. With auxiliary regularization it is easy to minimize constraint violation, but does so at the cost of overall learning as seen Figure
11

---

10P
Max constart viclation
Mizzn absoluto aroor
10-4
1=1
y=1
y 50
y 200
2-2
Epochs
Epachs
Figure 7: Comparison of neural networks using penalty stabilization with different strengths All neural networks are trained 100 samples with k 20 on the multi-body pendulum system described in 4.1 Each run has been repeated five times, the solid/dashed lines are the average based on validation/training data while the shaded regions around each line designate error bounds of one standard deviation.
10"
10a
Mox conotraine voletion (m)
Maai Lisouut antrr mm)
10
10-2
1=0
0=0
0=1
20
10
1000
9 1000
Eoocht
Epochs
Figure 8: Comparison of neural networks using auxiliary regularization with different strengths All neural networks
are trained on 100 samples with 20 on the multi-body pendulum system described 4.1 Fach run has been repeated five times, the solid/dashed lines are the average based on validation/training data while the shaded regions around each line designate error bounds one standard deviation.
12


Reviews:
[{"review_id": "hbOAbHmQ4-", "reviewer": "Reviewer_ku9X", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper investigates methods for including known constaints into neural ODEs. It presents new approaches based on projection and ODE augmentation, which are compared and evaluated on two k-step-ahead prediction tasks. The results show some benefits of including such constraints, in particular in the low-data regime.\n", "summary_of_the_review": "The paper considers the interesting and relevant setup of dynamical systems that satisfy known constraints, and aims at improving learned neural networks by including these constraints into the training regime, by using penalties and projections. \n\nUnfortunately, at the current state the paper does not fulfill its claim that the experimental results show \"that including constraints generally gives significant inference benefits\". This is due to (i) difficulties in understanding the \"significance\" of the results, in particular of table 1, and (ii) the lacking evaluation of the auxiliary+penalty combination, without which claims about the relevance of the projection methods is not possible. \n\nI therefore lean towards rejection. ", "strength_and_weaknesses": "## Strengths\n- The setup is well-presented and well-motivated, and the derivation of the multiple approaches follows very clearly from the initial presentation (eqs. 2 and 3).\n- The proposed methods seem conceptually simple and easy to implement, which increases their potential utility for practitioners that use neural ODEs.\n(The \"Clarity, Quality, Novelty and Reproducibility\" section also lists further strengths)\n\n\n## Weaknesses & Questions\nMy main concerns relate to the experimental evaluation:\n- Experiments - combining auxiliary and penalty loss: The experimental findings indicate that for the \"end constraints\" and \"smooth constraints\" methods, additional losses are necessary to get good performance, and that a combination of auxiliary loss and penalty gives the best result. But using _only_ this combination is not experimentally evaluated, which would be, I believe, necessary to make any claims about the performance of the \"end constraints\" and \"smooth constraints\" methods (since the performance could also come already from just the combination of the auxiliary loss and penalty). This is I believe a necessary evaluation to add, to be able to claim that \"the constraint projection methods [...] add a significant performance booster\". \n- Experiments - interpreting the results: I find it hard to evaluate the results presented in table 1. In table 2, I can observe that the error is roughly halved by using the penalty and/or constraints, compared to the two baselines, and I also observe that the scale of error does not change too much between the three setups. But in table 1, errors range from 66 to 2, which makes me question the relevance of the low-data experiments: If all methods perform \"badly\" for $n_t=100$ and $k=50$, as in they all perform close to random, can we really compare their performance and state that some perform better than others? Maybe this could be clarified by stating relative errors (to make the numbers shown more interpretable) and by adding error ranges (and according to Figures 4 and 5 experiments were performed for multiple seeds anyways).\n  Clarifying this point is particularly releavant since the paper claims that \"results in tables 1 and 2 [show] that including constraints generally gives significant inference benefits\".\n\n\n### Minor things\n- Section 2, \"Next, the residual network uses m layers with learnable parameter \\theta\": This sentence is not clear from the context as \"m\" does not appear anywhere. Please clarify.\n- Section 2.3 appeared to be a bit misleading. It performs a specific derivation with Lagrange multipliers which ends up at equation 9; only to then consider projections in general (eq. 11c & 12), which are then done in multiple ways. It might be helpful for clarity to mention projections earlier (which are already established in section 2.2) and instead discuss the choice of projection method more explicitly (since the method of choice seem to be gradient descent steps; this further motivates that deriving equation 9 might not be necessary since eq. 13 is used instead); and possibly compare this choice experimentally.\n- On \"Neural DAEs\": I am not quite convinced that this is an appropriate name for the proposed methodology. \"Neural ODEs\" can be understood as ODEs $y' = f_\\theta(y)$, where $f_\\theta$ is a neural network. A natural way to think of \"neural DAEs\" would then be to consider DAEs, e.g. of the form $y' = f_\\theta(y, z)$, $0 = g_\\theta(y, z)$, and have both $f_\\theta$ and $g_\\theta$ be neural networks; that is, the constaint would also be learned.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Dy8A-ZQ43B", "reviewer": "Reviewer_vpjY", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This study proposes a neural network for differential algebraic equation (DAE). It is an ordinary differential equation (ODE) with algebraic constraints. The proposal is based on projection, and composed of four implementations.\n", "summary_of_the_review": "The motivation is interesting, but the experiments do not adequately answer the motivation. The manuscript is not ready for publication.", "strength_and_weaknesses": "Strength:\n\nThe motivation is clear.\n\nWeaknesses:\n\nWhile four methods are proposed, the difference between them is unclear. They are compared only by using experiments, and the motivation to propose them separately is unclear.\n\nMany previous studies have tackled to satisfy constraints, but none of them is discussed or compared.\n\nThe definition of neural network architecture in (2) is awkward. $\\sigma$ is not defined explicitly. If it is a neural network or some nonlinear function, (2) is equivalent to Neural ODE sandwiched by two projections by $K_0$ and $K_c$. If it is just an activation function, I cannot understand what (2) does mean. Anyway, it is never a residual network. The relationship with previous studies is unclear.\n\nA detailed derivation of the first method is omitted and the reader is directed to the reference. In contrast, the detailed explanation about 4th-order Runge-Kutta method is provided. The organization is not reader-friendly.\n\nAccording to Table 1, in many cases, the case with no constraints demonstrated better performances. So, the contribution of the proposal is unclear. \n\nDetails about experiments are missing, so it is hard to reproduce experiments. Both experiments involve only constraints that come from links between masses (or molecules). This constraint is so-called holonomic constraint. However, in the introduction, the authors mentioned other constraints such as the divergence of the magnetic field. Hence, the experiments are insufficient or the introduction is excessive advertising.\n\nIt appears to violate the format (namely, smaller margins used).", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "4XcpLDG-YJ", "reviewer": "Reviewer_Ypc6", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper presents 4 different methods to explicitly adding auxiliary trajectory information to neural networks for dynamical systems rather than having the network implicitly learn it. Two experiments are conducted to verify that, depending on the amount of the data, these methods give significant inference benefit.", "summary_of_the_review": "In many simulations of dynamical systems, there are natural invariants that we can describe by constraints. Inherited from this fact, the authors present 4 different methods to explicitly adding auxiliary trajectory information to neural networks for dynamical systems rather than having the network implicitly learn it. The methods are: adding a term to the loss function to maintain the constraints, projecting the final state or every immediate state to the constraint manifolds, stabilizing by penalty. Two experiments with the multi-body pendulum and water molecule are derived. Depending on the amount of the data, these methods give significant inference benefits.\n\nHowever, the paper lacks theoretical justifications for the effects of the methods to the dynamical systems and training. The experiments are weak and the authors have not compared the current methods with previously known methods in the literature as mentioned in the Introduction. Therefore, I believe that the authors need to work further on the improvement of the papers to make it suitable for publication.", "strength_and_weaknesses": "*Strengths*\n- The motivation for the constrained neural network inherited from natural invariants of dynamical systems is clear and the experiments are interesting.\n- The methods are derived in detail and easy to understand.\n\n*Weaknesses*\n- The paper lacks empirical theory. Effects of the proposed 4 methods are discussed based on experiments and there is no theorem/proposition analysing any properties of the methods.\n- There is no comparison with any of the previous methods in the literature. Therefore, it is hard to evaluate the novelty of the results.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject