Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 7: Good paper, accept
review_content: summary: This paper presents a novel approach to model and capture the uncertainty of dynamical systems by learning from trajectory data. The method is composed of a Gaussian process smoother model, used for predictions, and a neural dynamics model, used for training. The novelty in the method comes from the use of the distributions over time-derivatives of the state vectors (via a Wasserstein-distance penalty between the distributions from the smoother and the dynamics model) during the training process. Basing the approach on integrals over the state space, instead of over the model parameters space, allows for computational efficiency and more robustness when compared to previous sampling-based approaches, as evidenced by experiments.
main_review: The paper is mostly well written, providing a good motivation, an appropriate use of citation and interesting insights in the derivation of the method. However, there are a few issues, especially regarding methodological details, which make it hard to understand the method and to assess its soundness.

Major issues:
- The description of the methodology is somewhat confusing in a few important aspects.
For instance, how does the set of "supporting gradients" play a role in the estimation of these two distributions? Are they random variables to be inferred as part of $\dot{\mathcal{X}}$? Or are they additional observations which will be used somehow in the calculation of the losses?
- Regarding $\mathbf{z}$ in Eq. 6, I get the idea, but it's a bit confusing at first-reading to use both $\mathbf{z}$ and the $(\mathbf{x}, t)$-tuples as inputs to the GP when the feature map "phi" is supposed to be part of the kernel. It'd be better to consider $k := k \circ \phi$ (composition), since the composition of a kernel with an input map is still a kernel, and then represent the inputs as simply $(\mathbf{x},t)$, without the extra $\mathbf{z}$ notation.
- In line 196, it is mentioned that the (instant) state noise vector in Eq. 12 is drawn only once per rollout. Does it mean that the noise term epsilon is drawn once and repeated for all time points along a trajectory? If so, the state noise is not independent, but actually correlated across time and that would possibly not allow the model in Eq. 12 to properly capture uncertainty in the state transitions.
- I'm not sure about the effects of this approximation in (14), since later on the difference between $p_D$ and $p_S$ will be minimised. Having $p_D$ as a function of $p_S$ might have negative side effects. Any ideas on the drawbacks?
- The certainty-equivalence approximation (see line 220) discards the uncertainty captured by the GP model in the expectation in Eq. 15. Why not using a quadrature approximation to compute Eq. 15, like the unscented transform, or anything else which could capture the states covariance matrix on the inputs?
- In the experiments with SGLD and SGHMC, how is the likelihood for these sampling-based methods computed? These methods don't directly produce a probability density, only empirical approximations.
- The experimental evaluation could be more complete if they included comparisons against methods which encode (inexact) prior knowledge about the dynamics, such as simulator-based inference approaches [e.g., A, B, C, below]. The relationship to these methods is also not discussed in the related work section. In my view, the proposed approach could be combined with simulator-based models by using a physics-informed mean function for the GP smoother, which is not explored/discussed in the paper.

Minor issues:
- Parametric vs. Non-parametric: Neural networks are "parametric" models, though in the main paper and in the appendix they are referred to as "non-parametric" models. The only non-parametric model in this work is the GP.
- The equation for $\mathbf{z}_i$ in line 164 is missing $\phi$.
- I couldn't find details of this weight sharing scheme in Sec. 4, despite the mention in line 200.
- Preliminary experiments: Are these results available? If not, it'd be better to add them to the appendix or give more details/references to back it up. For example, the Wasserstein GANs paper provides some experiments and theoretical justification on why the Wasserstein distance is more appropriate than the KL divergence for models over high-dimensional space, which usually have their data distribution concentrated on a (latent) lower-dimensional subspace.
- Experiments: At the first mention of SGLD, SGHMC and NUTS, please spell out and/or add references to these methods for readers who are unfamiliar with the literature. Also, please, consider adding a reference to NDP in line 265.
- Experiments: Line 306, any idea why SGLD failed with the Quadrocopter 1?
- Please, add (foot)note on what is an "Ansatz". At first, I thought it was a typo, as I was unfamiliar with the term.

