Augmenting Physical Models with Deep Networks for Complex Dynamics Forecasting
Authors: Yuan Yin
Vincent LE GUEN
Jérémie DONA
Emmanuel de Bezenac
Ibrahim Ayed
Nicolas THOME
patrick gallinari
Venue: ICLR 2021 Oral (2021)
Keywords: spatio-temporal forecasting
deep learning
physics
differential equations
hybrid systems

Abstract: Forecasting complex dynamical phenomena in settings where only partial knowledge of their dynamics is available is a prevalent problem across various scientific fields. While purely data-driven approaches are arguably insufficient in this context, standard physical modeling based approaches tend to be over-simplistic, inducing non-negligible errors. In this work, we introduce the APHYNITY framework, a principled approach for augmenting incomplete physical dynamics described by differential equations with deep data-driven models. It consists in decomposing the dynamics into two components: a physical component accounting for the dynamics for which we have some prior knowledge, and a data-driven component accounting for errors of the physical model. The learning problem is carefully formulated such that the physical model explains as much of the data as possible, while the data-driven component only describes information that cannot be captured by the physical model, no more, no less. This not only provides the existence and uniqueness for this decomposition, but also ensures interpretability and benefits generalization. Experiments made on three important use cases, each representative of a different family of phenomena, i.e. reaction-diffusion equations, wave equations and the non-linear damped pendulum, show that APHYNITY can efficiently leverage approximate physical models to accurately forecast the evolution of the system and correctly identify relevant physical parameters.

Published as a conference paper at ICLR 2021
AUGMENTING PHYSICAL MODELS WITH DEEP NET- WORKS FOR COMPLEX DYNAMICS FORECASTING
Yuan Yin'
Vinceet Le Guen 2.3
Jérémie Dona Emmanuel de Bézenac'
"lbrahim Ayed 1.4
Nicolas Thome
Patrick Gallinari
1 Sorbonne Université, CNRS, LIP6, Paris, France
2 Conservatoire National des Arts et Métiers, CEDRIC, Paris, France
3 EDF R&D, Chatou, France
4 Theresis Lab, Thales
5 Criteo AI Lab, Paris, France
ABSTRACT
Forecasting complex dynamical phenomena in settings where only partial knowl- edge of their dynamics is available is a prevalent problem across various scientific fields. While purely data-driven approaches arguably insufficient in this context, standard physical modeling based approaches tend to be over-simplistic, inducing non-negligible errors. In this work, we introduce the APHYNITY framework, a principled approach for augmenting incomplete physical dynamics described by differential equations with deep data-driven models. It consists in decomposing the dynamics into two components: physical component accounting for the dynamics for which we have some prior knowledge, and a data-driven component accounting for errors of the physical model. The learning problem is carefully formulated such that the physical model explains as much of the data as possible, while the data-driven component only describes information that cannot be captured by the physical model, no more, no less. This not only provides the existence and uniqueness for this decomposition, but also ensures interpretability and benefits generalization. Experiments made on three important use cases, each representative of different family of phenomena, i.e. tions and the non-linear damped pendulum, show that APHYNITY can efficiently leverage approximate physical models to accurately forecast the evolution of the system and correctly identify relevant physical parameters.
equations, wave equa-
INTRODUCTION
Modeling and forecasting complex dynamical systems major challenge in domains such as enviromment and climate (Rolnick et al., 2019), health science (Choi et al., 2016), and in many industrial applications (Toubeau et al., 2018). Model Based (MB) approaches typically rely on partial or ordinary differential equations (PDE/ODE) and stem from deep understanding of the underlying physical phenomena. Machine learning (ML) and deep learning methods are more prior agnostic yet have become state-of-the-art for several spatio-temporal prediction tasks (Shi et al., 2015; Wang et al., 2018; Oreshkin et al., 2020; Donà et al., 2020), and connections have been drawn between deep architectures and numerical ODE solvers, e.g. neural ODEs (Chen et al., 2018; Ayed et al., 2019b). However, modeling complex physical dynamics is still beyond the scope of pure ML methods, which often cannot properly extrapolate new conditions as MB approaches do.
Combining the MB and ML paradigms is an emerging trend to develop the interplay between the two paradigms. For example, Brunton et al. (2016); Long et al. (2018b) learn the explicit form of PDEs directly from data, Raissi et al. (2019); Sirignano & Spiliopoulos (2018) use NNs as implicit methods for solving PDEs, Sco al. (2020) learn spatial differences with graph network. Ummenhofer et al. (2020) introduce continuous convolutions for fluid simulations, de Bézenac et al. (2018) learn the
"qqual contribution, authors sorted by reverse alphabetical order.

---

