Constrained Physical-Statistics Models for Dynamical System Identification and Prediction
Authors: Jérémie DONA
Marie Déchelle
patrick gallinari
Marina Levy
Venue: ICLR 2022 Poster (2022)
Keywords: Deep Learning
Hybrid Models
Differential Equations

Abstract: Modeling dynamical systems combining prior physical knowledge and machine learning (ML) is promising in scientific problems when the underlying processes are not fully understood, e.g. when the dynamics is partially known. A common practice to identify the respective parameters of the physical and ML components is to formulate the problem as supervised learning on observed trajectories. However, this formulation leads to an infinite number of possible decompositions. To solve this ill-posedness, we reformulate the learning problem by introducing an upper bound on the prediction error of a physical-statistical model. This allows us to control the contribution of both the physical and statistical components to the overall prediction. This framework generalizes several existing hybrid schemes proposed in the literature. We provide theoretical guarantees on the well-posedness of our formulation along with a proof of convergence in a simple affine setting. For more complex dynamics, we validate our framework experimentally.

Published as a conference paper at ICLR 2022
CONSTRAINED PHYSICAL-STATISTIE MODELS FOR DYNAMICAL SYSTEM IDENTIFICATION AND PREDICTION
Jérémie Donå"1 11 Marie Déchelle" Marrnn Levy², Patrick Gallinari's
1 Sorbonne Université, CNRS, ISIR, F-75005 Paris, France
Sorbonne Université, CNRS, LOCEAN-IPSL, F-75005 Paris, France
*Criteo Labs, Paris, France fi firstname. 1ast nasinmmiisirr upmc. fr except marina. iissl........ fr
ABSTRACT
Modeling dynamical systems combining prior physical knowledge and machine learning (ML) is promising in scientific problems when the underlying processes are not fully understood, e.g. when the dynamics is partially known. A com- mon practice identify the respective parameters of the physical and ML ML compo- nents is to formulate the problem as supervised learning on observed trajectories. However, this formulation leads to an infinite number possible decompositions. To solve this ill- ill-posesssss, we reformulate the learning problem by introducing an upper bound on the prediction eиror of a physical-statistical model. This al- lows us to control the contribution of both the physical and statistical compo- nents to the overall prediction. This framework generalizes several existing hybrid schemes proposed in the literature. We provide theoretical guarantees the well- posedness of our formulation along with a proof of convergence in a simple affine setting. For more complex dynamics, we validate framework experimentally.
INTRODUCTION
Dynamical systems prediction and identification find crucial applications ranging from medicine and the study of tumors (Hanahan & Weinberg, 2011; Lu & Fei, 2014) oceanic and climate fore- casting (Oreskes et al., 1994; Caers, 2011). The modeling of such systems traditionally rely on ordinary or partial differential equations (ODE/PDE) (Madec, 2008; Marti et al., 2010), and their resolution via numerical solvers and data assimilation (Ghil & Malanotte-Rizzoli, 1991). In real world applications, two main pitfalls occur: first the dynamics may only be partially known and thus do not fully represent the studied phenomena (Rio & Santoleri. 2018); second, the system state may only be partially observed as in ocean models (Gaultier et al., 2013). Machine learning (ML) has be- come complementary approach to traditional physics based models (denoted MB for model based) (Reichstein et al., 2019; Dueben & Bauer, 2018). Both offer advantages: whereas MB approaches generalize and extrapolate better, ML high expressivity approaches benefit from the ongoing growth of available data such as satellite observations, with reduced costs compared to data assimilation.
In that perspective, recent lines of work tackle the learning of hybrid models relying on prior physi- cal knowledge and machine learning (Yin et al., 2021; Mehta et al., 2020). Efficiently learning such decompositions actually means solving two different tasks: system identification, i.e. estimating the parameters of the physical model, and prediction, i.e. recovering the trajectories associated to the dynamics. Both are essential for hybrid MB/ML models of dynamical systems. Whereas predic- tion aims at robust extrapolation, identification accounts for physical interpretability of the MB/ML model. While solving both problems using model-based formulation admits well-known numeri- cal solutions, for example using the adjoint method (Le Dimet & Talagrand, 1986; Courtier et al., 1994), the combination of physical models and deep learning is still an open area of research. In this context, ML applications mainly focus on the prediction task, the expense of the system identifica- tion: Willard et al. (2020) underlines the lack of generalizability of black-box ML models and their
"Equal contribution

---

Published as a conference paper at ICLR 2022
inability to produce physically sound results. Indeed, Ayed et al. (2020) show that without any prior knowledge, the recovered estimates of a dynamical system states are not physically plausible despite accurate predictions. Moreover, as noted by Yin et al. (2021), learning a linear MB/ML decompo- sition with the sole supervision on the system trajectories is ill-posed and admits an infinite number of decompositions. Such observations highlight the need to incorporate physically motivated straints in the learning of bybrid models, c.g. through regularization penalties, and several works propose additional constraints guide the model towards physical solutions (Jia et al., 2019; Yin et al., 2021; Linial et al., 2021). Finally, to complete prior dynamical knowledge with data-driven component and ensure interpretability of the decomposition, we work out principled framework that generalizes previous attempts in the regularization of hybrid models. Our contributions are:
In section 3.1, we introduce a novel way to recover well-posednness and interpretability in the learning of hybrid MB/ML models via the control of upper bound. We extend our framework to incorporate auxiliary data when available to handle complex real-world data.
section 3.2, we propose a novel alternate-ottimization algorithm to learn hybrid models.
In section 3.3, provide analysis of the convergence of the proposed algorithm simplified case and experimentally evidence the soundness of our approach on more com- plex settings of increasing difficulty including challenging real world problems (section 4).
BACKGROUND AND PROBLEM SETUP
We consider dynamical system with state time i denoted Z. Z(t). Z, might be fully only partially observed: we write Z, = (X,, Yi), where X, is the observed component and Y, the unobserved one. The evolution of Z is governed by differential equation with dynamics:
dz,
X1 Y,
(fx
(r
(1)
dt
dt
The objective is to predict trajectories of X, i.e. to model the evolution of the observable part fol- lowing dx. (Z,). For simplicity, we omit the index X in fx and write f(.) A fx).
Dynamical Hypothesis
We assume partial knowledge of the dynamics of the observed X₃:
dX, dt
5(Z)= fZ,) fu,
where f& E Hl is known operator with unknown parameters 0*, and fa H. is the unknown residual dynamics. Ht and H. denote function spaces, see discussion in appendix B.
Learning Problem Our objective is to approximate f with function h learned from the observed data. According to eq. (2), we assume = 14 + hgg hęg E H4. i.c. belongs to the same hy- pothesis space as fy: has the same parametric form. Its parameters are denoted 0g. Note that (r) fk. h- E Ha is represented by free form functional with parameters 0,, e.g. neural network. The learning problem is to estimate from data the parameters of he so that they match the true physical ones and h, approximate at best the unknown dynamics f. In this regard, an intuitive training objective is to minimize a distance d between h hı + hi and f:
d(h, f) Ez~pz [I(Z)) f(Z)ll2
where Pz is the distribution of the state Z that accounts for varying initial states. Each Z defines training sample. Minimizing eq. (3) with h = h& h- enables to predict accurate trajectories but may have an infinite number solutions and h, may bypass the physical hypothesis hg. Thus, interpretability guaranteed. We now develop our proposition to overcome this ill-posedness.
METHOD
In hybrid modeling, two criteria essentials: 1. identifiability, i.e. the estimated parameters of hg
should correspond to the true physical ones; 2. prediction power, i.c. the statistical component h. should complete hy that h hi h, performs accurate prediction over the system states. To
2

---