References:

[A] Cranmer, Kyle, Johann Brehmer, and Gilles Louppe. 2020. “The Frontier of Simulation-Based Inference.” Proceedings of the National Academy of Sciences.

[B] Ramos, Fabio, Rafael Carvalhaes Possas, and Dieter Fox. 2019. “BayesSim : Adaptive Domain Randomization via Probabilistic Inference for Robotics Simulators.” In Robotics: Science and Systems (RSS). Freiburg im Breisgau, Germany.

[C] Okada, Masashi, and Tadahiro Taniguchi. 2019. “Variational Inference MPC for Bayesian Model-Based Reinforcement Learning.” In 3rd Conference on Robot Learning (CoRL 2019).
limitations_and_societal_impact: Some of the scalability limitations are explored in the appendix, and no negative societal impact is foreseen by the authors.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: This paper develops a method called distributional gradient matching for learning unknown differential equations from data. The method is based on a combination of a neural ODE model and a Gaussian process model (with a deep covariance kernel) which is encouraged to produce solutions which follow the neural ODE model by introducing a regularisation term based on Wasserstein loss. Numerical examples show that the proposed method outperforms some existing alternatives.
main_review: I am not an expert on learning ODEs, but the literature reviews in Sections 1 and 5 do make it seem like the approach is new. It would be interesting to have more commentary on how the use of the GP model in this paper differs from prior work - I presume that most prior work does not use a deep GP model.

The numerical examples seem impressive, the quadrocopter system in particular being highly non-trivial. In these examples the proposed method clearly outperforms some existing methods, both in terms of accuracy and speed. I lack expertise to say if the two MC-based algorithms, SGLD and SGHMC, that the comparisons are made to are representative of the state-of-the-art in this setting.

In my opinion the main weakness of this paper is that the presentation of the method is often quite confusing and the notation not always properly explained:
- Line 141: What does \dot{\mathcal{X}} stand for?
- p4: It is not told what the GP model is supposed to be modelling, just that "we define a Gaussian process".
- p4: The use of a mapping \phi is not standard in GP regression (at least yet), so it would be good to make it clearer that this part of the GP model comes from Wilson et al. (2016).
- p4: To call \phi a "feature map" seems potentially confusing if one is used to defining covariance kernels via inner products of feature maps.
- p4: Is the GP covariance kernel supposed to be potentially different for each k (and hence the notation \mathcal{K}_k)? If so, is it intentional that there is only a single prior mean function?
- Remarks on p6 that it has "trainable parameters" and in Appendix B.2 that it is "deep" appear to be the extent of exposition on the GP prior mean function \mu. More should perhaps be said.
- What are t and t_{supp} in the definitions on p4?
- \mathcal{L}_{data} defined in Equation (8) is never used afterwards; it is denited p_{GP}(y | \varphi) later.
- Is \mu_{S,supp} on p6 equal to \mu_S in Equation (7)?
- It would be helpful to tell explicitly that the Wasserstein loss in Equation (16) will be used as \mathcal{L}_{dynamics}.

I am not sure if the paper contains sufficient details to replicate the experiments.

Other minor comments:

- The sentence on lines 119-121 is quite awkward. There is also an extra closing parenthesis on line 120.
- Line 131: "smother-dependent"
- I believe the dot over \mu in Equation (7) is larger than that in the equation between lines 166 and 167.
- The GP prior mean needs to be subtracted from second term in Equation (10).
- Transpose notation is not consistent: E.g., Equations (10) and (11) use "T" while (6) and (7) use "\top".
- \mathcal{D} in Equation (9) should probably have subscript k.
- That the noise variance matrix is \Sigma_\epsilon given that Equation (12) has a different random variable \epsilon.
- Line 196: \epsilon should probably be bold.
- Lines 216-17: "regularizaton"
- Line 223: "indepdent"
- "Lahdesmaki" - "Lähdesmäki" in one of the references