Published as a conference paper at ICLR 2021
Mean Square Error(MSEE=1155 10-1
MSE=7.6 10 Error 12.9..
MSE=1.9 10 Error 30=0.3%
10-
0.5-
0.5
0.0
0.0
0.5
-05
-10-
GT trajectory
trojectory
-15-
GT tray trajectory prediction
-2.0
prodiction
prodiction
10
20
30 (c) Our APHYNITY framework
(a) Data-driven Neural ODE
(b) Simple physical model
Figure I: Predicted dynamics for the damped pendulum vs. ground truth (GT) trajectories d=0/ar²
w/ sin 0 + add/at: 0. We show that in (a) the data-driven approach (Chen al., 2018) fails to properly learn the dynamics due to the lack of training data, while in (b) an ideal pendulum cannot take friction into account. The proposed APHYNITY shown in (c) augments the over-simplified physical model in (b) with data-driven component, APHYNITY improves both forecasting (MSE) and parameter identification (Error To) compared to (b).
velocity field of an ddlvection-diffuiiii system, Greydanus et al. (2019); Chen et al. (2020) enforce conservation laws in the network architecture or in the loss function.
The large majority of aforementioned MB/ML hybrid approaches assume that the physical model adequately describes the observed dynamics. This assumption is, however, commonly violated in practice. This may be due to various factors, e.g. idealized assumptions and difficulty to explain processes from first principles (Gentine al., 2018), computational constraints prescribing a fine grain modeling of the system (Ayed et al., 2019a), unknown external factors, forces and sources which present (Large & Yeager, 2004). In this paper, aim leveraging prior dynamical ODE/PDE knowledge in situations where this physical model incomplete, i.e. unable to represent the whole complexity of observed data. To handle this case, introduce principled learning framework to Augment incomplete Piysical models for ideNtIfying and forecas Ting complex d'naaiies (APHYNITY), The rationale of APHYNITY. illustrated in Figure I on the pendulum problem, is to augment the physical model when- and only when-it falls short.
Designing general method for combining MB and ML approaches still a widely open problem, and clear problem formulation for the latter is lacking (Reichstein et al., 2019). Our contributions towards these goals are the following:
We introduce simple yet principled framework for combining both approaches. We decom- pose the data into physical and a data-driven term such that the data-driven component only models information that cannot captured by the physical model. We provide existence and uniqueness guarantees (Section 3.1) for the decomposition given mild conditions, and show that this formulation ensures interpretability and benefits generalization.
We propose a trajectory-basee training formulation (Section 3.2) along with an adaptive optimization scheme (Section 3.3) enabling end-to-end learning for both physical and deep learning components. This allows APHYNITY to automatically adjust the complexity of the neural network to different approximation levels of the physical model, paving the way to flexible learned hybrid models.
We demonstrate the generality of the approach on three use cases (reaction-diffusion wave equations and the pendulum) representative of different PDE families (parabolic, hyper- bolic), having a wide spectrum of application domains, e.g. acoustics, electromagnetism. chemistry, biology, physics (Section 4). We show that APHYNITY is able to achieve perfor- mances close to complete physical models by augmenting incomplete ones, both in terms of forecasting accuracy and physical parameter identification. Moreover, APHYNITY can also be successfully extended to the partially observable setting (see discussion in Section 5).
2

---

Published as a conference paper at ICLR 2021
2
RELATED WORK
Correction data assimilation Prediction under approximate physical models has been tackled by traditional statistical calibration techniques, which often rely on Bayesian methods (Pernot & Cailliez, 2017). Data assimilation techniques, e.g. the Kalman filter (Kalman, 1960; Becker et al., 2019), 4D-var (Courtier et al., 1994), prediction errors are modeled probabilistically and a correction using observed data is applied after each prediction step. Similar residual correction procedures commonly used in robotics and optimal control (Chen, 2004; Li et al., 2014). However, these sequential (two-stage) procedures prevent the cooperation between prediction and correction. Besides, in model- baased reinforcement learning, model deficiencies are typically handled by considering only short-term rollouts (Janner et al., 2019) or by model predictive control (Nagabandi et al., 201 The leverage model-based prior knowledge by augmenting with neurally
originality of APHYNITY parametrized dynamics. It does so while ensuring optimal cooperation between the prior model and the augmentation.
Augmented physical models Combining physical models with machine learning (gray-box or hybrid modeling) was first explored from the 1990's: Psichogios & Ungar (1992); Thompson & Kramer (1994); Rico-Martinez et al. (1994) use neural networks to predict the unknown parameters of physical models. The challenge of proper MB/ML cooperation was already raised as a limitation of gray-box approaches but not addressed. Moreover these methods were evaluated applications with a residual targeted to the form of the equation. In the last few years, there has been a renewed interest in deep hybrid models bridging data assimilation techniques and machine learning to identify complex PDE parameters using cautiously constrained forward model (Long al., 2018b; de Bézenac et al., 2018), as discussed in introduction. Recently, some approaches have specifically targetted the MB/ML cooperation. HybridNet (Long et al., 201 20) and PhICNet (Saha et al., 2020) both use data-driven networks to learn additive perturbations or source terms to a given PDE. The former considers the favorable context where the perturbations can be accessed, and the latter the special case of additive noise on the input. Wang et al. (2019); Mehta et al. (2020) propose several empirical fusion strategies with deep neural networks but lack theoretical groundings. settings, but with
specific
PhyDNet (Le Guen & Thome, 2020) tackles augmentation in partially-observe specific recurrent architectures dedicated to video prediction. Crucially, all the aforementioned approaches do not address the issues of uniqueness of the decomposition or of proper cooperation for correct parameter identification. Besides, we found experimentally that this vanilla cooperation is inferior to the APHYNITY learning scheme in terms of forecasting and parameter identification performances (see experiments in Section 4.2).
THE APHYNITY MODEL
In the following, we study dynamics driven by equation of the form:
dX, F(X.)
(I)
dt
defined over finite time interval [0, (), where the state X is either vector-valued, i.e. have X₂ E IRd for every f, (pendulum equations in Section 4), or X, is d-dimensional vector field over a spatial domain Ω C 1 with sraction-diffusion and wave equations in Section 4). We suppose that we have access to set of observed trajectories D {X. : [0. + A Vt E [0, T) dx..at F(X,)), where is the set of X values (either Rdd or vector field). In our case, the unknown F has A as domain and we only assume that F E F, with (F, |I |I) a normed vector space.
(2,3), i.e. Xx(x) gd for every (f, E [0, 7'] x Z
3.1 DECOMPOSING DYNAMICS INTO PHYSICAL AND AUGMENTED TERMS
As introduced in Section 1, we consider the common situation where incomplete information is available on the dynamics, under the form of a family of ODEs or PDEs characterized by their temporal evolution F2 Fp F. The APHYNITY framework leverages the knowledge of F, while mitigating the approximations induced by this simplified model through the combination of physical and data-driven components, F being vector space, we write:
F F, + F.
3

---

Published as a conference paper at ICLR 2021
where Fp Fp encodes the incomplete physical knowledge and Fa E F is the data-driven augmenta- tion term complementing Fp. The incomplete physical prior is supposed to belong to a known family, but the physical parameters (e.g. propagation speed for the wave equation) are unknown and need to be estimated from data. Both Pp and Fa parameters are estimated by fitting the trajectories from D.
The decomposition F Fg + Fa is in general not unique. For example, all the dynamics could be captured by the F. component. This decomposition is thus ill-defined, which hampers the interpretability and the extrapolation abilities of the model. In other words, one wants the estimated parameters of F3 to be as close as possible to the true parameter values of the physical model and Fa to play only a complementary role w.r.t Fp, so as to model only the information that cannot be captured by the physical prior. For example, when F E Fp. the data can be fully described by the physical model, and in this case is sensible to desire Fa to nullified; this of central importance in setting where one wishes to identify physical quantities, and for the model to generalize and extrapolate to new conditions. In a more general setting where the physical model is incomplete, the action of Fa on the dynamics, as measured through its norm, should be as small as possible.
This general idea embedded the following optimization problem:
dX, |all subject to VX E D, Vt, dt
min P FEF
= Fa)
(2)
The originality of APHYNITY is to leverage model-based prior knowledge by augmenting it with neurally parametrized dynamics. It does so while ensuring optimal cooperation between the prior model and the augmentation.
A first key question is whether the minimum in Eq. (2) is indeed well-defined, in other words whether there exists indeed decomposition with a minimal norm Fa. The answer actually depends on the geometry of Fp. and is formulated in the following proposition proven in Appendix B:
Proposition 1 (Existence of a minimizing pair). If Fp is proximinal set', there exists decomposi- tion minimizing Eg. (2).
Proximinality is mild condition which, shown through the proof of the proposition, cannot be weakened. It is a property verified by any boundedly compact set. In particular, it is true for closed subsets of finite dimensional spaces. However, if only existence is guaranteed, while forecasts would be expected to be accurate, non-uniqueness of the decomposition would hamper the interpretability of Pp and this would mean that the identified physical parameters are not uniquely determined.
is then natural to ask under which conditions solving problem Eq. (2) leads to a unique decom- position into physical and a data-driven component. The following result provides guarantees on the existence and uniqueness of the decomposition under mild conditions. The proof is given in Appendix B:
Proposition 2 (Uniqueness of the minimizing pair). F, is Chebyshev set', Eq. (2) admits unique minimizer. The Fg in this minimizer pair is the metric projection of the unknown F onto Fr-
The Chebyshev assumption condition is strictly stronger than proximinality but is still quite mild and necessary. Indeed, in practice, many sets of interest are Chebyshev, including all closed convex spaces in strict normed spaces and, if F L F, can be any closed convex set, including all finite dimensional subspaces. In particular, all examples considered in the experiments are Chebyshev sets.
Propositions and 2 provide, under mild conditions, the theoretical guarantees for the APHYNITY formulation to infer the correct MB/ML decomposition, thus enabling both recovering the proper physical parameters and accurate forecasting.
3.2 SOLVING APHYNITY WITH DEEP NEURAL NETWORKS
In the following, both terms of the decomposition are parametrized and are denoted as F and - - Solving APHYNITY then consists in estimating the parameters 6p and 0,. Op are the physical parameters and are typically low-dimensional, e.g. 2 or 3 in our experiments for the considered physical models. For Fa, we need sufficiently expressive models able to optimize over all F; we
A proximinal set is one from which every point of the space has at least one nearest point. A Chebyshev set from which every point of the space has unique nearest point. More details in Appendix A.
4

---

Published as a conference paper at ICLR 2021
thus use deep neural networks, which have shown promising performances for the approximation of differential equations (Raissi et al., 2019; Ayed et al., 2019b).
When learning the parameters of F and = we have access to finite dataset of trajectories discretized with given temporal resolution At: Pirin {(X RAt 1o5 Iisign. Solving Eq. (2) requires estimating the state derivative dX+/at appearing in the constraint term. One solution is to approximate this derivative using e.g. finite differences as in (Brunton et al., 2016; Greydanus al., 2019; Cranmer et al., 2020). This numerical scheme requires high space and time resolutions in the observation space in order to get reliable gradient estimates. Furthermore it is often unstable, leading to explosive numerical errors as discussed in Appendix D. We propose instead to solve Eq. (2) using an integral trajectory-based approach: we compute K from an initial state X using
the current r- + - dynamics, then enforce the constraint RAC.X. = XkAe This leads to our final objective function on (0p, 0a):
F2- subject to Vi, Vk, X(1 8
min
= x
(3)
RAt
kdr
where X Rirt is the approximate solution of the integral
g
+kAt (K + + dX, obtained
by a differentiable ODE solver.
In our setting, where we consider situations for which F only partially describes the physical phenomenon, this coupled MB ML formulation leads different parameter estimates than using the MB formulation alone, as analyzed more thoroughly in Appendix C. Interestingly, our experiments show that using this formulation also leads to a better identification of the physical parameters 9, than
when fitting the simplified physical model F alone (Section 4). With only an incomplete knowledge on the physics, 8p estimator will be biased by the additional dynamics which needs to be fitted the data. Appendix F also confirms that the integral formulation gives better forecasting results and more stable behavior than supervising over finite difference approximations of the derivatives.
3.3 ADAPTIVELY CONSTRAINED OPTIMIZATION
The formulation in Eq. (3) involves constraints which are difficult enforce exactly in practice. We considered a variant of the method of multipliers (Bertsekas, 1996) which uses a sequence of Lagrangian relaxations Lx,(6 (6p, Oa):
Lx, 10 0.) 11- + Aÿ- (iray 0a)
(4)
where (oray(0p. 0.) E I
1x hár X
(i)
This method needs an increasing sequence (Az), such that the successive minima of Lx, converge to so- lution (at least a local one) of the constrained prob- lem Eq. (3). We select (Az)ı by using an iterative strategy: starting from a value An. we iterate, mini- mizing Lx, by gradient descent?, then update Xy with: Aj+1 = Xj + T2Lirag (ftraa(00+1)) (0j+1). where 72 is chosen hyper-parameter and
Algorithm 1: APHYNITY
Initialization: Ao 2 0, > 0, 0;
for epoch = 1 N Nepoche do for iter in 1 Natr do
for batch in 1: B do
1j+1 = 0g
(0p. 9a). This procedure summarized in Algorithm 1. This adaptive iterative procedure allows us to obtain stable and robust results, in reproducible fashion, as shown in the experiments.
V [Ag Ltros (0,) + [|all]
Az+1 = Ayt T2 T2 (0,+1)
4 EXPERIMENTAL VALIDATION
We validate our approach on 3 classes of challenging physical dynamics: reaction-diffusion, wave propagation, and the damped pendulum, representative of various application domains such as chemistry, biology or ecology (for reaction-diffusion) and earth physic, acoustic, electromagnetism
Convergence local minimum isn't at't necessary, few steps often sufficient for successful optimization.
5

---

Published as a conference paper at ICLR 2021
even neuro-biology (for waves equations). The two first dynamics are described by PDEs and thus in practice should be learned from very high-dimensional vectors, discretized from the original compact domain. This makes the learning much more difficult than from the one- one-iieesiional pendulum case. For each problem, we investigate the cooperation between physical models of increasing complexity encoding incomplete knowledge of the dynamics (denoted Incomplete physics in the following) and data-driven models. We show the relevance of APHYNITY (denoted APHYNITY models) both in terms of forecasting accuracy and physical parameter identification.
4.1 EXPERIMENTAL SETTING
We describe the three families of equations studied in the experiments. In all experiments, F = (A) where A is the set of all admissible states for each problem, and the L norm is computed on Drrain by: [|=1 = Ei, (| )ll All considered sets of physical functionals Fp are closed and convex in F and thus are Chebyshev. In order to enable the evaluation on both prediction and parameter identification, all our experiments are conducted on simulated datasets with known model parameters. Each dataset has been simulated using an appropriate nigh-precision integration scheme for the corresponding equation. All solver-based models take the first state X, input and predict the remaining time-steps by integrating F through the same differentiable generic and common ODE solver (4th order Runge-Kutta), Implementation details and architectures are given in Appendix E.
Reaction-diffusion equations We consider 2D SizzHugh-Nagumo type model (Klaasen & Troy, 1984). The system is driven by the PDE du aAu + (u, n; k), 8t a and are respectively the diffusion coefficients of u and v, A is the Laplace operator. The local reaction terms are R,, k) = u u3 - k v, R, v) =
de bAv + (u, 2) where
St
v. The state is X (u,v) and is defined over a compact rectangular domain Q with periodic boundary conditions. The considered physical models are: Param PDE (a,, b). with unknown (a, b) diffusion terms and without reaction terms: Jp = {Fr. (u, H (aAu, bAv) | a 2 Amin > 0, 2 buin > 0}; Param PDE (a, b, k), the full PDE with unknown parameters: Fp = {Fp.b. (u,,) + (adu + Ru(u, (u, v; k), dAu + Ry v) amin > 0, 2 buin > 0, k 2 kmin > 0}.
Damped wave equations We investigate the damped-war wave PDE: Rw ca + k du 0 where 8
812
is the damping coefficient. The state is X = (w, dey and we consider a compact spatial do- ø
main Ω with Neumann homogeneous boundary conditions. Note that this damping differs from the pendulum, as its effect is global. Our physical models are: Param PDE (c), without damp- ing term: Fp UP (u, + (v, of Au) c E [e, +00) with e 0}; Param PDE (c, k): F.p = [Fo. (u, v) + (v - ku) c, k E [e, +00) with e > 03.
Damped pendulum The evolution follows the ODE d///ae" + w/ sin A + ad#/at = 0, where (t) is the angle, wo the proper pulsation (To the period) and a the damping coefficient. With state X = (0, dBfat), the ODE is Fyo.a X + (d0/az, -wi sin adf/ar). Our physi- cal models are: Hamiltonian (Greydanus et al., 2019), a conservative approximation, with F. =FF (u, v) + (0,H(u,,) v), - -2 H(u e)) H H (R²)), H (IR2) is the first order Sobolev space. Param ODE (wa), the frictionless pendulum: Fp (1200 wo E (e,+ with > 0} Param ODE (Wi, ot), the full pendulum equation: F,, = (Frlia wo; o [e, +00) with 0}.
Baselines As purely data-driven basclines, we use Neural ODE (Chen et al., 2018) for the three problems and PredRNN++ (Wang et al., 2018, for -nitionndiffusion only) which are competitive models for datasets generated by differential equations and for spatio-temporal data. As MB/ML methods, in the ablations studies (see Appendix F), we compare for all problems, to the vanilla MB/ML cooperation scheme found in (Wang et al., 2019; Mehta et al., 2020). We also show results for True PDE/ODE, which corresponds to the equation for data simulation (which do not lead to zero error due to the difference between simulation and training integration schemes). For the pendulum, we compare to Hamiltonian neural networks (Greydanus et al., 2019; Toth et al., 2020) and to the the deep Galerkin method (DGM, Sirignano& & Spiliopoulos, 2018). See additional details in Appendix E.

---

Published as a conference paper at ICLR 2021
Table 1: Forecasting and identification results on the (a) reaction-diffusion (b) wave equation, and (c) damped pendulum datasets. We for (a) 1 x 10
5 10-3 for (b) 330, & 50 and for (c) To = 6, a 0.2 as true parameters. log MSEs are computed respectively over 25, 25, and 40 predicted time-steps. %Err param. averages the results when several physical parameters are present. For each level of incorporated physical knowledge, equivalent best results according to a Student t-test are shown bold. n/a corresponds to non-applicable cases.
3,6 10
10
Dataset
Method
log MSE
&Er param. |Fa
Data- driven
Neural ODE PredRNN++
.3.6610.02
4..00+00.01
n/a
n/a
Incomplese physics
Param PDE (a,
-1.26=+0000
67.6
n/a
Reaction diffusion
APHYNITY Param PDE (a, b)
-5.10-1:0211
2.3
Param PDE (a, b,
-9.34:10.20
0.17
Complete physics
APHYNITY Param PDE (a, k)
-9.35++0.02
0.096
1.5c-6
True PDE
-8.81+0005
n/a
n/a
APHYNITY True PDE
-9.17-10.02
n/a
1.4e-7
Data-driven
Neural ODE
-2.51+0.229
n/a
n/a
Incumplete physics
Param PDE (c)
0.51
10.4
n/a
APHYNITY Param PDE (c)
.4.64:10.25
0.31
71.
equation
Param PDE (c,
4.68:10.55
1.38
n/a
Complese physics
APHYNITY Param PDE (c,
-.099+00.2
0.70
4.54
True PDE
4.66::10:30"
n/a
APHYNITY True PDE
-5.24
n/a
0.14
Data driven
Neural ODE
-2.84:10.70
n/a
n/a
Hamiltonian
-0.35110.10
n/a
APHYNITY Hamiltonian
13.974:1.20
n/a
623
Incomplete physies
Param ODE (we)
-0.14--0.10
13.2
n/a
Damped
pendulum
Deep Galerkin Method (we)
-3.1000.44
22.1
APIΣNITY Param ODE (we))
-7.86.10.60
4.0
Param ODE (wo,
-8.28.10.40
0.45
n/a
Deep Galerkin Method (wo, a))
-3.14.10.40
Complese
APHYNITY Param ODE (wo, a))
-8.31-10.30
0.39
8.5
physics
True ÖDE
-8.58.10.20
n/a
n/a
APHYNITY True ODE
-8.44:10.20
2.3
4.2 RESULTS
We analyze and discuss below the results obtained for the three kind of dynamics. We successively examine different evaluation quality criteria. The conclusions are consistent for the three problems, which allows us to highlight clear trends for all of them.
Forecasting accuracy The data-driven models do not perform well compared to True PDE/ODE (all values are test errors expressed as log MSE): -4.6 for PredRNN+ vs. -9.17 for reaction-diffusine -2.51 vs. -5.24 for wave equation, and -2.84 -8.44 for the pendulum in Table 1. The Deep Galerkin method for the pendulum in complete physics DGM (wo, a), being constrained by the equation, outperforms Neural ODE but is far inferior to APHYNITY models. In the incomplete physics case, DGM (We) fails to compensate for the missing information. The incomplete physical models, Param PDE (a,b) for the reaction- difeiiinniffusion. Param PDE (c) for the wave equation, and Param ODE (wo) and Hamiltonian models for the damped pendulum, have even poorer performances than purely data-driven ones, as can be expected since they ignore important dynamical components, e.g. friction in the pendulum case. Using APHYNITY with these imperfect physical models greatly improves forecasting accuracy in all cases, significantly outperforming purely data-driven models, and reaching results often close the accuracy of the true ODE, when APHYNITY and the true ODE models are integrated with the same numerical scheme (which is different from the one used for data generation, hence the non-null errors even for the true equations), e.g. -5.92 -5244 for wave equation in
This integration scheme is then different from the one used for data generation, the rationale for this choice heing that when training model does not know how exactly the data has been generated.
7

---

Published as a conference paper at ICLR 2021
t- 1.0
2.0
3.0
4.0
1-10
2.0
3,0
(a) Param PDE (a, b), diffusion-only (b) APHYNITY Param PDE (a, 6)
(c) Ground truth simulation
Figure 2: Comparison of predictions of two components u (top) and e (bottom) of the reaction-
diffusion system. Note that /. = is largely beyond the dataset horizon (L = 2.5).
3
(a) Neural ODE
(b) APHYNITY Param PDE (c)
(c) Ground truth simulation
Figure 3: Comparison between the prediction of APHYNITY when is estimated and Neural ODE for the damped wave equation. Note that + 32, last column for (a, b, c) is already beyond the training time horizon (t + 25), showing the consistency of APHYNITY method.
Table This clearly highlights the capacity of our approach to augment incomplete physical models with learned data-driven component.
Physical parameter estimation Confirming the phenomenon mentioned in the introduction and detailed in Appendix C, incomplete physical models can lead to bad estimates for the relevant physical parameters: an error respectively up to 67.6% and 10.4% for parameters in the reaction-diffusion and wave equations, and an еттоr of more than 13% for parameters for the pendulum in Table 1. APHYNITY is able significantly improve physical parameters identification: 2.3% error for the reaction-diffusion. 0.3% for the wave equation, and 4% for the pendulum. This validates the fact that augmenting simple physical model to compensate its approximations not only beneficial for prediction, but also helps to limit errors for parameter identification when dynamical models do not fit data well. This is crucial for interpretability and explainability of the estimates.
Ablation study We conduct ablation studies to validate the importance of the APHYNITY aug- mentation compared to a naive strategy consisting in learning F Fp + F, without taking care on the quality of the decomposition, as done in (Wang et al., 2019; Mehta et al., 2020). Results shown in Table 1 of Appendix F show a consistent gain of APHYNITY for the three use cases and for phys- ical models: for instance for Param ODE (a, b) in reaction-diffusion, both forecasting performances (log MSE 10 vs. -4.56) and identification parameter (Error= 2.33% vs. 6.39%) improve. Other ablation results are provided in Appendix F showing the relevance of the the trajectory basedddedd approach described in Section 3.2 (vs supervising over finite difference approximations of the derivative F).
Flexibility When applied to complete physical models, APHYNITY does not degrade accuracy, contrary to a vanilla cooperation scheme (see ablations in Appendix F). This is due the least action principle of our approach: when the physical knowledge is sufficient for properly predicting the observed dynamics, the model learns to ignore the data-driven augmentation. This is shown by the norm of the trained neural net component Fa, which is reported in Table 1 last columm: expected, |Fa| diminishes as the complexity of the corresponding physical model increases, and, relative to incomplete models, the norm becomes very small for complete physical models (for example in the pendulum experiments, we have Fa || = 8.5 for the APHYNITY model to be compared with 132 and 623 for the incomplete models). Thus, we see that the norm of F, is good indication of how imperfect the physical models F, are. It highlights the flexibility of APHYNITY to successfully adapt to very different levels of prior knowledge. Note also that APHYNITY sometimes slightly improves over the true ODE, as it compensates the error introduced by different numerical integration methods for data simulation and training (see Appendix E).
Qualitative visualizations Results in Figure for reaction -iffiiinn show that the incomplete diffusion parametric PDE in Figure 2(a) is unable properly match ground truth simulations: the
8

---

Published as a conference paper at ICLR 2021
behavior of the two components in Figure 2(a) is reduced to simple independent diffusions due to the lack of interaction terms between u and v. By using APHYNITY in Figure 2(b), the correlation between the two components appears together with the formation of Turing patterns, which is very similar to the ground truth. This confirms that Fa can learn the reaction terms and improve prediction quality. In Figure 3, we for the wave equation that the data-driven Neural ODE model fails at approximating du/dee as the forecast horizon increases: it misses crucial details for the second component dw/dt which makes the forecast diverge from the ground truth. APHYNITY incorporates Laplacian term well as the data-driven Fa thus capturing the damping phenomenon and succeeding in maintaining physically sound results for long term forecasts, unlike Neural ODE.
Extension to non-stationary dynamics We provide additional results in Appendix G to tackle datasets where physical parameters the equations vary in each sequence. To this end, we design an encoder able to perform parameter estimation for each sequence. Results show that APHYNITY accommodates well to this setting, with similar trends as those reported this section.
Additional illustrations We give further visual illustrations to demonstrate how the estimation of parameters in incomplete physical models is improved with APHYNITY. For the faction-diffusion equation, we show that the incomplete parametric PDE underestimates both diffusion coefficients. The difference is visually recognizable between the poorly estimated diffusion (Figure 4(a)) and the true one (Figure 4(c)) while APHYNITY gives a fairly good estimation of those diffusion parameters shown in Figure 4(b).
(a) 0.33 10-3 0.94 x 10-3, diffusion estimated with Param PDE (a, b)
a 0.97 10-3,
(e) 10 bb 5.0x
4.75 x 10-3, diffusion estimated true diffusion
with APHYNITY Param PDE (a, 5)
Figure 4: Diffusion predictions using coefficient learned with (a) incomplete physical model Param
PDE (a, b) and (b) APHYNITY-augmentd Param PDE(a, b), compared with the (c) true diffusion
5 CONCLUSION
In this work, we introduce the APHYNITY framework that can efficiently augment approximate physical models with deep data-driven networks, performing similarly to models for which the underlying dynamics are entirely known. We exhibit the superiority of APHYNITY over data-driven, incomplete physics, and -fate-of-the-rt art approaches combining ML and MB methods, both in terms of forecasting and parameter identification on three various classes of physical systems. Besides, APHYNITY is flexible enough to adapt to different approximation levels of prior physical knowledge.
An appealing perspective is the applicability of APHYNITY on partially-observall video prediction. Besides, hope that the APHYNITY framework will open up the way to the design of a wide range of more flexible MB/ML models, e.g. in climate science, robotics or reinforcement learning. In particular, analyzing the theoretical decomposition properties in a partially setting is an important direction for future work.
settings, such as
ACKNOWLEDGEMENTS:
Funding (P. Gallinari), Chaires de recherche et d' enseeggnmmen en intelligence artificielle (Chaires IA), DI 4CClmm project.

---

John R Dormand and Peter J Prince. A family of embedded runge-kutta formulae. Journal of Published as a conference paper at ICLR 2021
REFERENCES
Ibrahim Ayed, Nicolas Cedilnik, Patrick Gallinari, and Maxime Sermesant. Ep-net: Learning cardiac clectrophysiology models for physiology-base
constraints in data-driven predictions. In Yves Coudière, Valéry Ozenne, Edward J. Vigmond, and Nejib Zemzemi (eds.), Functional Imaging and Modeling of the Heart 10th International Conference, FIMH 2019, Bordeaux, France, June 6-8, 2019, Proceedings, volume 1504 of Lecture Notes in Computer Science, pp. 55-5555666 63. Springer, 2019a.
Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, Julien Brajard, and Patrick Gallinari. Learning dynamical systems from partial observations. arXiv preprint
136,, 2019b.
Philipp Becker, Harit Pandya, Gregor Gebhardt, Cheng Zhao, James Taylor, and Gerhard Neu- mann. Recurrent kalman networks: Factorized inference in high-dimensional deep feature spaces. International Conference on Machine Learning (ICML), 2019.
Dimitri Bertsekas. Constrained Optimization and Lagrange Multiplier Methods (Optimization and Neural Computation Series). Athena Scientific, edition, 1996.
Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz. Discovering governing equations from data by sparse identification nonlinear dynamical systems. Proceedings of the National Academy of 2016.
Sciences, 1 113(15):3932--3
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K. Duvenaud. Neural ordinary differential equations. In Advances in neural information processing systems (NeurIPS), pp. 6571-6583, 2018
Wen-Hua Chen. Disturbance observer based control for nonlinear systems. IEEE/ASME transactions mechatronics, 94):706-710, 2004.
Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks. International Conference on Learning Representations (ICLR), 2020.
Edward Choi, Mohammad Taha Bahadori, Jimeng Sun, Joshua Kulas, Andy Schuetz, and Walter Stewart. RETAIN: An interpretable predictive model for healthcare using reverse time attention mechanism. In Advances in Neural Information Processing Systems (NeurIPS), pp. 3504-3512, 2016.
Philippe Courtier, J-N Thépaut, and Anthony Hollingsworth. strategy for operational implemen- tation of 4d-var, using incremental approach. Quarterly Journal of the Royal Meteorological Society, 13761197 1994.
Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. ICLR 2020 Deep Differential Equations Workshop, 2020.
Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical processes: Incorporating prior scientific knowledge. International Conference on Learning Representations (ICIR), 2018.
Jérémie Donå, Jean-Yves Franceschi, Sylvain Lamprier, and Patrick Gallinari. Pde-driven spa- tiotemporal disentanglement. International Conference Learning Representations (ICLR), 2020.
computational and applied mathematics, 6(1):19-26, 1980.
James Fletcher and Warren Moors. Chebyshev sets. Journal of Australian Mathematical Society, 98:161-223 04 2014. doi: 10.1017/S144678877 711167777774000066 14000661.
P. Gentine, M. Pritchard, S. Rasp, G. Reinaudi, and G. Yacalis. Could machine learning break the convection parameterization deadlock? Geophysical Research Letters, 45(11):5722-571,, 2018.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Advances in Neural Information Processing Systems (NeurIPS), pp. 15353-15363, 2019.
10

---

Published as a conference paper at ICLR 2021
Michael Janner, Justin Fu, Marvin Zhang, and Sergey Levine. When to trust your model: Model- based policy optimization. In Advances Neural Information Processing Systems (NewrIPS), pp- 12519-12530, 2019.
Gordon G Johnson. A nonconvex which has the unique nearest point property. Journal of Approximation Theory, 51(4):289 - 332, 1987.
Rudolph Emil Kalman. A new approach to linear filtering and prediction problems. 1960.
Gene A. Klaasen and William C. Troy. Stationary wave solutions of system of reaction-diffusion equations derived from the terhugh-nagumi equations. SIAM Journal on Applied Mathematics, 44(1):96-110, 1984.
doi: 10.1137/0144008.
William Large and Stephen Yeager. Diumal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies, 05 2004.
Vincent Le Guen and Nicolas Thome. Disentangling physical dynamics from unknown factors for unsupervised video prediction. In Computer Vision and Pattern Recognition (CVPR). 2020.
Shihua Li, Jun Yang, Wen-Hua Chen, and Xisong Chen. Disturbance observer-hased control:
methods and applications. CRC press, 2014.
Yun Long, Xueyuan She, and Saibal Mukhopadhyay. Hybridnet: integrating model-based and data-driven learning to predict evolution of dynamical systems. Conference on Robot Learning (CoRL), 2018a.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-Net: Learning PDEs from data. In International Conference on Machine Learning (ICML), 2018b.
Viraj Mehta, Ian Char, Willie Neiswanger, Youngseog Chung, and Jeff Schneider. Neural dynamical systems. ICLR 2020 Deep Differential Equations Workshop. 2020.
Anusha Nagabandi, Gregory Kahn, Ronald S Fearing, and Sergey Levine. Neural network dynam- ies for model-based deep reinforcement learning with model-free fine-tuning. In 2018 IEEE International Conference on Robotics and Automation (ICRA), pp. 7559-7566. IEEE, 2018
Boris N. Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-BEATS: Neural basis expansion analysis for interpretable time series forecasting. International Conference Learning Representations (ICLR), 2020.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, nigh-performance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d' Alché Buc, E. Fox, and R. Garnett (eds.), Advances in Neurul Information Processing Systems 32, pp. 8024-8035. 2019.
Pascal Pernot and Fabien Cailliez. A critical review of statistical calibratiion/prediition models handling data inconsistency and model inadequacy. AIChE Journal, 63(10):4642 4664446624466.,, 2017.
Dimitris C Psichogios and Lyle H Ungar. A hybrid neural network-first principles approach to process
modeling. AICHE Journal, 38(10): (990)):1999955000 1511, 1992.
Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Physics-informed informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 473:686-707, 2019.
Markus Reichstein, Gustau Camps-Vall Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, and Prabhat. Deep learning and process understanding for data-driven Earth system science. Nature, 566:195 2019.
R Rico-Martinez, JS Anderson, and IG Kevrekidis. Continuous-time nonlinear signal processing: neural network based approach for gray box identification. In Proceedings of IEEE Workshop on Neural Networks for Signal Processing, pp. 596-605. IEEE, 1994.
11

---

Published as a conference paper at ICLR 2021
David Rolnick, Priya L Donti, Lynn H Kaack, Kelly Kochanski, Alexandre Lacoste, Kris Sankaran, Andrew Slavin Ross, Nikola Milojevic-Dupon Natasha Jaques, Anna Waldman Waldman-Brown, et al. Tackling climate change with machine learning. In NeurlPS 2019 workshop on Climate Change with Machine Learning, 2019.
Priyabrata Saha, Saurabh Dash, and Saibal Mukhopadhyay. PHICNNet Physics-incorporaaeed convolu- tional recurrent neural networks for modeling dynamical systems. arXiv preprint arXiv:2004.06243, 2020.
Sungyong Seo, Chuizheng Meng, and Yan Liu. Physics-aware difference graph networks for sparsely- observed dynamics. International Conference on Learning Representations (ICLR), 2020.
Xingjian Shi, Zhourong Chen, Hao Wang, Dit- Yeung, Wii-Kin Wong, and Wang-chun Woo. Convolutional LSTM network: A machine learning approach for precipitation nowcasting. In Advances in neural information processing systems (NeurlPS), pp. 802-810, 2015.
Justin Sirignano and Konstantinos Spiliopoulos. Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339-1364 2018.
Michael Thompson and Mark A Kramer. Modeling chemical processes using prior knowledge and neural networks. AICHE Joumal, 40(8))1328-15)) 1340, 1994.
Peter Toth, Danilo Jimenez Rezende, Andrew Jaegle, Sébastien Racanière, Aleksandar Botev, and Irina Higgins. Hamiltonian generative networks. International Conference on Learning Represen- fations (ICIR), 2020.
Jean-Free Sran-Fraccos
Toubeau, Jérémie Bottieau, François Vallée, and Zacharie De Grève. Deep learning- based multivariate probabilistic forecasting for short-term scheduling in power markets. IEEE Transactions on Power Systems, (22):1203-1215, 2018.
Benjamin Ummenhofer, Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian fluid simulation with continuous convolutions. International Conference on Learning Representations (ICLR), 2020.
Qi Wang, Feng Li, Yi Tang, and Yan Xu. Integrating model-driven and data-driven methods for power system frequency stability assessment and control. IEEE Transactions on Power Systems, 34(6):4557-4568, 2019.
Yunbo Wang, Zhifeng Gao, Mingsheng Long, Jianmin Wang, and Philip S. Yu. PredRNN++: Towards a resolution of the deep-in-time dilemma in spatiotemporal predictive learning. In International Conference Machine Learning (ICML), 2018.
12

---

Published as a conference paper at ICLR 2021
A REMINDER ON PROXIMINAL AND CHEBYSHEV SETS
We begin by giving definition of proximinal and Chebyshev sets, taken from (Fletcher & Moors, 2014):
Definition 1. proximinal set of a normed space (E, - II) is subser C C such that every x E E admits at least a nearest point in C.
Definition 2. A Chebyshev set of normed space (E, II) is subset C C such that every x E E admits unique nearest point in C.
Proximinality reduces a compacity condition in finite dimensional spaces. In general, it is weaker one: Boundedly compact sets verify this property for example.
In Euclidean spaces, Chebyshev sets are simply the closed convex subsets. The question of knowing whether it is the case that all Chebyshev sets are closed convex sets in infinite dimensional Hilbert spaces is still an open question. In general, there exists examples of non-convex Chebyshev sets, famous one being presented in (Johnson, 1987) for non-complete inner- inner-product space.
Given the importance of this topic in upproximation theory, finding necessary conditions for a to be Chebyshev and studying the properties of those sets have been the subject of many efforts. Some of those properties are summarized below:
The metric projection on boundedly compact Chebyshev set continuous.
If the norm is strict, every closed convex space, particular any finite dimensional subspace is Chebyshev.
In Hilbert space, every closed convex set is Chebyshev.
B PROOF OF PROPOSITIONS 1 AND 2
We prove the following result which implies both propositions in the article: Proposition The optimization problem:
dX.
min F.FFF
|[aall subject VX E D, Vt,
+ (a)(.)
(5)
di
equivalent metric projection onto Fp-
If Fy is proximinal, Eq. (5) admits a minimizing pair.
F,,K is Chebyshev, Eq. (5) admits a unique minimizing pair which F, is the metric projection.
Proof. The idea is to reconstruct the full functional from the trajectories of D. By definition, A is the of points reached by trajectories in so that:
A (i E Rd EX X, =]
Then let us define a function FD in the following way: For a E A, we can find X. E D and to such that X a. Differentiating X to which possible by definition of D, take:
dX, di
FD
t=to
For any (Fp, Fa) satisfying the constraint in Eq. (5), then have that (F, Fa)(a) dx./dtja (a( for all e A. Conversely, any pair such that (Fp, Fa) E Fp x and Fp + F. FD, verifies the constraint.
Thus we have the equivalence between Eq. (5) and metric projection formulated as:
minimize Fg E Fp
F Fpl
(6)
13

---

Published as a conference paper at ICLR 2021
If Fp is proximinal, the projection problem admits solution which denote Fp. Taking F -D-F- FF. we have that F + F 1r so that (4:s 42) verifies the constraint of Eq. (2). Moreover, there is (Fp. Fa) satisfying the constraint of Eq. (2), we have that Fg + Fg = F2 by what was shown above and |a = I Fyll 2 IFD Ppl by definition F7. This shows that (F7,
is minimal.
Moreover, Fp is a Chebyshev set, by uniqueness of the projection, if P. + FP then [Pal 1 Thus the minimal pair is unique.
C PARAMETER ESTIMATION IN INCOMPLETE PHYSICAL MODELS
Classically, when set Fp C F summarising the most important properties of system available. this gives simplified model of the true dynamics and the adopted problem is then to fit the trajectories using this model as well as possible, solving:
minimize Fp E Fp
L .X)
subject to NgEz, X = and VI, dxp
F (X:)
dt
where L is a discrepancy measure between trajectories. Recall that X is the result trajectory of an ODE solver taking Xo as initial condition. In other words, we try to find function F, which gives trajectories close possible to the ones from the dataset. While estimation of the function becomes easier, there is then a residual part which left unexplained and this can be a non negligible issue in at least two ways:
When F & Fp. the loss is strictly positive at the minimum. This means that reducing the space of functions J, makes us lose in terms of accuracy.
The obtained function F, might not even be the most meaningful function from F, as would try to capture phenomena which are not explainable with functions in Fpp. thus giving the wrong bias the calculated function. For example, if one considering dampened periodic trajectory where only the period can be learned in Fp but not the dampening, the estimated period will account for the dampening and will thus be biased.
This is confirmed in the paper in Section 4: the incomplete physical models augmented with APHYNITY get different and experimentally better physical identification results than the physical models alone.
Let us compare our approach with this one on the linearized damped pendulum to show how estimates of physical parameters can differ. The equation is the following:
d²0
d0
dr + who
dt
We take the same notations as in the article and parametrize the simplified physical models as:
db
F X
-00)
di
where a > corresponds to w.. The corresponding solution for an initial state Xo. which we denote X can then written explicitly as:
g = fg cos Vat
Let us consider damped pendulum solutions X written
0g Obe cost t
which corresponds to:
d0 d/
d8
F: X
2(0 +
di
"This theory, although necessarily in practice when F overfits small datasct.
14

---

Published as a conference paper at ICLR 2021
It is then easy to see that the estimate of with the physical model alone can be obtained by minimizing:
J. le- cost cos v
This expression depends on and thus, depending on the chosen time interval and the way the integral is discretized will almost always give biased estimates, In other words, the estimated value of a will not give us the desired solution + cost.
On the other hand, for given a, in the APHYNITY framework, the residual must be equal to:
F : (0, (a - 2)8 -2 dt
in order to satisfy the fitting constraint. Here corresponds to 1 + w not to wg as in the simplified case. Minimizing its norm, we obtain a = 2 which gives us the desired solution:
0₂ Oue cost coost
with the right period.
D DISCUSSION ON SUPERVISION OVER DERIVATIVES
In order to find the appropriate decomposition (Fp, Fa), trajectory- based error by solving:
minimize
||Fall
F3 E F E
dx, = (F, + Fx)()),
(8)
subject to
Ag E I, Xg 9 and Vt,
dt
VX E D. L(X,X*))=0
In the continuous setting where the data is available at all times 4, this problem is in fact equivalent to following one:
minimize Fp F
EX-D
/ dX,, dt
Fx((
where the supervision is done directly over derivatives, obtained through finiteddifference schemes. This echoes the proof in Section the Appendix where F can be reconstructed from the continuous data.
However, in practice, data only available at discrete times with a certain time resolution. While Eq. (9) is indeed equivalent Eq. (8) in the continuous setting, in the practical discrete one, the way error propagates is not anymore: For Eq. (8) it is controlled over integrated trajectories while for Eq. (9) the supervision is over the approximate derivatives of the trajectories from the dataset. We argue that the trajectory-bant based approach is more flexible and more robust for the following reasons:
Eq. (8), if F, is appropriately parameterized, is possible to perfectly fit the data trajectories at the sampled points.
The use of finite differences schemes to estimate F as is done in Eq. (9) necessarily induces non-zero discretization error.
This discretization error is explosive in terms of divergence from the trajectories.
This last point is quite important, especially when time sampling is sparse (even though we do observe this adverse effect empirically in our experiments with relatively finely time-sampled trajectories). The following gives a heuristical reasoning as to why this is the case. Let F = F + be the function estimated from the sampled points with an enror such that ||e]loo S a. Denoting X the corresponding trajectory generated by F, we then have, for all X E D:
a(x-x) F(X.) - - c(x,)
Vt,
dt
15

---

Published as a conference paper at ICLR 2021
Integrating over [0, T) and using the triangular inequality as well as the mean value inequality, supposing that F has uniformly bounded spatial derivatives:
[ IX
WE E [0, I(X x
|VFl
Xell + al
which, using variant of the Grönwall lemma, gives us the inequality:
Vt E [J, |X X,Il - VF
1)
When a tends to 0, we recover the true trajectories X. However, a is bounded away from by the available temporal resolution, this inequality gives rough estimate of the way X diverges from them, and it can an equality in many cases. This exponential behaviour explains our choice of a trajectory- -basedd optimization.
E IMPLEMENTATION DETAILS
We describe here the three cases studied in the paper for validating APHYNITY. All experiments are implemented with PyTorch (Paszke et al., 2019) and the differentiable ODE solvers with the adjoint method implemented in torchdi f ffeq."
E. REACTION-DIFFUSION EQUATIONS
The system is driven by FitzHugh-Naguma type PDE (Klaasen & Troy, 1984)
du adu R2 (u, v; Øt
de Øt bAv + (u,
where a and b are respectively the diffusion coefficients of u and v, A is the Laplace operator. The (u,
local reaction terms (u, u; k)
The state X = (u, v) is defined over a compact rectangular domain Ω = 1| with periodic boundary conditions. Ω spatially discretized with 32 x 32 2D uniform square mesh grid. The periodic boundary condition is implemented with circular padding around the borders. A is systematically estimated with a x 3 discrete Laplace operator.
Dataset Starting from randomly sampled initial state X E 10, 122332x32 we generate states by integrating the true PDE with fixed a, and a dataset (a 1x 10-3 bb 5x 10-30-3, k 5x 10-3). We firstly simulate high time-resolution (5taim = 0.001) sequences with explicit finite difference method. We then extract states every ô' = 0.1 to construct our low time-resolution datasets.
We set the time of random initial state to f -0.5 and the time horizon to 2.5. 1920 sequences are generated, with 1600 for training/validation and 320 for test. take the state at = 0 as Xo and predict the sequence until the horizon (equivalent to 25 time steps) in all praction-diffusii experiments. Note that the sub-sequence with t < 0 are reserved for the extensive experiments in Appendix G.I.
Neural network architectures Our Fa here a 3-layer convolution network (ConvNet). The two input channels u) and two output ones are (vu dm). The purely data-driven Neural ODE uses such ConvNet as its F. The detailed architecture is provided Table 2. The estimated physical parameters 0p in Fp are simply a trainable vector (a, b) R2 or (a, k) E R
ag
htttps://githin com http://///m/tgccce////eeeeteo
16

---

Published as a conference paper at ICLR 2021
Table 2:
ConvNet architecture in reaction-diffusion and wave equation experiments, used as data- driven derivative operator in APHYNITY and Neural ODE (Chen et al., 2018).
Module
Specification
2D Conv.
3 3 kernel, 2 input channels, 16 output channels, 1 pixel zero padding
2D Batch Norm.
No average tracking
ReLU activation
2D Conv.
kernel, 16 input channels, 16 output channels, pixel zero padding No average tracking
2D Batch Norm.
ReLU activation
2D Conv.
3 kernel, 16 input channels, 2 output channels, 1 pixel zero padding
Optimization hyperparameters We choose to apply the same hyperparameters for all the reaction- diffusion experiments: Niler = A0 = 1, = x 10-3, x 10
E.2
WAVE EQUATIONS
The damped wave equation defined by
82w
dw Ot
2
Aw + k
0
0t
where e is the wave speed and k is the damping coefficient. The state is X (w, 8t )
We consider a compact spatial domain Ω represented as 64 x 64 grid and discretize the Laplacian operator similarly. A is implemented using 5 5 discrete Laplace operator in simulation whereas in the experiment is a 3 x 3 Laplace operator. Null Neumann boundary condition are imposed for generation.
Dataset öt was set to 0.001 to respect Courant number and provide stable integration. The simu- lation was integrated using a 4th order finite difference Runge-Kutt Kutta scheme for 300 steps from an initial Gaussian state, i.e for all sequence at t O, we have:
(a mo1 ty 20
w(x,y, v, 0) = C. x exp
(10)
The amplitude fixed to and (Io. (o) (32, 32) to make the Gaussian curve centered for all sequences. However, a is different for each sequence and uniformly sampled in [10, 100). The same öt was used for train and test. All initial conditions are Gaussian with varying amplitudes. 250 sequences are generated, 200 are used for training while 50 are reserved as test set. In the main paper setting, e 330 and k 50. As with the reaction diffusion case, the algorithm takes as input a state Xz, (w, due
de (4o) and predicts all states from to + öt up to to + 255t.
Neural network architectures The neural network for F& is a 3-layer convolution neural network with the same architecture as in Table 2. For Fz. the parameter(s) to be estimated is either scalar E E vector (c, k) E R?. Similarly, Neural ODE networks are build as presented in Table 2.
Optimization hyperparameters We use the same hyperparameters for the experiments: Niter 3, Ao = 1, 10 .2 x 103.
E.3 DAMPED PENDULUM
We consider the non-linear damped pendulum problem, governed by the ODE
d²0
d0
+ w sin 0
d/2
dt
where 0(1) is the angle, wy coefficient. With the state X = (0, d0) ), the ODE can be written as dX, F(X,) with F X + de
25 the proper pulsation (To being the period) and the damping To
( der - sin 0
(d0
d0 az
17

---

Published as a conference paper at ICLR 2021
Dataset For each train / validation / test split, simulate dataset with 25 trajectories of 40 timesteps (time interval [0, 20), timestep öt = 0.5) with fixed ODE coefficients (Th = 12, a = 0.2) and varying initial conditions. The simulation integrator is Dormand- -Piinnce Runge-Kutte Kutta method of order (4)5 (DOPRIS, Dormand & Prince, 1980). We also add a small amount of white gaussian noise (a 0.01) to the state. Note that pendulum dataset is much challenging than the ideal frictionless pendulum considered in Greydanus et al. (2019).
Neural network architectures We detail in Table 3 the neural architectures used for the damped pendulum experiments. All data-driven augmentations for approximating the mapping X, F(X;) are implemented by multi-layer perceptrons (MLP) with 3 layers of 200 neurons and ReLU activation functions (except the last layer: linear activation). The Hamiltonian (Greydanus et al., 2019: Toth et al.. 2020) is implemented by a MLP that takes the state X, and outputs a scalar estimation of the Hamiltonian H. of the system: the derivative is then computed by an in-graph gradient of H. with ar 0M
respect to the input: F(X:)
0(d0/ de)
Table 3: Neural network architectures for the damped pendulum experiments. n/a corresponds to non-applicable cases.
Method
Physical model
Data-driven model MLP(in=2, units=200, layers=3, out=
Neural ODE
n/a
Hamiltonian
MLP(in=2, units=200, layers=3, out=l) MLP(in=2, units=200, layers=3, out=1)
n/a
APHYNITY Hamiltonian
MLP units=200, layers=3, out
Param ODE (wag)
trainable parameter wb
n/a
APHYNITY Param ODE (wlg)
trainable parameter 40
MLP(in=2, units=200, layers=3, out=2)
Param ODE (wby o:)
2 trainable parameters wh, À
n/a
APHYNITY Param ODE (wo,
2 trainable parameters wp,
MLP units=200, layers=3, out=2)
Optimization hy hyperparameters The hyperparameters the APHYNITY optimization algorithm (Niter, Åo Ti, T2) were ross-validated on the validation set and are shown in Table 4. All models were trained with a maximum number of 5000 steps with early stopping.
Table 4: Hyperparameters of the damped pendulum experiments.
Method
Niter
Ao
APHYNITY Hamiltonian
5
1
0.1
APHYNITY ParamODE (wo)
5
1
10
APHYNITY ParamODE (wo, A)
5
1000
1
100
F ABLATION STUDY
We conduct ablation studies to show the effectiveness of APHYNITY's adaptive optimization and trajectory-based learning scheme.
F.1 ABLATION TO VANILLA MB/ML COOPERATION
In Table 5, we consider the ablation case with the vanilla augmentation scheme found in Le Guen Thome (2020); Wang et al. (2019); Mehta et al. (2020), which does not present any proper decomposition guarantee. We observe that the APHYNITY cooperation scheme outperforms this vanilla scheme all case, both in terms of forecasting performances (e.g. log MSE= -0.35 vs. -3.97 for the Hamiltonian in the pendulum case) and parameter identification (c.g. Err Param=8.4% vs. 2.3 for Param PDE (a, for reaction-dinn -iffuiiodi. confirms the crucial benefits of APHYNITY's principled decomposition scheme.
18

---

Published as a conference paper at ICLR 2021
Table 5: Ablation study comparing APHYNITY to the vanilla augmentation scheme (Wang et al., 2019; Mehta et al., 2020) for the deactionddiffusion equation, wave equation and damped pendulum.
Dataset
Method
log MSE
%Er Param.
| |
Param. PDE (a, b) with vanilla aug.
4.56-+0.52
8.4
(7.5+:144)el
APHYNITY Param. PDE (a, b)
-5.10:10.21
2.3
(6.7:00.4eel
Reaction-
Param. PDE (a, b, k) with vanilla aug.
-8.04:10.033
25.4
(1.5:=0.2)e-2
diffusion
APHYNITY Param. PDE (a, k)
-9.35:+0.02
0.096
(1.5::0.4)e-66
True PDE with vanilla aug.
8.12++0.05
n/a
(6.1=2.3)e-4
APHYNITY True PDE
-9.17=00.02
n/a
(1.4::08))-7
Param PDE (c) with vanilla aug.
-3.90 0.27
0.51
88.66
Wave APHYNITY Param PDE (c)
4.64+0.25
0.31
71.0
equation
Param PDE (c, k) with vanilla aug.
-5.9 0.10
0.71
25.1
APHYNITY Param PDE (c, k)
-6.0910.28
0.70
4.54
Hamiltonian with vanilla aug.
-0.35+0.11
n/a
837=117
APHYNITY Hamiltonian
-3.97+1.22
n/a
623-168
Param ODE (w²) with vanilla aug. APHYNITY Param ODE (wa)
-7.02:11.7
4.5
148=:49
Damped
-7.8640.6
4.0
132:111
pendulum
Param ODE (We, Q:) with vanilla aug.
-7.60+0..6
4.65
35.5
APHYNITY Param ODE (wo, a)
-8.3140.3
0.39
8.5+2.0
Augmented True ODE with vanilla aug.
-8.40-1.2
n/a
3.4
APHYNITY True ODE
-8.44+0.2
n/a
2.3+0.4
F.2 DETAILED ABLATION STUDY
conduct also two other ablations in Table 6:
derivative supervision: in which F, F, is trained with supervision over approximated derivatives on ground truth trajectory, as performed in Greydanus et al. (2019); Cranmer et al. (2020). More precisely, APHYNITY's Luzj is here replaced with Laeriv = =| dxx - F(X,)l as in Eq. (9), where dX,, is approximated by finite differences on X,- dt
non-adaptive optim.: in which we train APHYNITY by minimizing Fall without the adaptive optimization of À shown in Algorithm 1. This case is equivalent to X T2
We highlight the importance algorithm described in paper) compared to non-adpative optimization: for example in the reaction- diffusion case, log MSE= -4.55 vs. -5.10 for Param PDE (a. Finally, when the supervision occurs on the derivative, both forecasting and parameter identification results are systematically lower than with APHYNITY's trajectory based approach: for example, log MSE=-1I.1 vs. -4.64 for Param PDE (c) in the wave equation. It confirms the good properties of the APHYNITY training scheme.
principled adaptive optimization algorithm (APHYNITY
19

---

Published as a conference paper at ICLR 2021
Table 6: Detailed ablation study on supervision and optimization for the reaction-diffusion equation. wave equation and damped pendulum.
Dutaset
Method
log MSE
%Er Puram.
|Fa||
Augmented Param. PDE (a, b) derivative supervision
-4,42:1.225
12.6
(6.8:10.6)e
Augmented Param. PDE (a, b) non-adaptive optim.
-4,55:0..11
7.5
(7.6::1.0)e
APHYNITY Param. PDE (e,
-5.10:0.21
2.3
(6.7:00.4)el
Bostion
Augmenied Puram. PDE (a, b, derivative supervision
-4.90=10.06
11.7
(1.9::.3)e-1
dilscin
Augmented Param. PDE (a, b, non-adaptive optim.
-9.00=0.02
0.21
(5.5+229))-7
APHYNITY Param. PDE (a, b,
-9.35+00.2
0.096
(1.5+00.4)-6
Augmented True PDE derivative supervision Augmented True PDE non-adaptive optim.
-6.03+00001
n/a
(3.1+00.8e-3
-9.01+0.01
n/a
(1.5:
APHYNITY True PDE
99.17:10.02
n/a
(144008))-7
Augmented Param PDE (c) derivative supervision
-1.16410.88
12.1
0.00024
Augmented Param PDE (c) non-adaptive optim.
-2.57=0.21
3.1
43.6
APHYNITY Param PDE (c)
.4.64:1.25
0.31
71.0
Wine equation
Augmented Param PDE (c, derivative supervision
-4.19=00.36
7.2
0.00012
Augmented Param PDE (c, non-adaptive optim.
-4.93+00.51
1.32
0.054
APHYNITY Param PDE (c,
-6.09-+0.28
0.70
4.54
Augmented True PDE derivative supervision
4.42 = 0.33
n/a
6.02e-5
Augmented True PDE non-adaptive optim.
4.97:10.49
n/a
0.23
APHYNITY True PDE
-5.2410.45
n/a
0.14
Augmented Hamiltonian derivative supervision Augmented Hamiltonian -daativee optim.
+0.83+1..3
n/a
642:t:21
-0.49:10.58
n/a
165:1:00
APHYNITY Hamiltonian
-3.97.1122
n/a
623:168
Augmented Param (wip derivative supervision
-1.02:1..04
5.8
136-133
Augmented Param (44p) non-adaptive optim.
4.30+-1.3
4.4
90.4::27
Dassped peodalum
APHYNITY Param ODE (wig)
-7.6640.6
4.0
132+11
Augmented Param ODE (wpp ()) derivative supervision
-2.1100.2
5.0
3.2+1.7
Augmented Param ODE (wp, (r) non-adaptive optim.
-..69+1.3
1.65
4.8+7.7
APHYNITY Param ODE (wa, Oi)
-8.31+..33
0.39
8.542.0
Augmented True ODE derivative supervision Augmented True ODE non-adaptive optim.
2.1410..3
n/a
4.1110.6
-8.3410.4
n/a
1.410.3
APHYNITY True ODE
-8.4410.2
n/a
2.310.4
G ADDITIONAL EXPERIMENTS
G.1 REACTION-DIFFUSII SYSTEMS WITH VARYING DIFFUSION PARAMETERS
We conduct an extensive evaluation on setting with varying diffusion parameters for reaction- diffusion equations. The only varying parameters are diffusion coefficients, i.c. individual and for each sequence. We randomly sample a E x 10-3 2 10-3 and b E [3 10-3 10 is 10-3 across the dataset.
still fixed
In order to estimate a and b for each sequence, we use here ConvNet encoder E to estimate parameters from reserved frames (t < 0). The architecture of the encoder is similar to the one in Table 2 except that F. takes 5 frames (10 channels) as input and F. outputs a vector of estimated (ā. b) after applying sigmoid activation scaled by 1 x 10-2 (to avoid possible divergence). For the baseline Neural ODE, we concatenate a and each sequence as two channels.
In Table 7, we observe that combining data-driven and physical components outperforms the pure data-driven one. When applying APHYNITY to Param PDE (a, the prediction precision is significantly improved (log MSE: -1.32 vs. -4.32) with and respectively reduced from 55.6% and 54.1% to 11.8% and 8... For complete physics cases, the parameter estimations are also improved for Param PDE (a, b, k:) by reducing emror of a and & (resp. 1.55/0.59 vs. 1.29/0.39).
60% of the error of (3.10 1.23) and 10% 20% of the
The extensive results reflect the conclusion shown in the main article: APHYNITY improves the prediction precision and parameter estimation. The same decreasing tendency of FF, is also confirmed.
20

---

Published as a conference paper at ICLR 2021
Table 7: Results of the dataset of reaction-diffue diffusson with varying (a. across the dataset.
5 10-3 is shared
Method
log MSE
%Er
%Err b
%Err k:
Neural ODE (Chen al., 2018)
-3.61
n/a
n/a
n/a
n/a
driven
Incomplese
Param PDE (a,b)
-1.32:1:0.2
55.6
54.1
n/a
n/a (4.6:
physics
APHYNITY Param PDE (a. b)
-4.32:1..32
11.8
18.7
n/a
Param PDE (a, b, k)
0..54
1.55
3.10
0.59
n/a
APHYNITY Param PDE (a, b, k)
-5.72-1..22
1.29
1.23
0.39
(5.9
(5.9944.3)
Complete physics
True PDE
-8.8610.02
n/a
n/a
n/a
n/a
APHYNITY True PDE
08.5
n/a
n/a
n/a
(1.820.6)e-5
G.2 ADDITIONAL RESULTS FOR THE WAVE EQUATION
We conduct an experiment where each sequence is generated with different wave celerity. This dataset challenging because both and the initial conditions vary across the sequences. For each simulated sequence, an initial condition is sampled as described previously, along with a wave celerity e also sampled uniformly in [300, 400j. Finally our initial state integrated with the same Kunge-Kutta scheme. 200 of such sequences are generated for training while 50 are kept for testing.
For this experiment, we also use a ConvNet encoder to estimate the wave speed e from 5 consecutive Dw ). The architecture of the encoder E' is the same in Table 2 but with 10 input reserved states (w, 06
channels. Here also, & is fixed for all sequences and & = 50. The hyper-parameters used these experiments are the same than described in the Section E.2.
The results when multiple wave speeds c are in the dataset are consistent with the one present when only is considered. Indeed, while prediction performances slightly hindered, the parameter estimation remains consistent for both e and k. This extension provides elements attesting for the robustness and adaptability of our method to more complex settings. Finally the purely data-driven Neural-ODE fails to cope with the increasing difficulty.
Table Results for the damped wave equation when considering multiple e sampled uniformly in [300, 400] in the dataset, k is shared across all sequences and
50.
Method
log MSE
EError
EError k
Fa 2
Neural ODE
0.05610.34
n/a n/a
n/a
driven
Param PDE (c)
-1.3
23.9
n/a
n/a
Incamplete physics
APHYNITY Param PDE (c)
-4.11-0.38
3.2
n/a
171
Param PDE kc
0.25
3.54
1.43
n/a
Complete physics
APHYNITY Param PDE kc
-..8
2.41
0.064
3.64
True PDE (c, k)
0..51
n/a
n/a
n/a
APHYNITY True PDE (c, k)
-4.49-0.22
n/a
n/a
0.0005
G.3 DAMPED PENDULUM WITH VARYING PARAMETERS
To extend the experiments conducted in the paper (section 4) with fixed parameters (To - < = 0.2) and varying initial conditions, evaluate APHYNITY much more challenging dataset where we vary both the parameters (To, a) and the initial conditions between trajectories.
We simulate 500/50/50 trajectories for the train/valid/test sets integrated with DOPRI5. For each trajectory, the period To (resp. the damping coefficient a:) are sampled uniformly in the range [3, 10] (resp. [0, 0.5)).
We train models that take the first 20 steps as input and predict the next 20 steps. To account for the varying ODE parameters between sequences, we encoder that estimates the parameters based
21

---

Published as a conference paper at ICLR 2021
on the first 20 timesteps. In practice. units. The output of the encoder fed as additional input to the data-driven augmentation models and to an MLP with final softplus activations to estimate the physical parameters when necessary (wo E 4 for Param ODE (wn). (wb, a) E R2 for Param ODE (wo, o)).
recurrent encoder composed of layer of 128 GRU
In this varying ODE context, also compare to the state-of-the-ar univariate time series forecasting method N-Beats (Oreshkin et al., 2020).
Results shown in Table 9 consistent with those presented in the paper. Pure data-driven models Neural ODE (Chen et al., 2018) and N-Beats (Oreshkin et al., 2020) fail to properly extrapolate the pendulum dynamics. Incomplete physical models (Hamiltonian and ParamODE (wa)) worse since they do not account for friction. Augmenting them with APHYNITY significantly and consistently improves forecasting results and parameter identification.
Table 9: Forecasting and identification results on the damped pendulum dataset with different parameters for each sequence. log MSEs are computed over 20 predicted time-steps. For each level of incorporated physical knowledge, equivalent best results according to a Student t-test are shown in bold. n/a corresponds non-applicable cases.
Method
log MSE
&Error To %Error
Neural ODE (Chen et al., 2018)
-4.35:10.99
n/a n/a
n/a
N NeBeets (Oreshkin et al., 2020)
-..5:=0.55
n/a n/a 5.6++0.6
n/a n/a
Hamiltonian (Greydanus al., 2019)
-1.31+00.4
n/a
n/a
Incompieeee
APHYNITY Hamiltonian
0.772+0.4
n/a 21.5-19
n/a
plysiss
Param ODE (wo)
-2.66+009
n/a 0.4910.1 n/a (5 n/a
................................
APHYNITY Param ODE (wo)
-5.94.10.7
1.0
n/a
Param ODE (Wp,
-5.71=10.4
4.08
152+:129
Cumples
APHYNITY Param ODE (co, a)
-6.22-+0.77
0.26 n/a
62--27 n/a
plasics
True ODE
1..58.1011
APHYNITY True ODE
-8.58.10.1
n/a
n/a 1.
22


Reviews:
[{"review_id": "urtQibqEUUa", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 7, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\nThis paper outlines a method for forecasting and parameter estimation when you have a partial physics model (possibly with unknown parameters) and time series data. This is a hybrid approach where the data-driven (deep learning) approach only learns the parts not accounted for by the physical model. A key feature is being able to decompose the problem in such a way that the data-driven model only models what cannot be captured by the physical model. The parameters of these two models must be fit jointly so that the physical model's parameters are more correct. They prove existence and uniqueness for this decomposition. \n\nStrong points:\nI think this point of making the decomposition unique (and thus better fitting the physical model's parameters) is intuitive and a good one. To my knowledge, other work on hybrid models do not have theory showing this. (But I haven't extensively read the literature.)\n\nThe authors included extensive experiments with multiple versions of each of three test cases, plus several baselines. They also conducted ablation studies for three features of their method: \n- using a decomposition to ensure that the data-driven model only learns what the physical model cannot contribute\n- the trajectory (integral) approach instead of training on estimated derivatives\n- using adaptive weights in the objective function\n\nWeak points:\nI think that the claim \"We show that APHYNITY can perform similarly to complete physical models by augmenting incomplete ones, both in forecasting accuracy and physical parameter identification.\" is overstated. I compared APHYNITY for incomplete physics to the Param PDE or True PDE with complete physics in Table 1, and this claim seems true for the wave equation only.\n\nI think that the related work section should acknowledge that people have been using neural networks for hybrids with physical models since the 1990s. Example: [1]. This is sometimes called \"gray box\" modeling or \"hybrid\" modeling. Since this section only mentions very recent work, I'm not very confident in the paper's claim that they are the first to have a principled decomposition with existence & uniqueness.\n\nI think that \"Propositions 1 and 2 provide, under mild conditions, the theoretical guarantees for the APHYNITY formulation to infer the correct MB/ML decomposition, thus enabling both recovering the proper physical parameters and accurate forecasting.\" is overstated. The theory is about the decomposition existing & being unique, not about this method being able to find that decomposition. \n1. The theory makes no assumptions about the form of F_a (the data-driven component). Showing that the decomposition is possible is not the same as showing that there exists a neural network that could represent F_a. The different versions of the universal approximation theorem for neural networks that I've seen all require that the neural network is approximating something reasonably well-behaved (such as continuous). I could imagine that the parts of a dataset that cannot be modeled by an ODE/PDE might also be non-continuous, etc.\n2. Further, even if the universal approximation theorem applied here, that would show that such a neural network exists, not how well any existing training methods can find that neural network.\n\nI think \"as expected, ||F_a||^2 diminishes as the complexity of the corresponding physical model increases, and becomes almost null for complete physical models.\" is too strong. It does decrease, but I wouldn't call 0.14 and 2.3 \"almost null.\" \n\nI personally can't vouch for the theory (in Section 3.1 & Appendix A & B), as I'm having a hard time following it. I hope that another reviewer can. Some things that would help me:\n- The intuition would be easier if it's mentioned earlier what \\mathcal{F} is in this paper's examples. L^2? It would be nice to put it before Section 3.1, right when the notation is introduced. Similarly, it would help if an example of \\mathcal{F}_p was given early on. \n- Some examples are given in Section 3.1 & the appendix of conditions that make a set Chebyshev. However, since I was having a hard time knowing which sets & spaces were important in this paper, I was having a hard time understanding which conditions were important. \n\nSince reviewers aren't required to read the appendix, and it's quite long, I read it less closely. However, it contains the proofs and backs up many of the nice claims in the paper. This means that quite a bit of content in the paper is less vetted. \n\n\n\nOther clarification questions:\n\nI may have missed it, but does the paper describe how parameters are fit for the PDEs (physical models)? Example: ParamPDE(a,b) for reaction-diffusion. \n\nDoes the requirement of \\mathcal{F}_p being Chebyshev mean that the parameters of the ODE/PDE need to be closed on one end? In other words, we couldn't search for a parameter in (-infinity, infinity)? \n\nAre the results in Table 1 test errors? (as opposed to training or validation?)\n\nIt's not clear to me what you mean by the \"partially observable\" setting.\n\nFor the ablation studies, we're presumably comparing these variants to the full method's results, as described in Table 1, right? However, there are some discrepancies between Table 1 and Tables 5 & 6. Some of the numbers don't match and one of the wave equation cases is not included in Table 5. The \"Augmented True PDE derivative supervision\" case is also missing for the wave equation in Table 6. \n \n\nMinor points:\n\nIt seems that the definition of the set of observed trajectories (in the beginning of Section 3) is not right. \n- Your observations are for a finite set of values t in [0,T].\n- Your observations are not across all solutions X for the dynamics: just at a finite number of IC/BCs. \n- For the reaction-diffusion & wave equations, X is also a function of spatial x, and the data is collected at finitely many spatial points. \n\nI think instead of defining a set \\mathcal{A} on page 3, it would be clearer to use the space \\mathbb{R}^d? I don't think F maps into just the set A, where A is defined as the set of values that X takes. Couldn't the derivatives have values that F doesn't reach?\n\nOn the bottom of Page 3, the question of the minimum being well-defined is explained as the existence question. However, I'm used to \"well-defined\" in various mathematical contexts meaning roughly \"a unique answer.\" Is there a different meaning of \"well-defined\" in this context?  \n\nSummary: I tentatively recommend accepting this paper. Their method has quite a few empirical results showing improvements. People are interested in combining physical models with machine learning in a variety of scientific application areas, so I could see this being a well-used method. However, I have listed some points above that I would like fixed or clarified. I also hope that someone else can vouch for the proof.  \n\n\n[1] Rico-Martinez, R., J. S. Anderson, and I. G. Kevrekidis. \"Continuous-time nonlinear signal processing: a neural network based approach for gray box identification.\" In Proceedings of IEEE Workshop on Neural Networks for Signal Processing, pp. 596-605. IEEE, 1994.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "kLGd-T6of1", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 9, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "I really enjoyed this paper. It's one of my favorite papers from 2020. Authors: thank you for writing it!\n\nThe paper proposes a framework for jointly fitting the parameters of a physical model (such as the parameters of an ODE or PDE) and learning a neural network to model the error or residual of this physical model. The idea is to find a physical dynamics model $F_p \\in \\mathcal{F}_p$ (e.g. where $\\mathcal{F}_p$ is a set of PDEs with different parameter values) and a neural residual dynamics model $F_a \\in \\mathcal{F}_a$ (e.g. where $\\mathcal{F}_a$ is a hypothesis class of neural networks) which minimize the norm of $F_a$ while constraining the composed dynamics $F = F_a + F_p$ to agree with observed data. Interestingly and importantly, the paper proves the minimum-norm decomposition of the observed dynamics into physical model dynamics and neural residual dynamics is unique, given a condition on the geometry of $F_p$. This condition is that $\\mathcal{F}_p$ should be a Cheybshev set: a sufficient condition is that it is a closed convex set in a strict normed space. To me, this condition seems very mild.\n\nThe paper goes on to show that this method, termed APHYNITY, produces significant gains in predictive accuracy over purely learned methods, purely physics-driven methods, and other forms of combining physical models and learned models. The gain is most significant when the physical model is incomplete. When it is complete there may be some small advantage due to the neural residual accounting for some discretization error; regardless there is no harm, and the learned neural residual is very small. \n\nThe paper also shows that APHYNITY produces better parameter estimates for the physical model in the presence of incomplete physics. This seems important. Our physics models are always approximations and in many interesting applications (climate/atmosphere, bioengineering, mechanics of materials, etc) commonly used physics models may have an interestingly-sized gap with reality. When the identified parameters, not the predictions, are needed for some downstream task such as decision making, APHYNITY could help with better parameter ID.\n\nI thought this paper was clear and well written. The main paper presents an easy-to-follow story, the appendices contain plenty of detail, and when while reading the main paper I wanted more detail on specific points, it was usually easy to follow links to the correct appendices. (Should be true of all papers, but often isn't). \n\n---\nThoughts and feedback: \n- It would be useful to see some simple visual demonstration of the effect of APHYNITY on parameter estimation in incomplete models. (I should note the paper is already quite long and thorough, though).\n- The integral trajectory-based approach (section 3.2, 3.3, motivated in Appendix D) used to fit the parameters of both NN and physics model seems like the \"right\" way to do this to me. Nonetheless it would be interesting to see numerical comparison with the alternative (supervision over derivatives).\n- I wonder if the integral trajectory-based approach, vs the supervision-over-derivatives approach, can be related to traditional methods for parameter ID in ODE?\n- It seems to me that the authors do a good job explaining and relating to prior work on learning physical systems. However, with both the density of recent literature in this space and the decades of work combining ODE solvers and function approximators, if there are missing references (papers on similar work that the current submission does not cite), it's quite likely I would not have noticed.\n- It would be interesting to know if this could help figure out *in what way* a physical model is misspecified, and guide design of an approximate physical model that better captures reality. I suspect yes, although it might simply boil down to this method having a better estimate of the residual than if one just does a least-squares fit with the physical model.\n\n---\nTypos:\n- \"bayesian\" -> \"Bayesian\" (end page 2)\n- \"si\" -> \"is\" (page 3)\n\n---\nI think this is a well written paper likely to be of interest to a large number of ICLR attendees, with some important novel contributions, and that the method proposed has a good chance of being widely adopted in the subfield of ML+physical simulation. ", "extra_scores": {}, "extra_text": {}}, {"review_id": "foEvNOmCUso", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 8, "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": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "### Summary of my understanding\n\nThe authors propose a method of function fitting for differential equations. They premise that a model $F$ is the additive combination of $F_p$ and $F_a$, which denote physics and augmenting parts, respectively. The functional form of the physics part, $F_p$, is given in accordance with prior knowledge, whereas the augmenting part, $F_a$, is modeled by neural nets. The proposed method follows the principle of least action of $F_a$, and the authors suggest solving a constrained optimization problem via a method of Lagrange multipliers. They show numerical results on three PDE/ODE-governed systems.\n\n### Evaluation\n\nI really enjoyed reading the paper, which is well written. The motivation is clearly presented. The related work can be more detailed given the recent active studies on physics + ML but seems sufficient from the viewpoint of ODE/PDE fitting. The proposed method is simple yet reasonable, and the experiments are enough supportive to see the superiority of the method. The ablation study in Appendix F is also very interesting. Possible improvement, which the authors would be aware of, can be found in the lack of experiments on real-world datasets. (I understand the difficulty of finding illustrative real-world examples in this kind of problem, but I couldn't stop pointing out it.) Overall, I think this is certainly sound work.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors propose a method for modeling dynamical systems that balances theoretically derived models, which may be grounded in domain knowledge but subject to overly strict assumptions, with neural networks that can pick up the slack. All reviewers were enthusiastic about this work, appreciating its balance of mathematical rigor and experimental assessment. One concern was that this paper follows on decades of related work, which was difficult to adequately summarize. However, changes made throughout discussion phase did address these concerns.
Decision: Accept (Oral)