Published as a conference paper at ICLR 2022
control the contribution of each term h2 and fuu, work upper bounds out of eq. (3) (section 3.1). We then propose to minimize d(h, f) while constraining the upper bounds, which provide us with well-posed learning framework (section 3.2). Besides, show that several previous works that introduced constrained optimization to solve related problems are specific cases of our formulation (Yin et al., 2021; Jía al., 2019; Linial et al., 2021). Finally, we introduce an alternate optimization algorithm which convergence is shown in section 3.3 for a linear approximation of f.
3.1 STRUCTURAT. CONSTRAINTS FOR DYNAMICAL SYSTEMS
To ensure identifiability, we derive regularizations on h2 and h, flowing from the control of an upper bound of d(h, f). In particular, to minimize d( ft)) would enable us to accurately interpret hy as the true fk. and h. as the residual dynamics fa. However, since we do not access the parameters of fu. computing d f&) is not tractable. We then consider two possible situations. In the first one, the only available information on the physical system is the parametric form of fx (or equivalently of ha), training thus only relies observed trajectories (section 3.1.1). In the second one, we consider available auxiliary information about fk that will be used to minimize the distance between h& and fR (section 3.1.2). While the first setting is the more general, the physical prior it relies on is often insufficient to effectively handle real world situations. The second setting makes use of more informative priors and better corresponds to real cases shown in the experiments (section 4.2).
3.1.1 CONTROLLING THE ML COMPONENT AND THE MB HYPOTHESIS
We propose general approach to constrain the learning of hybrid models when one solely access the functional form of hg- In this case, to make hę accountable in our observed phenomena, a solution is to minimize d(ha, f). Following the triangle inequality we link up both errors d(h, and d(b f) (computations available appendix C.I):
d(h, f) I d(h, h&) + d(hx. f) = d(hu, 0) + d(hx,
(4)
We want the physical-statistien model h = hy + ha to provide high quality forecasts. Minimizing the sole upper bound does not ensure such aim, h, is only penalized through d(h,,, 0) and is not optimized to contribute to predictions. We thus propose to minimize d(h, J) while controlling both dh 0) and d(hx; J). Such control of the upper bound of eq. (4) amounts to balancing the contribution of the ML and the MB components. This will be formally introduced in section 3.2.
Link to the Literature The least action principle on the ML component i.e. constraining d(bu 0) is invoked for a geometric argument in (Yin et al., 2021). and appears as co-product of the intro- duction of d(hg, f) in eq. (4). Optimizing (h,,, f) to match the physical model with observations is investigated in (Forssell & Lindskog, 1997).
The general approach of eq. (4) allows us to perform prediction (via h) and system identification (via h₂) on simple problems (see section 4.1). The learning of real-world complex dynamics, via data-driven hybrid models, often fails at yielding physically sound estimation, as illustrated in section 4.2. This suggests that learning complex dynamics requires additional information. In many real-world cases, auxiliary information available in the form of measurements providing comple- mentary information
fk- Indeed, common issue in physics is to infer an unobserved variable of interest (in our case fk parameters 8*) from indirect or noisy measurements that we refer to as proxy data. For instance, one can access a physical quantity but only at coarse resolution, in (Um et al., 2020; Belbute-Peres et al., 2020) and in the real world example detailed in section 4.2. We show in the next subsection how incorporate such an information in order to approximate d(hk, fu).
3.1.2 MATCHING THE PHYSICAL HYPOTHESES: INTRODUCING AUXILIARY DATA
We here assume one accesses a proxy of fk. denoted f E H4. Our goal is to adapt our framework to incorporate such auxiliary information, bringing the regularization induced by f within the scope of the control of an upper bound. This enables us to extend our proposition towards the solving of real world physical problems, still largely unexplored by the ML community. We have:
d(h, f) S d(h, h2) + d(hk. f(")) + I = d(hu, 0) + d(hu, JK") ) +
(5)
where I is a constant of the problem that cannot be optimized (see appendix C.2). In that context, benefit from auxiliary information providing with coarse estimates of 0*, denoted 0Pr.
3

---

Published as a conference paper at ICLR 2022
such that fx" = hss,,, (T)) 2 ft. To use the available 9rr guide estimation towards the true parameters 6* of fu. a simple solution is to directly enforce the minimization of d(hg, f in the parameter space by minimizing ||0, l9112. where 0g are the parameters of hg. Indeed, because SK and SK" have identical parametric forms (as both belong to the same functional space H&), minimizing |0g- 1 will bring byg closer to f and thus to fe- As above, we propose minimize dih, f) while controlling both d(tt, 0) and d( fk as described in section 3.2.
Link to the Literature In (Linial et al., 2021) Jk" stands for true observations used to constrain a learned latent space, minimizing d(hx, Jk") Jia et al. (2019) uses synthetic data as f to pre-train their model which amounts the control an upper bound, see appendix C.3. Finally, this setting finds an extension, when the model fø" is learned model, for example trained using eq. (4), leading to self-supervision approach described in appendix C.4.
3.2 LEARNING ALGORITHM AND OPTIMIZATION PROBLEM
From the upper bounds, we first recover the well-posedness of the optimization and derive a theo- retical learning scheme (section 3.2.1). We then discuss its practical implementation (section 3.2.2).
3.2.
WELL-POSEDNESS AND ALTERNATE OPTIMIZATION ALGORITHM
Recovering Well-Posedness We reformulate the ill-posed learning of minky, EHe xH. d(h, f), by instead optimizing d(h, f) while constraining the upper bounds. Let define Sa and S., as
Sk i E Hp ((hx)) p
S, hh E Hu d(hu, S Mu
(6)
where Pk: Ma Pa are two positive scalars and E(h.)) d( f) in the case of section 3.1.1 and ((hx) = d(hg, f in the case of section 3.1.2. Our proposition then amounts to optimizing d(h. over the Minkowski-sum S₁. S.
he + h. h1 E Sk, hn E )
min d(h, f), heS
This constrained optimization setting enables recover the well-posedness of the optimization problem under the relative compactness of the family of function Hr (proof in appendix D.3).
Proposition 1 (Well-poxednees). Under the relative compactness of Sy, eq. (7) finds a solution that writes h h4 + hu St + Sy. Moreover, this solution is unique.
Alternate Optimization Algorithm As the terms in both upper bounds of eqs. (4) and (5) specif- ically address either h&g or hug we isolate losses relative to hą and hy and alternate projections of h2 on S₃ and hn on S,,, as described in Algorithm 1. Said otherwise, we learn h by alternately optimizing hy (hy being fixed) and h, (hg being fixed). practice, we rely on a dual formulation (see section 3.2.2 and the SGD version of Algorithm in Appendix F).
Algorithm Alternate estimation: General Setting
Result: Converged ha and hi Set Ah =0, 0, Ag minnyeHH d(hu, f), tol R while d(h, f) tol do
hg+1 arg min dlhk + h" f); ha+1 arg min d hy,
(8)
hxES2
ES
4- n +++1
end
The convergence of the alternate projections is well studied for the intersection of convex sets or smooth manifolds (von Neumann, 1950; Lewis & Malick, 2008) and has been extended our setting of Minkowski-sum of convex sets (Lange et al., 2019). Because d as defined in eq. (3) is convex, S, and S₁ are convex sets as soon as H4 and H. are convex (Appendix A). Thus, if d(..) f) is strongly convex, eq. (8) finds one and only solution (Boyd et al.. 2004). However, neither the convexity of H., nor of H4 is practically ensured. Nonetheless, we recover the well coll-posedsees of eq. (7) and show the convergence of Algorithm in the simplified case where / is affine function
4

---

Published as a conference paper at ICLR 2022
of X5 (see section 3.3). For complex PDE where convexity may not hold, we validate approach experimentally and we evidence in section 4 that this formulation enables us to recover both an interpretable decomposition h h& h, and improved prediction and identification performances.
3.2.2 PRACTICAL OPTIMIZATION
Equation (6) involves the choice of μak and juu- In practice, we implement the projection algorithm by descending gradients on the parameters of hy and hus with respect to the following losses: (h(h)) = d,, f) + An E(hs)
(9) where Ab, Abg Ah. are positive real values, dynamically increased/decreased during training. In- deed, d(hu, 0) can be interpreted as a stability loss, preventing the neural networks to trump the physical component. On the other hand, d( f) can be interpreted has an initializution loss yield a first estimate of 0g explaining the dynamics.
Cu(hu) = dh, f) + Xn. d(hu, 0)
Yet. f being unknown: dihh, f) is not tractable. To estimate d(hh, f). rely the trajectories associated to the dynamics. We minimize the distance between the ODE flows on and Φy defined Xo:
by h and f, do(ón, 25), over all initial conditions
dolon. 01) Exo [ [ (T)
(10)
We have: (ds, (f) = 0 + d(h, f) O. Definitions of flows for ODE and in depth considera- tion on these distances are available appendix A. The gradients of (o 2y) with respect to the parameters of h& or fn can be either estimated analytically using the adjoint method (Chen et al., 2018) or using explicit solvers, e.g. Rk45, and computing the gradients thanks to the backpropaga- tion, (Onken & Ruthotto, 2020). To compute eq. (10), we rely a temporal sampling of X: our datasets are composed of n sequences of observations of length N, X =Xi ig...... X No+N where each sequence X follows eq. (2) and corresponds to one initial condition Xi. We then sam- ple the space of initial conditions X1, to compute a Monte-Carlo approximation to dø (1hn, 18)). Let ODESolve be the function integrating any arbitrary initial state 2t, up to time with dynamics h, so that Te ODESolve(IIe, h t). The estimate of do(ón, Q5) then writes as:
N
dslhn, 11) ~
1
ΣΣ
aDEsolve(XX, 1,,44, -xx)x
1=1
Note that the way to compute ODEsolve differs across the experiments (see section 4),
3.3 THEORETICAL ANALYSIS FOR A LINEAR APPROXIMATION
We investigate the validity of our proposition when approximating an unknown derivative with an affine function (interpretable first guess approximators). We here consider hę as a linear function. We do not assume any information on f. thus relieving this section from the need of an accurate prior knowledge fk. In this context, we show the convergence of the learning scheme introduced in Algorithm with é == J), demonstrating the validity our framework this simplified setting. For more complex cases, for which theoretical analysis cannot be conducted, our framework validated experimentally in section 4. All proofs of this section are conducted using the distance dg. Let X be the unique solution the initial value problem:
dXx
dt f(x)) with X = Xo
(11)
With hx( AX and hu (x) A the affine approximation of f writes as:
dX, AX, + DA with X Xo
(12)
dt
where A E (R,p(R), DA E R2. We write X the solution to eq. (12) and X the solution to eq. (12) when DA 0. The alternate projection algorithm with the distance de writes as:
Å arg min [
(((())xp(-7]l,6)) dr AA - x x
(13)
to
DA = arg min
[[))() -x- x dr ADll
(14)
D
5

---

Published as a conference paper at ICLR 2022
where AD. AA > O. As the optimization of eq. (13) is not convex on A, the solution existence and uniqueness is not ensured. The well-posednesss w.r.t A can be recovered by instead considering a simple discretization scheme, e.g. X (AX, + D + X, and solving the associated least square regression, which well-posedness is guaranteed, see details in appendix D.2. Such strategy is common practice in system identification. Theoretical considerations on existence and uniqueness of solutions to eqs. (13) and (14) are hard to retrieve. If A is an invertible matrix:
Proposition 2 (Existence and Uniqueness). If Å is invertible, There exists a unique DA, hence unique X solving eq. (14). (proof in appendix D.4)
Finally, formulating Algorithm as a least square problem in an affine setting (see appendix D.5), prove the convergence of the alternate projection algorithm (appendix D.6)
Proposition 3. For ÀD and A sufficiently high, the algorithm that alternaess benween the estimation of A and the estimation of D following eqs. (13) and (14) converges.
4 EXPERIMENTS
We validate Algorithm 1 on datasets of increasing difficulty (see appendix E), where the system state either fully or partially observed (resp. section 4.1 and section 4.2). We longer rely on an affine prior and explicit fug and h. for each dataset. Performances are evaluated via standard metrics: MSE (lower is better) and relative Mean Absolute Error (rMAE, lower is better). We assess the relevance of our proposition based on eqs. (4) and (5), against NeuralODE (Chen et al., 201 10), Aphynity (Yin et al., 2021) and ablation studies. We denote Ours eq. (4) (resp. Ours eq. (5)) the results when == d(hk. f) i.c eq. (4), (resp. f d(h,, f(") i.c. eq. (5)) When d(hg, f) (resp. d(hu, 0)) is not considered in the optimization, we refer to the results as ddh, +) 0) (resp. d(h, d f))- When h is trained by only minimizing the discrepancy between actual and predicted trajectories the results are denoted <'Only d(h,, f), We report between brackets the standard deviation of the metrics over 5 runs and refer to Appendices F and for training information and additional results.
4.1 FULLY OBSERVABLE DYNAMICS
To illustrate the learning scheme induced by eq. (4), we focus on fully observed low dimensional dy- namics: a simple example emerging from Newtonian mechanics and a population dynamics model.
Damped Pendulum (DPL.) Now a standard benchmark for hybrid models, we consider the motion of pendulum of length L damped due viscous friction (Greydanus et al., 2019; Yin al., 2021). Newtonian mechanics provide an ODE describing the evolution of the angle æ of the pendulum:
#
g/ I. sin (x) + k: =
(15)
We suppose access to observations of the system state Z = (x, t). We consider as physical motion hypothesis hy (ac, 0&) Og (i)...). The true pulsation 0" g/ L of the pendulum has be estimated with 0g. The viscous friction term ki remains to be estimated by hup
Population Dynamics (LV) Voltera ODE system models a prey/predator population dy- namies describing the growth of the preys (z) without predators (y), and the extinction of predators without preys, the non linear terms expressing the encounters between both species:
t = 0x0 - Sry, and U = -y + 8xy
(16)
We observe the system state = (x, and set as prior knowledge: ha y) = (0).z. -08y) 8" (a, has to be estimated by 02 (0). 02). ). h accounts for the non linear terms (Bry, dry).
Experimental Setting For both DPL and LV experiments, we consider the following setting: we sample the space of initial conditions building 100/50/50 trajectories for the train, validation and test sets. The sequences share the same parameters: respectively (7. k), for DPL, and (a, 8, 7, 8) for LV. The parameter 0 is set to a neuron (of dimension 1 in the pendulum and 2 for LV) and hy is a 2-layer MLP. Further experimental details are available in appendices 1, E.2 and F.

---

Published as a conference paper at ICLR 2022
Table 1: Experimental Results for PDL and LV data. The presented metric for parameter evaluation is the 1MAE reported in %. Pred. columms report the prediction log MSE on trajectories on test set.
PDL
LV
Model
rMAB(O;, 0*)
Pred. logMSE
rMAB(O,, 0*)
Pred. logMSE
Ours eq. (4)
1.56 (0.009)
-13.7 (0.84)
7.80 (0.011)
-9.28 (0.75)
Only dh,,
9.35 (0.04)
-13.3 (0.65)
24.5 (0.017)
-9.21 (0.91)
d(h, f) + d(hr::
1.82 (0.01)
-13.4 (0.56)
7.91 (0.02)
-9.01 (0.99)
d(h, f) + d(ha,
11.1 (0.03)
-12.9 (0.29)
9.80 (0.098)
-9.45 (0.55)
Aphynity
6.15 (0.009)
-12.2 (0.13)
21.1 (0.016)
-9.89 (0.53)
NeuralODE
-10.1 (0.32)
-9.11 (1.1)
Identification and Prediction Results Table shows that despite accurate trajectory forecasting, the unconstrained setting <OOll d(h, f) fails at estimating the models parameters, showing the need for regular- of the
ization for identification. Constraining the MI. component can be insufficient: for LV data, both Aphynity and d(h, f) d(h,, 0) do accurately es- timate the model parameters. However, the control of d(hx: f), following eq. (4), significantly improves the parameter identification for both datasets. Indeed, in the PDL case, fAg and f are (pseudd)-periodic of the same period, hence the gain in the performances. Fi- nally, our proposition based on eq. (4) is able to identify the parameters of DPL and LV equation with preci- sion of respectively 1.56% and 7.8% beating all consid- ered baselines. Regarding prediction performances, in under-constrained settings ( <<Only d(h, f' in Table 1), h learns to corrects the inaccurate h4- Table 1 and figs. and (appendix G.1) show that our proposition provides more consistent prediction performances. These experiments confirm that the constraints on h4 and la arising from the control of the upper bound of eq. (4) increase interpretability and maintain prediction performances.
1.54
10-
gradient iterations (x2000)
5
10
28
Figure 1: Affine Case Evolution of the MSE between estimated dynamics (Á, D) and the true one (A, D) with the number of gradients steps for linearized DPL.
Throwback to the Affine Case We verify the convergence proved in section 3.3 using the damped pendulum (eq. (15)) linearized the small oscillations regime (see appendix E.1). Making an affine hypothesis following eq. (12), we apply our alternate projection algorithm and optimize A and D alternately using SGD. Figure shows that we are able to accurately estimate and D using our proposition, recovering both the oscillation pulsation and the damping coefficient.
4.2 HIGH DIMENSIONAL DYNAMICS
We now address the learning of transport equations, describing wide range of physical phenomena such as chemical concentration, fluid dynamics or material properties. We evaluate the learning setting induced by eq. (4) and (5) on two physical datasets depicting the evolution of the temperature advected by time-dependent velocity field U and subject forcing S, following:
8T 0t
+ V.(TU S(U)
(17)
The system state Z -T, is partially observed, we only access Every quantities, observed or to estimate, are regularly sampled on spatiotemporal grid: at each timestep t, the time varying velocity field U, writes as U, (u,, v,) and u, u, T, and the forcing term S, are all of size 64 64.
Experimental Setting We consider physical prior the advection i.e (T, 02) -V. (TOk). Thus, 8% is time-dependent, as we learn it to approximate 8* U. We identify the velocity field fx from observations of T, learning mapping between T and U parameterized by neural network Gg, so that 0₂ = (T, ..., T.) x U,, which is common practice in oceanography (Béréziat & Herlin, 2015). G is optimized following eq. (9). S remains be learned by ha- hy implements

---

Published as a conference paper at ICLR 2022
Table 2: Results for Adv+S and Natl data. We report the MSE (x 100) the predicted observations T, the velocity fields and the source term S over time steps on test set.
Adv+S
Natl
Models
T
U
S
T
U
S
Ours eq. (4)
0.74 (0.05)
1.99 (0.13)
0.17 (0.01)
8.27 (0.06)
11.72 (0.07)
6.01 (0.08)
Ours eq. (5)
6.86 (0.12)
6.81 (0.07)
4.35 (0.11)
Aphynity
0.85 (0.35)
3.07 (0.74)
0.18 (0.05)
8.18 (0.16)
11.75 (0.49)
6.02 (0.02)
NeuralODE
1.35 (0.02)
8.83 (0.98)
differentiable semi-Lagrangian scheme (Jaderberg et al., 2015) (see appendix E.3) and h, is ResNet. G. is a UNet. Training details and a schema our model are to be found in appendix F.
Synthetic Advection and Source (Adv+S) To test the applicability of the learning setting in-
v Iruth
duced by eq. (4) on partially observed settings, first study synthetic setting (denoted Adv+S) of eq. (17) by generating velocity fields U, simu- lated following (Boffetta et al., 2001) and adding a source term S inspired by (Frankignoul, 1985). The simulation details are given in appendix E.3.
v essimated
u trutt
colimated
Real Ocean Dynamics (Natl) We consider dataset emulating real world observations the North ATLantic ocean (denoted Natl) (Ajayi et al., 2019). Modeling the evolution of in Natl is challenging its dynamics is chaotic and highly non-linear. This simulation is representative of the complexity encountered in real world data. The principled approach of eq. (4) is insufficient here and one must resort to additional physical informa- tion. We illustrate section 4.2 and make use of aux- iliary data: satellite observations provide a coarse estimate of surface velocity fields (appendix E).
Testimated
truth
Figure Best viewed in color, Estimations of S, and = (u, v) on Adv+S. Prediction ranges from 1 to 20 half-days.
The goal to refine the approximated velocity fields to the ocean dynamics. We proceed as described in eq. (5) and enforce d(hx, (²") supervising G with the proxy data (appendix E.3).
Identification and Prediction Results Table 2 indicates that for Adv+S dataset, we estimate ac- curately the unobserved velocity fields. Qualitatively, Figure shows that controlling our proposed upper bound eq. (4) facilitates the recovery of truthful velocity fields U along with an accurate pre- diction of T. For the highly complex Natl, Table 2 shows that the introduction of auxiliary data following the formulation in eq. (5) significantly helps identification, as the dynamics is too com- plex to able to recover physically interpretable velocity fields using the bound of eq. (4). Regarding prediction performances estimates of U, our model provides more precise prediction than NODE and Aphynity. For real world data, thanks to the proxy data model recovers better velocity fields terms while providing a better estimate for T. Besides, adding prior knowledge in the prediction systems improves pre- diction performances: appendix G shows that NODE minimizes d(h, f) by predicting average and blurred frames. This shows the need for regularization when learning on structured physical data.
the Adv+S data, Table shows that thanks to truthful
Ablation Study We present in Table an ablation study on the Adv+S dataset evidencing the influence of our learning choices on the resolution of both identification and prediction tasks (see appendix G for detailed results). "Joint" rows of Table 3 indicate that the learning of hú and h2 is done simultaneously. Table 3 shows that the sole optimization of d(h, f) fails at estimating phys- ically sounded U. This evidences the ill ool-posedness in such unconstrained optimization. Table indicates that all introduced regularizations improve the recovery of U w.r.t. the <OOnl d(h, f)>
8

---

Published as a conference paper at ICLR 2022
baseline, while adding d(by,, 0) significantly improves both prediction performances and velocity fields estimation. We highlight that the altematt optimization performs better compared to optimiz- ing jointly all parameters of h& and h. Notably, our proposition to optimize h& and h, altemately beats all baselines on both T prediction and U identification (Table 3, Joint rows). Finally, jointly trained models fail estimating U in Table 3, forcing hy to capture the whole dynamics.
Table 3: Ablation Study on Adv4S. We report the MSE (x 100) on the predicted observations T, the velocity fields U and the source S over 6 time steps. "Joint" rows refer the simultancous optim. of hy and hu-
Training
Models
T
U
S
Ours (U known)
0.52
n/a
0.19
Ours eq. (4)
0.74 (0.05)
1.99 (0.13)
0.17 (0.01)
Only d(h, f)
1.02 (0.16)
4.08 (0.23)
0.19 (0.06)
Alternae
d(h, f) + - d(kx, f)
1.02 (0.09)
3.66 (0.15)
0.19 (0.03)
dhh, f) + d(hu: 0)
0.77 (0.06)
2.38 (0.17)
0.19 (0.01)
Ours eq. (4) Only dhh, f)
1.44 (0.08)
3.30 (0.18)
0.30 (0.03)
Joint
1.38 (0.19)
6.96 (0.21)
0.39 (0.08)
5 RELATED WORK
Grey-box or hybrid modeling, combining ODE/PDE and data based models, has received an increas- ing focus in the machine learning community (Rico-Martinez et al., 1994; Thompson & Kramer, 1994; Raissi et al., 2020b). Hybrid approaches allow for alleviated computational costs for fluid simulation (Tompson et al., 2017: De Avila Belbute-PPres. et al., 2020; Wandel et al., 2021), and show better prediction performances through data specific constraints that preserve physics (Raissi al., 2020a; Jia al., 2019). They offer increased interpretability via constraints on convolutional filters (Long et al., 2018; 2019) or on learned residual (Geneva & Zabaras, 2020). Physical knowl- edge, introduced through ODE/PDE regularization (Psichogios & Ungar, 1992: Bongard & Lipson, 2007: de Bézenac et al., 2018) or Hamiltonian priors (Greydanus et al., 2019: Lee et al., 2021), increases generalization power w.r.t pure ML approaches. Closer to our work, (Mehta et al., 2020; San & Maulik, 2018; Young et al., 2017; Saha et al., 2020) study the learning of a physical model augmented with a statistical component. Yin et al. (2021) tackle the same task, ensuring the unique- ness in the decomposition by constraining the norm of the ML component. We generalize latter approaches and address the well-posedness in the learning of hybrid ML/MB models through ad- ditional regularization phenomena relying on hybrid MB/ML models, one major task lies in the estimation of the MB part parameters. This can be done using neural networks (Raissi et al., 2019; Mehta et al., 2020). How- ever, identification tasks being intrinsically ill-posed (Sabatier, 2000), imposing prior knowledge or regularization is necessary to ensure sound estimations (Stewart & Ermon, 2017). Yet, using only prediction as supervision, the recovered parameters are not physically interpretable (de Bézenac al., 2018: Ayed et al., 2020). To ensure uniqueness of the estimation solution, Ardizzone et al. (2018) use invertible neural networks. Linial et al. (2021); Tait & Damoulas (2020): Saemundsson al. (2020) combine variational encoding (Kingma & Welling, 2013) and PDE model, sampling the space of initial conditions and parameters to solve both identification and prediction. However, such methods only deal with low-dimensional dynamics. Besides low dimensional systems, we also show the soundness of our approach on complex high dimensional and partially observed dynamics.
the estimated parameters of the physical part. Indeed, to describe natural
6 DISCUSSION
We propose in this work an algorithm to learn hybrid MB/ML models. For interpretability purposes, we impose constraints flowing from an upper bound of the prediction error and derive a learning algorithm in general setting. We prove its well posedness and its convergence in linear ap- proximation setting. Empirically, we evidence the soundness of our approach thanks to ablation studies and comparison with recent basclines several low and high dimensional datasets. This work can see several extensions: considering non uniform 2 or 3-D grid for climate models, further considerations on the investigated upper bounds, or different decomposition hypothesis.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
We would like to thank all members of the MLIA team from the ISIR laboratory of Sorbonne Uni- versité for helpful discussions and comments. We acknowledge financial support from ANR AI Chairs program via the DLACLIM ANR-19-CHIA- 0018-01 project, the LOCUST ANR project (ANR-15-CF23-0027 and the European Union's 's Horizon 2020 research and innovation program under grant agreement 825619 (AI4EU), The Nat660 data were provided by MEOM research team, from the IGE laboratory from the Université Grenoble Alpes.
REFERENCES
Adekunle Ajayi, Julien Le Sommer, Eric Chassignet, Jean-Marc Molines, Xiaobiao Xu, Aurelie Albert, and Emmanuel Cosme. Spatial and temporal variability of north atlantic eddy field at scale less than 100km. Earth and Space Science Open Archive, pp. 28, 2019. doi: 10.1002) essoar. 10501076.1.
Lynton Ardizzone, Jakob Kruse, Carsten Rother, and Ullrich Köthe. Analyzing inverse problems with invertible neural networks. In International Conference on Learning Representations, 2018.
Ibrahim Ayed, Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari. Learning the spatio- temporal dynamics of physical processes from partial observations. In ICASSP 2020 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3232-3236, 2020.
Filipe de Avila Belbute-Peres, Thomas Economon, and Zico Kolter. Combining differentiable pde solvers and graph neural networks for fluid flow prediction. In International Conference on Ma- chine Learning. pp. 2402-2411 2211. PMLR, 2020.
Dominique Béréziat and Isabelle Herlin. Coupling dynamic equations and satellite images for mod- clling ocean surface circulation. Communications in Computer and Information Science, 550: 191-205, 2015. doi: 10. 00.0007978.33311--5117--1
URL https://hal.inii hnw.aaaal.iiiiaaff////aa frrria.ff/aal-0111 0014h369.
Guido Boffetta, G Lacorata, G Redaelli, and A Vulpiani. Detecting barriers to transport: a review of different techniques. Physica D: Nonlinear Phenomena, 159(1-2):58-70 2001.
Josh Bongard and Hod Lipson. Automated reverse engineering of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 104(24):9943-9948, 2007.
Stephen Boyd, Stephen Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge uni- versity press, 2004.
Jef Caers.
Modeling uncertainty in the earth sciences. John Wiley & Sons, 2011.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary dif- ferential equations. In Samy Bengio, Hanna Wallach, Hugo Larochelle, Kristen Grauman, Nicolò Cesa-Biancha and Roman Garnett (eds.), Advances in Neural Information Processing Systems Curran Associates, Inc., 2018.
31, pp. 6571-6583
PHILIPPE Courtier, J-N Thépaut, and Anthony Hollingsworth. A strategy for operational imple- mentation of 4d-var, using an incremental approach. Quarterly Journal of the Royal Meteorolog- ical Society, 1377519) 1994.
Filipe De Avila Belbute-Peres. Thomas Economon, and Zico Kolter. Combining differentiable PDE solvers and graph neural networks for fluid flow prediction. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 9402-2411. PMLR, 13-18 Jul 2020.
Emmanuel de Bézenac, Arthur Pajot, and Patrick Gallinari. Deep learning for physical processes: Incorporating prior scientific knowledge. In International Conference on Learning Representa- fions, 2018.
10

---

Published as a conference paper at ICLR 2022
Michail Diamantakis. The semi-lagrangian technique in atmospheric modelling: current status and future challenges. In Seminar on Recent Developments in Numerical Methods for Atmosphere and Ocean Modelling, 2-5 September 2013, pp. 183-200, Shinfield Park, Reading, 2014. ECMWF, FCMWE.
Jérôme Droniou. Intégration et Espaces de Sobolev à Valeurs Vectorielles. working paper or www.aarciiieeoouueetssrraaaaaaaa
preprint, April 2001. URL
Peter Dueben and Peter Bauer. Challenges and design choices for global weather and climate models based machine learning. Geoscientific Model Development, 11(10):3999-44) 4009, 2018.
U. Forssell and Lindskog. Combining semi-physical and neural network modeling: An ex-
ample ofits usefulness. IFAC Proceedings Volumes, 10(11):767-77, 1997. ISSN 1474-6670. doi: https://doi.oom http://doi.or//1O.. 100.00001///4444660000000 tion (SYSID' 97), Kitakyushu, Fukuoka, Japan, 8-11 July 1997.
IFAC Symposium on System Identifica-
Claude Frankignoul. Sea surface temperature anomalies, planetary waves, and air-sea feedback in the middle latitudes. Reviews of geophysics, 23(4):357-390, 1985.
Lucile Gaultier, Jacques Verron, Jean-Michel Brankart, Olivier Titaud, and Pierre Brasseur. On the inversion of submesoscale tracer fields to estimate the surface ocean circulation. Journal of Marine Systems, 126:33 42, 2013.
Nicholas Geneva and Nicholas Zabaras. Modeling the dynamics of pde systems with physics- constrained deep auto-regressive networks. Journal of Computational Physics, 403:109056, 2020. ISSN 0021-9991. doi: https://dii.or//1..1000 2ttp://000.11001//...00111111100
Michael Ghil and Paola Malanotte-Rizzoli. Data assimilation in meteorology and oceanography. Advances in geophysics, 33:141-266 1991.
Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks. In Hanna Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' Alché Buc, Emily Fox, and Roman Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 15379- Curran Associates, Inc., 2019.
Douglas Hanahan and Robert A Weinberg. Hallmarks of cancer: the next generation. cell, 144(5): 646-674, 2011.
Phillip Isola, Jan Zhu, Tinghui Zhou, and Alexei A Efros. Image-to-image translation with conditional adversarial networks. CVPR, 2017.
Max Jaderberg, Karen Simonyan, Andrew Zisserman, and koray kavukcuoglu. Spatial transformer networks. In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett (eds.), Advances in Neurul Information Processing Systems, volume 28. Curran Associates, Inc., 2015. URL https: Paper.pdf.
Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan Read, Jacob Zwart, Michael S Steinbach, and Vipin Kumar. Physics guided rnns for modeling dynamical systems: A case study in simulating lake temperature profiles. In SIAM International Conference on Data Mining, SDM 2019, SIAM International Conference on Data Mining, SDM 2019, pp. 558-566. Society for Industrial and Applied Mathematics Publications, 2019. doi: 0..1137719781161197755
Diederik Kingma and Max Welling. Auto-encoding variational bayes. arXiv:1312.6114. 2013.
arXiv preprint
Kenneth Lange, Joong-Ho Won, and Jason Xu. Projection Minkowski sums with application to constrained learning. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 3642-3651. PMLR, 09-15 Jun 2019.
Laançois-Xavier Le Dimet and Olivier Talagrand. Variational algorithms for analysis and assimi- lation of metcorological observations: theoretical aspects. Tellus A: Dynamic Meteorology and Oceanography, 3812)977110,, 1986.
11

---

Published as a conference paper at ICLR 2022
Seungjun Lee, Haesang Yang. and Woojae Seong. Identifying physical law of hamiltonian systems via meta-learning, In International Conference on Learning Representations, 2021. URL https: Sopppmmciiionn//ommmm//////////-
Adrian S. Lewis and Jéröme Malick. Alternating projections on manifolds. Mathematics of Opera- tions Research, 33(1):216-255 234,, 2008. doi: 10. 1128//moor.107.00222
Ori Linial, Neta Ravid, Danny Eytan, and Uri Shalit. Generative ode modeling with known un- knowns. In Proceedings of the Conference on Health, Inference, and Learning, pp. 79-94, 2021.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-Net: Learning PDEs from data. In Jennifer Dy and Andreas Krause (eds.), Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 3208-3216, Stockholmsmässan, Stockholm Sweden, July 2018. PMLR.
Zichao Long. Yiping Lu, and Dong. PDE-Net 2.0: Learning PDEs from data with a numeric- symbolic hybrid deep network. Journal of Compulational Physics, 199:108925, 2019.
Guolan Lu and Baowei Fei. Medical hyperspectral imaging: a review. Journal of biomedical optics, 19(1):010901, 2014.
Gurvan Madec. NEMO ocean engine. Note du Pôle de modélisation, Institut Pierre-Simon Laplace (IPSL), France, No 27, 2008.
Olivier Marti, Pascale Braconnot, J-L Dufresne, Jacques Bellier, Rachid Benshila, Sandrine Bony, Patrick Brockmann, Patricia Cadule, Arnad Caubel, Francis Codron, et al. Key features of the ipsl ocean atmosphere model and its sensitivity atmospheric resolution. Climate Dynamics, 34 (1):1-26, 2010.
Viraj Mehta, Ian Char, Willie Neiswanger, Youngseog Chung, Andrew Oakleigh Nelson, Mark D Boyer, Egemen Kolemen, and Jeff Schneider. Neural dynamical systems: Balancing structure and flexibility in physical prediction. arXiv preprint arvi:22006. /2682, 2020.
Derek Onken and Lars Ruthotto. Discretize-optimizz vs. optimize-discretize for time-series regres- sion and continuous normalizing flows. arXiv preprint arXiv:2005. 2020.
Naomi Oreskes, Kristin Shrader-Frechette, and Kenneth Belitz. Verification, validation, and confir- mation of numerical models in the earth sciences. Science, 263(5147):641-646, 1994.
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 De Vito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, high-performance deep learning library. In Hanna Wallach, Hugo Larochelle, Alina Beygelzimer, Florence d' Alché Buc, Emily Fox, and Roman Gamett (eds.), Advances in Neural Information Processing Systems 32, pp. 8026-8037. Curran Associates, Inc., 2019.
Dimitris C. Psichogios and Lyle H. Ungar. hybrid neural network-first principles approach to process modeling. AIChE Journal, 38(10)::499-1511 1992. doi: http://doi.or//10 1002/aic. 690381003, URL
Maziar Raissi, Paris Perdikaris, and George E Kamiaddkis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Maziar Raissi, Alireza Yazdani, and George Em Karniadakis. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, (67(6481):10266) 10266611000, 2020a. ISSN 0036-8075. doi: 10. 1116/science..aw.47111
Maziar Raissi, Alireza Yazdani, and George Em Karniadakis. Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations. Science, (67 2020b.
12

---

Published as a conference paper at ICLR 2022
Markus Reichstein, Gustau Camps-Vall Bjor Stevens, Martin Jung, Joachim Denzler, Nuno Car- valhais, et al. Deep learning and process understanding for data-driven earth system science. 56677733366 204, 2019.
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.
M-H Rio and R Santoleri. Improved global surface currents from the merging of altimetry and sea surface temperature data. Remote sensing of Environment, 216:770-785, 2018.
Pierre C Sabatier. Past and future of inverse problems. Journal of Mathematical Physics, 41(6): 4082-411 4082-4124, 2000.
Steindor Saemundsson, Alexander Terenin, Katja Hofmann, and Mare Deisenroth. Variational inte- grator networks for physically structured embeddings. In International Conference on Artificial Intelligence and Statistics, pp. 3078-3087. PMLR, 2020.
Priyabrata Saha, Saurabh Dash, and Saibal Mukhopadhyay. Phicnet: Physies-incorporated convolutional recurrent neural networks for modeling dynamical systems. arXiv preprint arXiv:2004.06243 2020.
Omer San and Romit Maulik. Machine learning closures for model order reduction of thermal fluids. Applied Mathematical Modelling, 60, 04 2018. doi: 10. (0001///..pm.2000 8.03.037.
Russell Stewart and Stefano Ermon. Label-free supervision of neural networks with physics and domain knowledge. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 31, 2017.
Daniel J Tait and Theodoros Damoulas. Variational autoencoding of pde inverse problems. arXiv preprint arXiv:2006. 15641, 2020.
Michael L Thompson and Mark A Kramer. Modeling chemical processes using prior knowledge and neural networks, AICHE Journal, 40(8):1328-1344 1994.
Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Accelerating Eulerian fluid simulation with convolutional networks. In Doina Precup and Yee Whye Tch (eds.), Pro- ceedings of the 34th International Conference on Machine Learning. volume 70 of Proceedings of Machine Learning Research, pp. 3424 3433. PMLR, 06-11 Aug 2017.
Kiwon Um. Robert Brand, Fei Yun, Philipp Holl, and Nils Thuerey. Solver-in-the-loog Learning from differentiable physics to interact with itcrative pde-solvers. In Neural Information Process- ing Systems (NeurIPS), volume 33, pp. 611 1-622,, 2020.
John Neumann. Functional Operators (AM-22), Volume 2: The Geometry of Orthogonal Spaces. Princeton University Press, 1950. ISBN 18-11-0008-818- doi: https://doi.or//10.
1515/9781400888199
Nils Wandel, Michael Weinmann, and Reinhard Klein. Learning incompressible fluid dynamics from scratch towards fast, differentiable fluid models that generalize. In International Confer- ence on Learning Representations, 2021. URL
Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach, and Vipin Kumar. Integrating physics-based modeling with machine learning: A survey. arXiv preprint arXiv:2003.04919, 2020.
Yuan Yin, Vincent Le Guen, Jérémie Dona, Emmanuel de Bezenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting Physical Models with Deep Networks for Complex Dynamics Forecasting. In ICLR, 2021.
Chih-Chich Young, Wen-Cheng Liu, and Ming-Chang Wu. A physically based and machine learning hybrid approach for accurate rainfall-runoff modeling during extreme typhoon events. Applied Soft Computing, 53:205-216 2017.
13

---

Published as a conference paper at ICLR 2022
A DISTANCE
A. DISTANCE BETWEEN DYNAMICS
We here give the definition of the distance d. Let a and be two functions of (RP,, We consider the distance:
d(u, u) Ex || exX)l22
(18)
Naturally, eq. (18) verifies the triangle inequality, the symmetry and the positiveness, Moreover, in this case, for all functions f. d(.: f) is convex. Indeed, for u, two functions, and [0, 1]:
d(Au + (1 x).. = Ex~px [u( + (1 X)((X) - f(x0ll2
-xxxpx Auu() (X) Af(x) - (1-A)((X) + (1 (1-))o(X)
I AE x=xxx [lxx - ((x)ll2 - X)=
((X)ll
Hence the convexity of f), This consideration suffices to ensure the convexity of S₁ and S. defined in section 3.
A.2 DISTANCE BETWEEN FLOWS
Consider the ODE with X(t), Xo IR::
dx((t) = f(X(t), X(t = 0) = Xy
(19)
dt
Equation (19) admits unique solution soon as f is Lipschitz. We note X this solution. Then, we can defined the flow oy such ODE as
[0.7] xR" RP
(20)
t, Xo +
= (t)
With the definition of eq. (20), we can define the distance between two flows of ODE as:
ds (o( E Exomxx [ [|ou(t, of Xt)ll dt
(21)
de is positive and symmetric. Let du. du be two flows, we have the triangle inequality:
dolon, Q5) ExOmpx I' (ta,(t, st Xo)l dt
xxxxxx 110a11, do(1, Xo) (1, Xo) + (1, di
EXo~xx 110(11, - (t,,, Xo)ll + 110 (1, Xo) + (,, Xo)
S doldr: dv) dol 25)
Let Øy be fixed, we also have the convexity of do(,, 6r) with respect to the first argument. Indeed for X []:
I [t (t) +
do(2 + (1 - 11-A)) f) = Exoxxx,
(1 +11
A)o. - (t, dt
= =Fxorpxo ( (1, ) + +11 -))0 - At Xo) (1 - -))og(t, Xo)l de
S \ Sdd + (1- (1,))ds(Ornd))
However, in this case the convexity is not ensured with respect to and This is the reason why for theoretical investigations, we consider the distance d instead of dg.
Nonetheless, de do 22) 0 Du = Φ5
u=
14

---

Published as a conference paper at ICLR 2022
B REMARK ON ADDITIVE DECOMPOSITION
First, note that in the case of metric space the decomposition as defined in eq. (2) always exists.
We now detail an intuition for the well-posedness of such decomposition.
Let H₂ be a closed convex subset of functions of an Hilbert space (E, <,,,,>>, and I the function we want to approximate with partial knowledge (represented by the space of hypothesis Hz). Then, thanks to Hilbert projection lemma, we have the uniqueness of the minimizer of minge Il all, ie exists one unique hgg Hr such that: Vig E H4, |l h& II gll-
Finally, the additive decomposition hypothesis presents remarkable advantage in the case of vector space. Indeed, Hy (closed) vector space, let Hg be its supplementary E, then we have the uniqueness in the decomposition: f frLe + fH+ where fuz H4 and fH. E Hx-
The existence and uniqueness flowing directly from the additive decomposition hypothesis, this can explain why such assumption is common when bridging ML and MB hypothesis.
C
UPPER BOUNDS
C.1 DERIVATION OF EQUATION (4)
The first upper bound is simple use of the triangle inequality:
d(h, f) d(h,, f) + d(hu, f) - d(hu, f) S d(h;, f) + [d(h, f) d(hx,
S d(hh, f) + d(h,
C.2 DERIVATION OF EQUATION (5)
To derive the second upper bound, we assume that J comes from an overall dynamics f obeying the additive decomposition hypothesis of eq. (2) so that f and f verifies: f S8" + frr. First, with computations similar to eq. (4), we have:
d(h, S SP", + A(
(22)
Then:
d(h, fr)) = d(h, f d(hk; fb")
d(hx. ft
< d(hx, 1 + |d(h, fo
d fk()
d(hk: f + |d(h, fo
d(h, fh"
d(hk: f(r))
d(hx; 18" per) d(h, f(") ) + |a(hs, fk") - d(h, hts)
< d(hk: fh") + d(ftr ffh + d(h, ha)
(23)
Combining Equations (22) and (23), we retrieve eq. (5): d(h, f) S d(hs, ff") + d(h, h&) + d(fpr, .ft") + a(f",
(24)
and we have: I d(f") ff"") + d(for, . .f). I is constant of the problem that cannot be optimized.
C.3 UPPER-BOUND USING AUXILIARY DYNAMICS for
Let f be the dynamics of model data, dynamics f) and the error made by h on model data (with dynamics fPr) via:
link up the error made by h true data (following
d(h, f) < d(h, fP") + d(frr, .f)
(25)
Thus pre-training on auxiliary data of dynamics for amounts to control the term d(h, frr in the upper-bound of eq. (25).
15

---

Published as a conference paper at ICLR 2022
C.4 SELF- SULE-SUPERVIIION
Let - hu + hu be the function to learn and Ge the recognition network providing an estimate 0% of the parameters from an initial sequence (X: ................................ . -iu+kke) This learning setting corresponds to how velocity fields learned from consecutive measurements of the tracer fields 7 in section 4.2.
To compute d( fk") in the case where for = h" where h* = h th, is a learned model, we rely on the computed 6% associated to hg (thanks for example the algorithm of section 3.2 associated to eq. (4)) to generate a synthetic dataset with achievable supervision in the space the parameters Og.
From a real initial sequence (Xi, X tott) we can estimate the unknown parame- ter 0% associated to sequence i with the recognition network G learned with h*, i. 0g X", +kad)) Then, integrating from the initial condition Xi, we generate tra- Sampling the space of initial conditions, we obtain synthetic dataset: ((X) (1)))))).. (.),.... (Xm, (Tr) enabling us to perform self-supervision for Gys. Let 0% be the parameters estimated by G from the simulated (X) X x(ts))) make the following approximation:
jectory of known parameters of with dynamics h denoted by: X = (X "t,....,
d(hu, fk)
Σ
0%
0%
(26)
D PROOFS
D.1 NOTE ON THE CONVEXITY OF S₁ AND S..
Convexity of Sk
Proof. Let u, E Sy: d(tu + (1 f), = |tu + (1 t)e -|1= = 1|u - 1f + (1 - -)0 - (1 t)lll
< mmn + (1 1)1 m
Hence the convexity of Sg.
Convexity of S..
Proof. Let E [0, 1] and u, v E Sa.
d(hg, by + tu + (1 t)u) d(0, tu (1 fle)
< (d(u, 0) + (1 (l)(w,)))
S u12
Hence the convexity of S..
D.2 ODE IDENTIFICATION
Consider the following set: SA (X(t) e C³ ((O, RP) such that: 3A MR, AX), where T > 0.
SA not convex set. Consider and in SA, and consider A and A" that u'(t) A'u(t)
and et( = (t)( For X E|0, 1[: we have:
[Au + (1 - A)o}' Au' (1 1 -))o'
= AA" + (1 A)44'
In general the last term is equal A1u+(1-A)) (1 A)u), for some matrix A>u+(1-Jje
Thus SA is not a convex set. However, discretizing the trajectories and employing simple integra- tion scheme leads considering the following cost function:
((A) = Σ < 14)
+
(27)
16

---

Published as a conference paper at ICLR 2022
As least square regression problem, [(A) is convex with respect to A. A least square regres- sion setting can also be recovered using more complex integration schemes, or several time steps integration.
D.3 PROOF FOR POLL--SSEDNESS OF EQUATION (7)
We set ourselves in the Hilbert space of squared integrable functions with the canonical scalar prod- uet (R", R²), >,>} For further consideration on such functional space we refer to (Droniou, 2001).
We assume that Hy. hence Sş is convex and a relatively compact family of functions.
Convexity of S5 + Su Let S = S + S, = (f B3f E Su, E Sw, f = fk + fa).
Let f.g E S and e]0,
A)gu E Sk: + Su A (1 A)g Aʃk + N)gx Afa + +11 [1
Hence the convexity of S.
CLOSENESS OF S, show that S, is closed set. Indeed, S, ([0, pul), where g(u) ||u/|, Because is I-Lipschitz (using the triangle inequality). is continuous. Therefore Su is closed set as the inverse image closed set by continuous function.
Sequential Limit We now show that S is a closed set thanks to the sequential characterisation: let f" converging sequence of elements of S and denote its limit. We prove that f converges in S.
Because Vn, fn E S, we have: f = ff + fü where f E Sy and ft Sg.
Thanks to the relative compactness of Sp, from J so that l + fa Sz.
extract converging sub-sequence, of indexes nj,
Because f + ʃ, the sub-sequence f converges: g +
By definition, f" is sequence of S. and we also have that: f right member of the equation converges (as a sum converging functions), the left member of the equation converges i.e. for converges.
- f
Because the
Since S, a closed set für converges S... We write fu its limit. Therefore, f fd) fk + fk fa S,. Hence, f = fu + fo with fu E Su and fo E Sy.
Therefore S is closed set.
Finally, we can apply Hilbert projection lemma on the closed convex set S and retrieve the unique- ness of the minimizer eq. (7).
The relative compactness of family functions a common assumption in functional analysis. For example, in the study of differential equation Cauchy-Pean theorem provides the existence to the solution of an ODE under the assumption of relative compactness.
Remark
Also, Ascoli theorem provides the relative compactness of family of function F under the hypothesis of the equi-continuity of F and the relative compactness of the image space (xx) (f(x)) F).
D.4 PROOF OF PROPOSITION 2
We now set ourselves in the Hilbert space (C3 []] R²), <,>) of squared integrable functions, where <, > the canonical scalar product of [ Ti, R*).
Proof. Let A be given invertible matrix. We consider the following space Sp {X C. ([0, T]. RP) such that: 3D e R = AX + D and X(t = 0) Xo). where T > 0. We show that Sp is closed convex set.
17

---

Published as a conference paper at ICLR 2022
Convexity Indeed, let À eJo, 1 and u, Sp. Au + (1 - A)u is differentiable and:
[Au + (1 A)v)J' Au' + A)e' A(Au + (1 - 2)e) + D,
Where D = AD + a. Hence Au + (1 A)U E Sp.
Closeness via Affine-Space To prove the closeness of Sp. we prove that is an affine space of finite dimension.
Let g the application that to any vector Rd associate the solution X
Let Do E RR we show that ID. D - g(D + D) (D)) is linear application.
Naturally, for IDo (0gp) Oga. Then for D + Ogr have:
9Du (D) eAl(( (Xn + (D + D)) (Do e (D=) +A-1Do
EA -AAAA-1 A-1D
Therefore ID. is linear function and g is an affine function.
Moreover, g is injection. Indeed, if two functions equals, then they have most one inverse image by thanks to Cauchy-Lipschiz theorem. Therefore it defines a bijection of R4 in g(R)). Since, Sp g (R Sp is affine space of dimension p and g is continuous in particular for the canonical norm induced on ([0, T], RP). Therefore Sp is an affine space of finite dimension and a closed set.
Finding a Unique Minimizer We conclude by applying Hilbert projection lemma: our problem
of minimizing S"' [T) - amounts an orthogonal projection problem. Because Sp is a closed convex set, we have existence and uniqueness of such projection. Therefore, it exists unique function XD SD and unique vector D minimizing its distance the function X*. O
D.5 ALGORITHM LINEAR SETTING
We detail in Algorithm 2 the alternate projection algorithm in linear setting. We denote = (Xi bo+ A Xi (nrr) and X i Xitt0r)) For readability purposes we set At =
Algorithm 2 Alternate estimation: Linear Setting
Result: (p.,(()) E Rp
0, D = 0, A.
A = minalyy - XAII
while | D > and | 1 > do
D minp []] X D1 + >|DD1|
A = minA || - XA - D >|lY -
A k++
end
D.6 PROOF TO PROPOSITION 3
Naturally, one could estimate jointly and using least square regression. However, the idea is to verify the convergence of such alternate algorithm in simple case. We conduct the proof for the first dimension of y to lighten notations, meaning that we are regressing the first dimension Y against the X.
similar reasoning for the other dimension completes the proof.
Proof. We first give the analytical solution for D. Let A" be fixed.
Estimation of D Consider:
D12 CD I - +
(28)
where D (d....,d) E IRQ, For Q samples, find d so that
01
04 0:
18

---

Published as a conference paper at ICLR 2022
Q
ØL.
= =2*
U/
X, d) + 2\d
0d
Q
4 Qd + Ad =
Σ
(y) - X,A")
1=1
+ + Σ
(y) X,A*)
Y-XA
e 1+
where Y Y 1 Ei X,
Estimation of A Let D be fixed and consider:
=A= = Y XA DD12 + >/1x - XAI2
(29)
Similarly, we aim to cancel the first derivative of LA with respect to all parameters of A = (ai,.., ap):
θLA 0uj
0 + - 2
Tiij (1l 4000.0 - + Aptip d)
i=1
-2
Fig (U: a03400 aptipp)
-X XA D) 2 XA)
(1+7)) X'(Y D) -XX' =0
*(1+7)x'XA= X'(X (
xx - 1+y ((++)) - D)
(30)
where B X.X. Equation (30) indicates that as soon converges, A" converges. Thus, we now prove the convergence of (D*). Then, for n > 1 consider:
1
D
D*11
Y
X Y X A
T+ T/Q
1
X(A" An-1]
1 1 X/Q 1
XB-IX
Dd) D=-1)]
(1+ a/Q)(1
XB-IX Dn-1 Dt- D
(1+ \Q)((1
K
(1)\/)(11 1 17
where = 1I=---X-'XII
Therefore, for X, y. sufficiently large, (1+2/Q11+7) 1. | Dn-1ll converges as positive decreasing sequence. Finally, the sequence of (D.) converge and so the sequence of (An).
In conclusion, the proposed algorithm converges.
19

---

Published as a conference paper at ICLR 2022
E DATASETS
In this section, we provide exhaustive simulation details for the damped pendulum, Lotka- Volterra, and both geophysical datasets.
E. DAMPED- DAMPED-PENDULAN PENDULUM
For the damped pendulum data, eq. (15) is integrated with At 0.2s using a Runge-Kutta 4- 5 scheme from t 0 up to = 10s. Both the pulsation and the damping coefficient k are fixed across dataset. We generate 100/50/50 sequences respectively for train, validation and test sampling over the initial conditions so that (0, 0) (11 x [-0.1,0.1]..
Small Oscillations To linearize the pendulum, we consider the small oscillations regime and take the initial conditions so that (0,0) (1 / n /6) x 0.1.1 In that case eq. (15) writes as:
d di
i/4i
(31)
and following notations of section 3.3, we have: DA = 0 and A = (-> 0)
E.2 LOTKA- VOLTERAA
For Lotka-Volterra data, eq. (16) is integrated with At = 0.05 using a Runge-Kutta 4-5 scheme from up 20. All parameters a, 7, 8 set across the dataset. We generate 100/50/50 sequences respectively for train, validation and test sampling over the initial prey and predators populations so that (x, v) ~ 4 ((0,2 2)²).
Practical Issues and Adaptation Assuming that a and y have positive values makes the following problem arises: the trajectories defined by h4 for the prey are unbounded, whereas the trajectories defined by eq. (16) are. Minimizing d(h4; f) over long term horizon will lead in an underestimation of to match the bounded behaviour of true data. Therefore, we enforce d(hg. J') on the prey component as soon as the number of predator is small. In practice, we set this threshold to 0.15.
E.3 GEOPHYSICAL DATASETS
We present in this section introductory tools for the understanding of the fluid dynamics data pre- sented in section 4.2. We first introduce the physical modeling of ocean dynamics. Then, we outline the Adv+S dataset simulation which draws from ocean modeling. Finally, we introduce the Natl dataset and the proxy data used in the experiments.
Introduction To Ocean Modeling The increase in ocean observations thanks to satellites and floats enabled great development in Earth modeling over the last decades. The ocean circula- tion, that is the current velocity fields dynamies, now realistically modeled in tri-dimensional structured models such as NEMO (Madec, 2008).
Such models rely on in-depth physical knowledge of the studied system and its representation through partial differential equations. Integrated over depth, the equations associated to the transport of the Sea Surface Temperature T by varying horizontal velocity field U can be written as:
0T
01 V.(UU))T)))))) DT + FT
(32)
aU Ōt = (U.V)U + fAU g'vh+ DU + FU
(33)
where f is the Coriolis parameter, h the depth of the surface layer obtained from surface height (SSH) observations, g' the reduced gravity which takes the stratification in density the ocean into account such that g' R 1013 In two-dimensional setting. V(TU) refers to the advection of
20

---

Published as a conference paper at ICLR 2022
scalar quantity T by velocity field U = (u, v) and writes V(TU) terms, referred to as D and the forcings FT/6 are not known.
01 Ox u + 0T The mixing
In the context of the presented work, the physical state Ze - (T2, Ut), fx and fy from eq. (1) can be interpreted as follows: fx represents the dynamics of the observed r.c.c. fx = (TU) D" F in eq. (32). fy represents the dynamics of U in eq. (33), i.c. (U, -U. + fAU g'Vh+ DU + F
Whereas observed by satellites, U is not known. However, the Sea Surface Height (SSH) could be used to compute coarse estimates of U. Indeed, under hypothesis such as stationarity (BL 0), incompressibilityy (U.V)U = =)), forcings can be omited. this case, eq. (33) can be rewritten
into
JAU=-9'
(34)
When projected onto x and y axis, eq. (34) becomes
Oh -fo - dx
Oh Oy
fu=
(35)
Note that eq. (34) and eq. (35) not hold at fine scales as the stationarity and incompressibility assumptions only hold at large scale. In this case, the SSH h can be regarded as stream function.
Both datasets considered in the paper follow the same equations approximating the tracer equation (eq. (17)) inspired by eq. (32):
0T =
(36)
We study the equations 32 and 33 in an incremental approach. In the following parts, we describe how T, U and S are computed in both datasets Adv+S and Natl.
E.3.1 ADV+S
We first investigate a dataset generated following simplifying assumptions (Adv+S). We don't 't rely on true U and S, we instead build them so that they correspond to our hypothesis.
Building a Velocity Field U Under stationarity and incompressibility hypothesis, U can be ap- proximated from a stream function H. Note that, in this dataset, H. is not equal to the SSH h, is simulated following (Boffetta et al., 2001):
y B(I cos kx
H(s, y, 1) tanh|
V BC) x sin + cy:
(37)
As introduced precedently (see eq. (34)), eq. (33) can be simplified and we compute U = (u, u) so that it follows:
HH
TH
(38)
Oy
Or
Note that B varies periodically with time according to B Bo + cosswwt Φ). We compute 10 different velocity fields sampling random parameters Bo, ki C, m,,,,,
Building Source Term S In eq. (32), the diffusion term DJ is omitted. We generate the source term S so that it represents the forcing term FT in eq. (36). To illustrate heat exchanges, we draw from Frankignoul (1985). This source term is a non linear transformation of U (u, v) multiplied by the difference between the ocean temperature and a reference temperature:
S(U,T) w. x (T Te) where We = {
0 if 0H
< 10-1 otherwise.
at
where T. is the sequence mean image (computed without source).
Dataset Generation Using computed U and S, we integrate eq. (36) with At 8640s over 30 days, using a Semi-Lagrangian scheme (see explanations below). We generate 100/100/200 sequences respectively for train, validation and test sampling over the initial conditions, which are images of size 64 x 64 sampled from Natl dataset. Finally, for integration, we impose East- East-West periodic conditions, implying that what comes out the left part re-enters at the right, and reciprocally. We also impose velocity to be null on both top and bottom parts of the image.
21

---

Published as a conference paper at ICLR 2022
Semi-Lagrangian Integration Unlike Eulerian scheme, relying time discretization of the derivative, the semi Lagrangian scheme relies on the constancy of the solution a PDE along characteristic curve. Consider a solution to the advection equation, i.c. eq. (36) with S = 0. The method of characteristics consists in exhibiting curves (z((), (())) along which the derivative of the solution T is simple, i.e ar (a(s), t(s)) = For ID constant advection scheme, computations lead to:
as
dt
s=t
ds dr
x x= + U
ds
giving therefore, t(z,t) = To(x Ut), linking the value of the solution at all time to its initial condition. Therefore from a single observation at to, it suffices to estimate the original departure points 20 Ut to infer the prediction at
However, when is not constant time, the method remains doable, not along characteristic lines defined by: (ao + Ut), but along characteristic curves which are given by:
dt ds
dr ds
(x,
(39)
great deal the semi-Lagrangian literature involves solving correctly eq, (39). We use the con- ventional mid-point integration rule and the semi- aani--agrangin is implemented using Pytorch function gr following in (Jaderberg et al., 2015). Further developments can be found for exam- ple in (Diamantakis, 2014).
E.3.2 NATI.
This second dataset depicts the actual ocean circulation, i.c we consider both eq. (32) and eq. (33). In this case, assumptions made on U and S represents both diffusion terms D and forcing terms FT, We access daily data over year of ocean surface temperature of the North Atlantic observations model resulting from (Ajayi al., 2019) ! The dataset covers 2300km 2560km zone at 1.5km resolution, in the North Atlantic Ocean.
In this real-life dataset, sea surface height (SSH) partial derivative with respect to æ and y serves as proxies to the (unobserved) velocity fields U. Indeed, recall that simplifying hypotheses led us to eq. (35).
We divide the Natl zone into 270 patches of size 64 x 64. For each region, we extract sea surface temperatures, velocity fields, source terms and height variables. We sample 200/20/50 sequences of 1 year, for respectively train, validation and test. In this case, At = 86400s (1 day).
TRAINING DETAILS
All experiments were conducted on NVIDIA TITAN GPU using Pytorch (Paszke et al., 2019).
Interpretation From eq. (4), two independent terms appear justifying an al- ternate projections approach.
First, we highlight that strictly minimizing d(hu, f) biases our estimation of hg. However, it may yield good estimation of hgg provided that f&x contributes significantly to the prediction of f. Hence, we interpret this loss as an initialization loss. Thus, in most applications, we progressively decrease its magnitude along training as detailed in appendices F.1 F.3.
On the other hand, d(h,,, 0) aims at constraining the free form function hy to make its action as small possible. We interpret this loss stability penalty.
Details available at http://mmeemggggonggiiiiooooooot
22

---

Published as a conference paper at ICLR 2022
Finally, aiming to recover exact trajectories of observations, proceed suggested in (Yin al., 2021) progressively increasing the ayper-parameten associated to d(h,
The practical implementation is summarized in the following algorithm:
Algorithm Alternate estimation: Practical Setting
Initialization: 0% 0, 0% min d(hg, f'), Ana Anaa Ana
for epoch = 1 Nepoche do
for batch = 1: B3 do
0ht =0 - - -o, + Ana {(hx))
end for batch- Bv: : Bu do
gutt1 t" - -I Vo. Dn,h,, + An d(hu, 0))
end
A = 72i Anee
1 Ani An An
72
end
F.1 DAMPED PENDULUM
Architecture Details The physical parameters to be learned is scalar of dimension 1, and h, is a 1-hidden layer MLP with 200-hidden neurons with leaky-relu activation.
Optimization For this dataset we use RMSProp optimizer with learning rate 0.0004 for 100 epochs with batch size 128. We supervise the trajectories up to t At x 50, i.e we enforce dø over (to At, .., to + 50At). Overall the number of optimization subsequences for training is 17000. alternate projection on S₁ and S,, by descending the gradient 10-batches on hę then 10-batches on h..
Hyperparameters We initialize Ana 0.1 and decrease it geometrically down to Au. 0.001. We initialize A 0.1 and increase it geometrically to λ = 100. X,., is fixed through training 0.1.
The hyper-parameters were chosen by randomly exploring the hyper- parameters space by sampling them so that A ~ U11,0.1.... 10-4). We select the ones with the lowest prediction errors, i.e with lowest (h,
For the ablation study of Table 1, we set to 0 the prree-parameters associated to the non-considered loss.
The training time for this dataset hour.
F.2 KOTKA-VOLTEREE
Architecture Details The physical parameters to be learned is a vector of dimension 2 accounting for (a, 3) eq. (16), and ht is a 1-hidden layer MLP with 200-hidden neurons with leaky-relu activation.
Optimization We use Adam optimizer with learning rate 0.0005 for 200 epochs with batch size 128. Overall the number of sequences for training is 15000. We supervise the trajectories up to t At 25, i.e we enforce dy over (to At, .., to + 25tt). We alternate projection on Sş and S. by descending the gradient 10-batches on h2 then 10-batches on hy.
Hy Hyperparameters We initialize An 0.1 and decrease it geometrically down to Ane 0.001. We initialize X = 0.001 and increase it geometrically up to X 1. An, is fixed through training 0.001.
23

---

Published as a conference paper at ICLR 2022
The hyper-parameters chosen by randomly exploring the hyper-panne ryper-parameters space by sampling them so that À lowest d(,h f)).
10 "). We select the ones with the lowest prediction errors (i.e
For the ablation study of Table I, we set to 0 the hyper-parameteon associated to the non-considered loss.
The training time for this dataset 2 hours.
F.3 ADV+S
Architectures Details The physical parameters to be estimated are the velocity fields U, of dimen- sion (2, 64, 64). As U varies over time, follow data assimilation principles map a sequence of 4 consecutive measurements of the tracer field T to the associated velocity field (Gaultier et al., 2013). To do so, we parameterize a recognition network G by U-net with most 512 latent chan- nels also following the implementation of (Isola et al., 2017), taking as input a sequence of 4 time steps of T: (Tip.-.. +30+)))) The residual dynamics h, is learned by convolutional ResNet, with residual block taking entry the same sequence of T. We implement hgg via a semi-lagrangian scheme, taking as input T. and the estimated U, to predict T++1-
We use Adam optimizer with learning rate 0.0001 for epochs with batch size 32. We supervise the trajectories up to f = At x 6, i.e we enforce de on + tto+) Overall the number of sequences for training 36800. We alternate projection S₃ and Sy by descending the gradient 4-batches on h4 then 6-batches on ha..
Optimization
hu
Sio
Th,11
h4
Ti -sto
U
U.)
To
Figure 3: Best viewed color. Schematic view of our model the context section 5.2, applied on Adv+S dataset.
Hyperparameters, setting of eq. (4) We initialize An, = 0.1 and decrease it geometrically down to Any = 0.00001. We initialize A 0.01 and increase it geometrically every epoch up to Ap. = 1000. Ao. is fixed through training at 0.1. We select the hyperparameters with the lowest prediction emros (i.e lowest f)). For the ablation study of Table 1, we set to 0 the hyper-parameters associated to the non-considered loss.
The training time for this dataset 8 hours.
F.4 NATL.
Architecture Details The architectures this setting are identical to the ones described in ap- pendix E.3.
Optimization We use Adam optimizer with learning rate 0.00001 for 50 epochs with batch size
32. Overall the number of sequences for training is 67000. We enforce dg over 6 time-steps, i. we supervise the predictions timesteps: (1o + Al, to + 6A1). We dropout in both Gp and h...
24

---

Published as a conference paper at ICLR 2022
Hyperparameters, setting of eq. (4) For this setting, X geometrically increases from 0.01 up to 100. We initialize X = 0.1 and decrease it geometrically down to A = 0.00001. Xp. is fixed through training at 0.1. We alternate projection on S₁ and S,, by descending the gradient 10-batches on both h& and ha.
The selected model is the with lowest prediction errors validation set (i.e lowest d(h, f)), sampling uniformly the hyperparameters: X ~~
Hyperparameters, setting of eq. (5) Because the dynamics of Natl is highly non linear and chaotic, we follow (Jia et al., 2019) and first warm-up the parameters recognition network Ge on the velocity fields proxies for epochs. For this setting, Å geometrically increase from 0.01 up to 1. Abe is set equal to Xς. An.. is fixed through training at 0.01.
After warm-up, we alternate projection on S₃ and S, by descending the gradient 100-batches on by and 300 on hu. In this setting of eq. (5), the selected model is the one with lowest d(h, f) + d(hy: fk") error, sampling uniformly the hyperparameters: ~ ((1.0.1,...... .0.
The training time for this dataset is 12 hours.
Baselines We train NODE (Chen et al., 2018) and Aphynity (Yin et al., 2021) on both the Adv+S and Natl dataset. For the training of Aphinity, we the learning at 0.0001 and train 30 epochs. We initialize A = 0.01 and increase it geometrically every epoch up to A = 100. A fixed through training at 0.1. For the training of NODE, set the learning rate at 0.00004 and train on 50 epochs. To perform prediction, we first encode the 4-consecutive measurements T (as 64 x 64 state) then learn integrate this state in time thanks to network h. h is 3-layer convolutional networks, with 64 hidden channels. It is integrated using RK4 scheme available from https: / com/ wttpgriceeeeettccccccccccc ffeq.
G ADDITIONAL RESULTS AND SAMPLES
G.1 RESULTS FOR PENDULUM AND LOTKA- VOLTERRA DATASETS
We provide respectively in figs. 4 and 5 phase diagrams for the damped pendulum and Lotka-Volteraa experiments. Both graphs in the phase space indicate that the trajectories and their nature are well handled by the learned decomposition, providing periodic phase space for and a converging spiral for the damped pendulum (fig. 4).
(fig. 5),
G.2 RESULTS FOR ADV+S AND NATL
In this section, we provide additionial results on both Adv+S and Natl datasets. A thorough abla- tion study (table 4) gives results with constant hyperparameters A and db, (row Vanilla Optim), which validates our hypee-parameters interpretation. Indeed, the results are better when respectively increasing and decreasing Å and λb.- Besides, the row Ours eq. (5) refers to training performed as introduced in appendix C.4 with for ht trained on eq. (4). Figure 7 shows predictions up to days on the Adv+S data. Finally, figs. 9 and 11 provide results on Natl dataset associated to training relying on both eq. (4) and eq. (5) and with NODE (Chen et al., 2018).
25

---

Published as a conference paper at ICLR 2022
true prediction
0.05
0.04
0.03
0.02
0.01
0.00
-0.01
0000
-0.03
-0.4
-0.2
Figure 4: Damped Pendulum Phase Diagram. The true phase diagram (bluc) and learned (orange dashed) are close, indicating consistency in the prediction
true prediction
2.5
y (pred pop)
1.0
0.5
(prey pop)
0.5
2.0
Figure Lotka-Volterre Phase Diagram. The true phase diagram (blue) and learned (orange dashed) are close, indicating consistency in the prediction
26

---

Published as a conference paper at ICLR 2022
Table Ablation Study Adv+S, We report the MSE 100) the predicted observations T, the estimated velocity fields U and the residual source term S over 20 time steps from an initial datum to. Unlike alternate training, i.c. Algorithm I, "Joint" rows refer to the simultancous optimization of hg and ha-
tu
6
tu + 20
Training
Models
T
U
S
T
U
S
Ours (U known)
0.52
n/a
0.19
2.0
n/a
0.32
Ours eq. (4)
0.74
1.99
0.17
8,49
2.26
0.31
only d(h, f)
1.02
4.08
0.19
10.59
4.19
0.32
d(h, f) + d(h,, f)
1.02
3.66
0.19
11.42
3.84
0.34
Alternate
d(h, f) + d(h, ha)
0.77
2.38
0.19
9.5
2.45
0.34
Ours eq. (5)
0.75
2.77
0.17
8.36
2.84
0.29
Vanilla optim.
1.51
3.77
0.3
13.33
4.1
5.15
Ours eq. (4) only d(h, f)
1.44
3.3
0.3
12.82
3.5
0.5
Joint
1.38
6.96
0.39
11.9
7.09
0.54
1+1
t+2
1+3
1+4
1+5
I+5
1+1
1+2
1+3
1+4
1+5
truth
Ttruth
buth ostimation
buth estimation
estimated
estimated
fruh estimation
ful estimation
Figure 6: Best viewed color. Estimations, targets and differences between estimations and targets on T, = (u, and S for Adv4S. Each column refers to time step. ranging from I to half-days. On the left, true and estimated U (u, v) over 6 time steps, and differences between largets and estimations. On the right, prediction of T and S over 6 time steps, and differences between targets and estimations.
27

---

Published as a conference paper at ICLR 2022
1-4
13 12
M
122
114
1+B
IF2
1+6
1+B
Inputs
Outputs
Targets
Figure 7: Best viewed in color. Estimations and targets T. U (u, v) and S for Adv+S. Each column refers to time step. ranging from half-days. On the left, sequence of T inputs (4 time steps). middle, prediction T, (u, v) and S over time steps. On the right, true T, and S over time steps.
1+1
1+2
1+3
1+4
1+5
I+2
1+3
1+4
1+5
1+1 1+2
t+3
144
1+5
8
Targets
Outputs
Target-Output
Figure 8: Best viewed color. Estimations, targets and differences between estimations and targets on T. = (u,v) and S for Adv+S. Each column refers to time step, ranging from I to half-days. On the left, true T. U and S over time steps.. In the middle, prediction of T,, U (u, v) and S over 8 time steps. On the right, differences between targets and estimations.
t+1
1+2
t+3
estimated (oq (4))
v estimated (eq (4)
estimated (eq (5))
estimated (eq (5))
Figure 9: Best viewed in color. Sequence of estimations on U = (u. v) for the Natl data. The second and third respectively refer to training according to (4) and eq. (5). The loss term d(ha, f("") in eq, (5) enables our model to learn more accurate velocity fields than when only trained following eq. (4).
28

---

Published as a conference paper at ICLR 2022
t+1
t+2
t+3
t+4
t+5
t+6
True SST
Ours eq.5
Neural-ODE
Figure 10: Best viewed in color. Sequence of prediction T for the Natl data. Contrary to our model (row eq. (5)), NODE (row Neural-ODE) struggles predict any motion in T.
I+2
T
Targets
Outputs
Figure 11: Best viewed color. Sequence prediction on T, u, v, S the Nail across days trained
using proxy data according to eq. (5)
29


Reviews:
[{"review_id": "dO2QIo_CTq3", "reviewer": "Reviewer_YsB3", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "The paper is generally clear and well-written and the proposed approach is interesting and well supported by experiments.\n\nWith that said, I have a few comments:\nTheory:\n* Proposition 1: It is not clear how the uniqueness of the solution is connected to the accurate identification of h_k. My main concern here is that the paper does not place any condition on the relation between the function spaces H_k and H_u. Contrast this to the linear case where A is invertible and thus cannot realize a constant derivative without D_A.\n\n* It appears that Section 3.3 abandons the \"realizable setting\" assumption, where the true function is assumed to be in the hypothesis space (which is good). I suggest that the authors put more emphasis on this. Otherwise, the transition can be confusing.\n\n* Proposition 2: I am not sure of the implication of the proposition. It seems to only suggest that for each \\hat{A} there exists a unique D_A but does not suggest the uniqueness of A, which is what we would be interested in for interpretability.\n\nExperiments:\nA question that arises is whether alternating optimization is needed in practice or if joint gradient descent will be as good. It would be interesting to see this experiment.\n\nClarity:\n* Section 4.2 needs much larification:\nEquation 17: The del operator needs to be defined. What does TU mean when U is a tuple (u, v). Why do you use parenthesis for S(U) but not TU. Also, that does T_{t-k} mean? Doesn't k stand for \"known\"?\nThe fact that \\theta^* is a function in time (rather than a constant vector) can be confusing and should be spelled-out.\n\n* After eq (16) \"\\theta^* = (\\alpha, \\beta)\", I believe \\beta should be \\gamma?\n\n\n\n", "summary_of_the_paper": "The paper proposes a method for hybrid model-based/machine-leanring learning, wherein a model is decomposed into an interpretable parametric prior and a residual (typically represented by a neural net). The paper demonstrates that prediction error minimization does not accurately identify the parametric component and proposes an alternating optimization method that augments prediction error loss with component-specific losses. The paper also introduces a variant where auxiliary data can be used to define the loss.\n\nTheoretical guarantees are provided where auxiliary losses are given as upper-bound constraints and exact minimization is performed. Also, uniqueness and convergence guarantees are provided for the linear approximation case using a more practical algorithm.\n\nExperiments show the efficacy of the method in linear and non-linear settings, in both identification and prediction.", "summary_of_the_review": "The paper proposes an interesting method to train hybrid models for both identification and prediction quality. \nConnected two ideas from the literature in a single framework is also a welcome contribution. The proposed method is supported by experiments on different domains. Theoretical aspects are promising but still limited.", "strength_and_weaknesses": "", "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": "RIGnAUWMEF_", "reviewer": "Reviewer_1AX7", "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": "This paper is very well written and the problem being addressed is very relevant to the machine learning community, especially to its emerging applications in science and engineering problems. The proposed method is sound and quite easy to implement, and the performance is also good as presented. I think it will be of interest to the broad readership on the intersection of machine learning and scientific applications.\n\nBelow, I have some questions that the authors may wish to clarify.\n\n1. The overall approach: The key method introduced to ensure well-posedness is page 4 Eq (7), via constraining each part of the model by $\\mu_u,\\mu_k$. However, to me it is often more convenient in practice to consider the “dual” formulation of optimization with penalization, i.e. instead of (7) one can consider\n   $$\n   \\min_{h_u,h_k} d(f, h_u+h_k) + \\lambda_k \\ell(h_k) + \\lambda_u d(h_u,0)\n   $$\n   This is certainly very easy to implement, and has the same number of hyper-parameters as the proposed case, replacing the $\\mu$‘s by $\\lambda$’s. Can the authors comment on why the current constrained approach is chosen over this approach, and can some quantitative comparisons be made?\n\n2. Algorithm 1: while successive projections to convex sets can find the intersection eventually, the convergence rate may be slow. Partial projections can sometimes speed up the process. Do you observe speedups if equation (8) is not solved to completion in each step?\n\n3. Prop 3: $D$ is proved to converge (and thus $A$) via a contraction argument. Can anything be said about the convergence point? Can it be a local minimum but not global?\n\n4. Table 1: As I understand, the performance of the proposed method will depend on the choice of the hyper-parameters $\\mu_k,\\mu_u$. Can the authors comment on the sensitivity of tuning parameters with respect to the good results shown in this table? How is fair comparison ensured, e.g. how are the baseline method(s) (e.g. Yin et al) tuned to each problem?\n\n5. Section 4.2: This example is interesting as it is more realistic for scientific applications. I believe this part has some clarity issues in the main text\n\n   1. I believe that Eq (32-33) in the appendix should be included in the main text to allow for easier understanding of what is $f_k$ and $f_u$.\n   2. Looking at Alg 3 in the appendix, the alternate minimization is in fact done by gradient descent with soft constraints instead of the alternative projections introduced in Alg 1 in the main paper. Can the authors explain why this modification is made, and whether the prior experiments on toy examples followed Alg 1 or Alg 3?", "summary_of_the_paper": "This paper studies the identification and prediction of dynamical systems from observed data by combining statistical approaches with physical priors. The main idea is to decompose the approximation of unknown dynamics into two parts: a parametric physics-based model $h_k$ (used to incorporate prior knowledge, or to identify parameters), and a non-parametric statistical model $h_u$ (used to fit the unknown, or unobserved degrees of freedom). \n\nSince the statistical part has a priori no restrictions, the overall problem is ill-posed in the sense that the statistical part can fit all of the dynamics. To ensure well-posedness, the authors propose to convert the joint learning problem to learning $h_k$ and $h_u$ subject to constraints: $h_k$ is given an allowed error tolerance and $h_u$ is given an allowed deviation from 0. The authors show that this modification implies well-posedness, and subsequently design an alternative minimization algorithms to solve for $h_k$ and $h_u$ iteratively. It is proved that this algorithm converges in the affine setting (via sequential projection to find intersection of convex sets) and it is demonstrated experimentally in low and high dimensional systems that this method works well in practice outside of the convex setting.", "summary_of_the_review": "This is a good paper tackling a relevant problem. The algorithm introduced is easy to implement and has promising performance. Some theoretical justification on convergence is also presented. I believe that this is a meaningful contribution and should be accepted. There are some clarity issues that the authors can address in a revision.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and 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": "OnjXFZGpjA", "reviewer": "Reviewer_TpQa", "summary": "", "questions": "", "limitations": "", "rating": 6, "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": "Authors propose a unifying framework for hybrid models, generalizing recent methods. They derive theoretical guarantees very compelling theoretical guarantees for the framework, as well as a convergence proof for a simple affine model. However, I believe that the paper is quite hard to read, and propose several models without giving intuition for which model to use. The models also seem to need complex hyper-parameter choice, which again is not very intuitive and not discussed extensively in the supplementary material. \n\nHere are two points that I think need to be discussed more extensively in the paper:\n - The paper proposes two models (with and without auxiliary data). When should we use one or the other? Should we just use it whenever auxiliary data is available? Moreover, equation 5 is the equivalent of equation 4 for auxiliary data. However, the model uses alternate optimization to converge. What is the equivalent procedure for auxiliary data? \n - The \"practical optimization\" is not straightforward, and although it is well-detailed, the hyper-parameters chosen for the experiments are pretty complex, which is to me an important drawback of the framework. To give intuition into the model, how do the lambda values relate to the two scalars mu in equation 6? \n \nAlthough the section titled \"THEORETICAL ANALYSIS FOR A LINEAR APPROXIMATION\" is compelling, I think it could be included in the appendix to allow for more space to discuss the previous points. \n\nOther points that would make the paper easier to read : \n - Using \"Ours eq 4\" in the results is not intuitive at all, since the model optimizes a different objective. This makes it a bit confusing. \n - Figures are not rendered properly, many of them appear blurry and could be made bigger (there are lots of white space in the figure that could be removed). This greatly affects the quality of the paper. \n - I would like to see the differences between true and estimated processes in the figures. This would make the results easier to grasp.\n\nThe method is compared to Aphynity and Neural ODE but there are no comparisons to Neural ODE on the simulated datasets (I think it would be beneficial to show the predictive logMSE for Neural ODE in Table 1).\n\n", "summary_of_the_paper": "Authors propose a general framework for learning and ensuring identifiability with hybrid models. They derive strong theoretical guarantees as well as a proof of convergence in a simple affine setting. They validate their framework in several simulated experiments and the Natl dataset, and draw comparisons with state-of-the-art hybrid model Aphynity as well as Neural ODEs. Experiments confirm that the proposed framework increase interpretability and maintain high prediction performances", "summary_of_the_review": "Authors propose a unifying framework for hybrid models, generalizing recent methods. They derive theoretical guarantees very compelling theoretical guarantees for the framework, as well as a convergence proof for a simple affine model. However, I believe that the paper is quite hard to read, and propose several models without giving intuition for which model to use. The models also seem to need complex hyper-parameter choice, which again is not very intuitive and not discussed extensively in the supplementary material. \n\nThus, I do not recommend acceptance of this paper, but I am willing to increase my score if authors answer my concerns, as the theoretical guarantees derived here are very strong and surely beneficial for a large class of hybrid models.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and 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": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "MCNzADzNK8q", "reviewer": "Reviewer_GA76", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "Strength: The main strength of this paper is the experiments. In the experiment, the results is promising, with better performance in both synthetic and real-world experiments. Ablation studied are also performed to show the model performance with different objective functions.\n\nWeak: I think the weakness of this paper is the technical novelty of the algorithm and theorem. The idea of combining physical models and machine learning models exist for decades. The theoretical derivations in this paper is nothing but applying triangle inequalities here and there. The alternate minimization algorithm is not novel at all, and the convergence of such an algorithm for a linear dynamic system is also well-known. Therefore, the technical contribution of this paper seems trivial.", "summary_of_the_paper": "This paper investigate a novel method to recover well-posedness and interpretability for the problem of learning a physical and machine-learning hybrid model, by controlling the ML component and the MB hypothesis. The authors also study the case where auxiliary data are introduced. Experiments show the performance of the proposed method in both synthetic data and real-world ocean dynamics.", "summary_of_the_review": "The experiment results of this paper is very strong and comprehensive. However, the novelty of this paper is not very good. The motivation, theoretical analysis and algorithms are not very novel.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and 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:

Meta-review: The paper proposes a method for hybrid model-based/ML learning, where a model is decomposed into an interpretable parametric prior and a neural net residual.  In this case, the prediction error minimization does no identify the parametric component, and an alternating optimization method is proposed to augments prediction error loss with component-specific losses. Empirical and theoretical results are obtained.  Initial questions of several reviewers were addressed.
Decision: Accept (Poster)