limitations_and_societal_impact: Limitations caused by the cubic computational cost of GP regression are mentioned.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: The paper tackles the problem of learning nonlinear dynamical systems from historical trajectories, while estimating predictive uncertainties. They propose a new method based on a Gaussian process smoother which is trained with maximum likelihood to fit the data of all training trajectories jointly, while constrained to follow in some sense the underlying dynamics. The latter is performed thanks to a regularizer forcing the smoother’s gradients distribution to match the  distribution of gradients of the dynamics model over a fixed set of support points. The Wasserstein-2 distance is approximated to evaluate the distribution discrepancies. Then, instead of deriving uncertainty estimates over the predicted parameters, as done in standard Bayesian methodology, uncertainty is directly estimated over the predicted states. For this, a Gaussian distribution is used to model  probability of the states derivatives conditioned on the states. Thanks to these tricks and modeling choices, the authors avoid to have to integrate the complex dynamics model, unlike previous Monte Carlo methods. The proposed approach is show-cased and compared to sampling-based methods in several experiments with simulated data, including overparametrized linear dynamics, as well as single and multi-trajectory cases of 4 well-known dynamical systems. An ablation study is also presented showing the benefits of joint smoothing and gradients distribution matching regularization on one of the four simulated datasets.
main_review: Post-rebuttal additions
------------------------------

I would like to thank the authors for their detailed answers, which address all concerns I had.

-----------------------------------------------------------------------------------

Originality
--------------

The paper seems original to me and the related work seems quite well referenced to me.

Significance
----------------

The proposed approach is very complex, but not more than existing ones, so I believe it could be a good alternative to continuous-time nonlinear system identification requiring uncertainty estimates. It solves indeed a real bottleneck of existing approches when used with NODE models (numerical integration). Also, instead of delivering parameter uncertainties as previous methods, is focuses on direct estimation of state uncertainties, which is what matters at the end of the day. Hence I believe that the work present could be quite significant in the field, both for practitioners and for researchers who could build on these ideas.

Quality
----------

Overall quality of the explanations and experiments is quite good, as the method was show cased in varied settings and compared to relevant state-of-the-art alternatives. Also, the code seems clean (+1 for reproducibility). I do have some remarks though concerning some points that do not seem explained in the paper, as well as a few questions listed below:

1. While it is said in section 2 that trajectories can be sampled at non-equality spaced times, all simulated experiments use uniform time grids, which is a shame.
2. How are the supporting points chosen and how do they influence the estimates’ quality? For example, having a set of supporting points not well spread or too small should lead to degraded performance should it not?
3. In equation 12 and line 196, it is said that dynamics noise $epsilon$ is drawn one per trajectory, while in the problem statement, observations $y$ are supposed to have i.i.d additive noise (one per sample in each trajectory). While I realize these quantities are not the same (one is in the states and the other in their derivatives), I wonder why make these different modeling choices and not also have a different disturbance realization for each instant in the dynamics model.
4. Concerning the ablation study, you say on line 330 that by comparing the rows of Figure 1 it can be seen that joint smoothing of multiple trajectories helps. While I do see it clearly when comparing rows of the second column (i.e. with the dynamics regularizer), I see the opposite when looking at the first column. Namely, joint smoothing seems to lead a noisier and less accurate predicted trajectory, with a larger confidence region. As this is very curious, since it contradict the text, I think it should be commented and also that the sentence should be clarified (by saying that improvement is visible on the second column for example).
5. How many training and testing trajectories are used to compute results for figure 2?
6. Likewise, while it is said that 5 to 10 points per trajectory are used to train the models in the multi-trajectory experiment, how many trajectories are used? Is that the meaning of the numbers next to the systems considered in the first column of Table 2? If that’s the case, please add it to the legend or paragraph text.
7. Could you please comment the fact that SGLD and SGHMC methods outperform DGM for configurations 3.3 and 3.6.3 on figure 2?
8. Could you please develop on the reason SGLD failed on the quadcopter single trajectory experiment, leading to NaNs? Did you try to understand what these NaNs meant? 

Clarity
---------

1. Although the paper is well written, the proposed method is very complicated and difficult to grasp (I had to read the paper many times before I started to get it). I cannot say it is not well explained though, but it might help to put more structure in section 3 (replace paragraphs by subsections + parapgraphs for example). Also, adding a diagram summarizing how the many elements of the method interact with each other could help a lot I believe.
2. Furthermore, I found some explanation was missing concerning the certainty-equivalence approximation of the Wasserstein distance. Maybe add that the expectation you are approximating is the one in (15) and that this approximation is carried in the second line go (16).
3. In the paragraph Final loss (page 6), the first data fitting term of the loss is denoted quite differently than in equation (8). It took me some time to realize they were the same, so I think some referencing and more consistent notation could help with the clarity here.
4. In line 285, you denote the model used to fit the dynamics with the same notation as in line 282, where you explain the real dynamics that generated the data. I find this confusing and believe that less ambiguous notation here would help with the clarity (by using for example the dynamics function $f$ introduced in (1)).
5. I could not find any introduction of the abbreviations SGLD and SGHMC.
6. On line 164, I believe $z_i:=$ should not be there, as it is defined differently in the text and in equations below line 166.
7. Sentence at lines 119-120 has a problem as the word maps is repeated twice.
8. Line 244: another -> each other
limitations_and_societal_impact: I can't see any societal impact related to this publication.

Also, the authors don't seem to have addressed the limitations of the proposed approach.
review_title: null

final_score: rating: 5: Marginally below the acceptance threshold
review_content: summary: The paper proposes to skip the integration of continuous-time models and directly learn a model that predicts the posterior of the state conditioned on the initial state and the desired time. To learn such a model with 'good' uncertainty estimation the authors mix GP's, neural networks and the wasserstein metric (I did not understand how the authors combine these approaches together). Within the experiments the paper shows that the log-likelihood is better for the proposed model on uncontrolled toy systems, e.g. Lotka Volterra, chaotic Lorenz, double pendulum and quadcopter.

main_review: First of all, I am an expert on learning dynamics models but not an expert for Bayesian models. My main problem with the paper is that I cannot understand the proposed approach. The authors described the problem and their approach. However, once the authors mix all concepts of bayesian networks, GP's, Wasserstein distance and stir it all together. I personally do not have any understanding how all of this plays together. Therefore, the paper needs a major revision that restructures the writing such that the reader gets a good understanding of the approach. A good approach is to provide the intuition of the proposed approach and only afterwards deep dive into the details. However, the authors skip the intuitive approach and directly deep dive into the theory and lose the reader in the process. Besides the writing, the performed experiments are not really impressive and only try a few toy domains. Even for the toy domains the experimental evaluation does not provide much insights except my number is bigger than yours. It would be beneficial to incorporate more qualitative evaluations and explain to the reader what to look out for. Furthermore, the authors motivate their approach using control and RL but then only model uncontrolled dynamics and never scale the proposed approach to interesting control problems. How does the model work for high-dimensional systems with contacts? Therefore, it would be necessary to evaluate the proposed approach on more complex domains such as walker or humanoid.

As the writing and structure of the paper needs a major revision and the performed experiments do not evaluate whether this model actually improves learned models for control, I cannot recommend the paper for acceptance.
limitations_and_societal_impact: N/A
review_title: null

Meta-review: Even if there was some spread in the review scores which called for discussion among the reviewers, in the end, the reviewer consensus was in favour of accepting this paper. Please, make sure to address their concerns in the reviews in the camera-ready. The concerns related to clarity, presentation, and notation (which you also comment on in your response) are of particular importance as these help ensure impact of your work.

Decision: Accept (Poster)