Venue: NeurIPS 2022 Conference

Reviews:

final_score: rating: 7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
review_content: summary: The paper presents a way to algorithmically construct Gaussian Process
covariance functions in such a way, that realizations strictly follow a given
linear ordinary differential equation. The proposed construction via smith
normal form algorithms overcomes existing limitations on controllability of the
system, which exists in other approaches. The resulting LODE-GPs are compared to
classic GPs in three experiments. 

strengths_and_weaknesses: # Strengths
The presented method appears to be novel and I believe that it is a significant
contribution to research on Gaussian processes and generally machine learning in
the context of dynamical systems. The technical statements in the paper appear
sound. Related work is adequately cited (some comments below). Additionally, the
paper is very well-written and the method is presented in a clear way; the
bipendulum example used throughout the paper is very helpful, and the
experimental setup is discussed in much detail in the appendix.

# Weaknesses
The main criticism I have for the paper is about the comparison only to a
classic Gaussian process in the experiments. While I believe that this _might_ be
appropriate, it would be helpful to comment more extensively to why this is the
only reasonable comparison, and to define the problem setting more formally /
precisely.

For example, the experiment settings in section 5 can also be interpreted as
parameter / initial value inference problems, which could be approached with
e.g. least-squares regression with classic ODE solvers (e.g. [Biegler]) or
probabilistic numerics [Tronarp], or (probabilistic) gradient matching (e.g.
[Wenk]). Note that the former two approaches would further come with the
advantage of having linear O(n) scaling (due to not modeling covariances or to
using Gauss--Markov priors, respectively). Clarfying the problem setting further
could be helpful for readers to understand why these do not apply; otherwise
some of these approaches might be relevant baselines.

[Biegler] Nonlinear parameter estimation: A case study comparison, Biegler et al, 1986  
[Wenk] ODIN: ODE-informed regression for parameter and state inference in time-continuous dynamical systems, Wenk et al, 2020  
[Tronarp] Fenrir: Physics-Enhanced Regression for Initial Value Problems, Tronarp et al, 2022

questions: My most relevant question / suggestion is about the relation to parameter inference problems, stated above under "weaknesses".

Some questions / comments on probabilistic numerics (PN):
- l. 23 writes that probabilistic ODE solvers allow to "learn" non-linear ODEs,
  but it is not quite clear to me in which sense this holds (i.e. inferring the
  solution or inferring parameters). To my understanding, those methods "solve"
  or "estimate the solution of" ODE initial value problems (IVPs)
- Solving IVPs with non-linear ODEs _requires_ approximations, since they can, in
  general, not be solved in closed form. It might thus be misleading to state
  approximations as limitations of the solver, rather than a necessity of the
  general problem. In contrast, IVPs with linear ODEs can be solved in closed
  form, and therefore PN does not apply since no numerical error needs to be
  estimated.
- l. 70: What is meant with "complex calculations"? It could be helpful to
  clarify.
- Page 6, footnote: What is meant with "free functions in their solution set"?
  Depending on the meant object, [47] might possibly be in disagreement since it
  describes inference of time-varying system parameters (= unknown functions).

I want to highlight that these comments are only meant as helpful suggestions and not as necessary changes.
The resolution about the relation to parameter inference, commented on above, as
well as possibly the extension of the baseline for the experimental evaluation
are of higher relevance.

Minor comments:
- l. 18: "such _as_ in"
- l. 88 contains a lengthscale l, but it is not mentioned.
- Around Lemma 1 it would be helpful to point the reader to the proof in the
  appendix.

limitations: The authors adequately discussed the limitations of the proposed method.
review_title: null

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: The paper describes how to constrain realisations of a Gaussian process $f \sim GP(m, k)$ to satisfy $A f = 0$, where $A$ is a linear, matrix-valued differential operator with constant coefficients.
This is done by computing the Smith normal form of $A$, which is a (symbolic) decomposition of $A$ into
$U D V = A$, or equivalently $D = U^{-1} A V^{-1}$, and identifying the null space of $D$.

This is most closely related to the work by Hegermann (2018) and Hegermann (2021), who tackle a similar problem using a Gröbner basis.
The difference between the submitted work and the aforementioned papers is the use of a Smith normal form, which allows for non-controllable systems, and can be implemented with standard software. The Smith normal form also restricts the class of problems that can be tackled to linear, homogenous, ordinary differential equations with constant coefficients.
strengths_and_weaknesses: ## Strengths
Overall, the methodology seems to be sound, the related work seems to be discussed adequately, and the experiments support the claims.
I would expect that the method described in the paper will be useful to some people. The paper is at the interface of differential equation models and Gaussian processes, which seems to be a relevant subject for Neurips.


## Weaknesses
* One criticism of the method is that the algorithm is restricted to linear, homogenous ODEs with constant coefficients and that there may not be much hope to generalise to nonlinear/nonhomogenous problems or even partial differential equations.
* From a didactical viewpoint, I would appreciate it if the paper could provide more intuition for Lemma 1. Maybe a few equations would help.
* This is minor, but I would like to see samples from the constrained Gaussian process, perhaps in one of Figures 2, 4, 5, or 6?
* Why does only Figure 2 show the credibility intervals? The remaining figures only discuss the mean. Is there a specific reason?
questions: 
* Abstract: what does "novel algorithmic and symbolic construction" mean? It is written frequently, but I don't understand "algorithmic and symbolic construction".
* Line 87: Please explain the following notation in the paper $\mathbb{R}_\preceq$. I suppose it refers to some sort of positive-definiteness?
* Line 223f: "We additionally analyse the training runtime and see that...": are those results in some figure/table? I can't seem to find the results

limitations: It would be useful if the limitations of the method would be discussed more centrally. What do I mean by this? Line 68 states "but are restricted to ODEs". Lines 72f. describe the limitations as "restricted to linear ODEs". Line 140 states: "Consider a system of linear homogenous ODEs with constant coefficients". Line 172 says: "This is generally impossible for PDEs". Line 255 states: "follow a given system of linear homogeneous ODEs".
It feels like a reader has to search the entire document to learn what the limitations of the method are.
review_title: null

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: The authors introduce a new method which tackles the important and much discussed problem of learning the solutions to systems of ordinary differential equations from noisy data. The authors construct a Gaussian process prior for the solutions to the equation, which has been a common approach in previous work in the area. In particular the authors attempt to solve the problem in the case of linear differential operators for systems which are underspecified and uncontrollable. The proposed method relies on the symbolic calculation of the Smith normal form of the differential operator representing the system, which the authors use to construct a covariance over both the system outputs and inputs, ensuring that samples obey the differential equation in question. The proposed method is validated on 3 synthetically constructed physical systems, providing convincing evidence that the method works.

strengths_and_weaknesses: Strengths 
---------

Overall I enjoyed reading the paper, and think it has a number of strengths. The paper for the most part is very clearly written, and introduces the problem, and the pre existing solutions, well. I particularly like the use of explicit examples at key points in the exposition which help the reader to understand how the proposed method is applied in practice. The figures are mostly clear and useful. 

The idea presented in the paper is quite an elegant one in my opinion, and seems to be a sensible way to tackle the problem. The fact that the method can be used to construct priors which consist of exact solutions to the given system is compelling. 

The method seems to work on a number of different problems, and in some cases shows good performance relative to standard GPs, although I do have some concerns about the evaluation, as listed below. 

I am not expert in the theoretical tools such as the Smith normal form used in the paper, and so I can not vouch for the correctness of those parts with 100% certainty. That being said, the results given in the paper are correct to the best of my knowledge.

Weaknesses
----------

There are a few weaknesses with the paper, mostly related to the evaluation of, and motivation for, the proposed methods. 

**Motivation for uncontrollability**
The authors say that their methods is an improvement on [34] because it can deal with uncontrollable systems.  As the authors note, in 2 of the experiments which are controllable, the proposed method generates exactly the same results as the previous work of [34]. To me, there is not sufficient discussion of the reasons that the uncontrollable case is an important problem, and thus why the proposed method represents an important contribution. I am not an expert in control theory, so it may be that I am unaware, but I think it would be good to motivate this problem more clearly in the text, and explain why it is important to be able to model such systems. Additionally for the experiment where the proposed model is novel (three tanks), the performance seems to be worse than a vanilla GP in terms of it's representation of the data. 

**Comparison with related work**
Presumably the proposed method can work for problems that are well specified, as well as those discussed in the text that are underspecified. Currently there is no indication of how this method performs in that case, and how well it performs relative to other prior work. I think it would improve the paper considerably to include some comparisons. 

The authors claim that "No other state of the art in probabilistic ML exists for our class of differential equations" I am not entirely sure that this is correct. I believe that is is possible to solve some of the problems using the framework of latent force models (LFMs) [Alvarez et al. 2009] which represent the input functions as GPs, and use the Green function of the differential operator to map to exact outputs. Most work on LFMs deals ODE systems which are not coupled, i.e. the part of $A$ relevant to the outputs is diagonal, so it may be additional effort to get systems 1 & 2 running in the LFM framework, but I believe system 3 can easily be recast in terms of LFMs. System 3 is equivalent to the model in equation 3 of [Alvarez et al. 2009], when $D_q = 0$, $B_q=0$ and $S_{rq}=\\begin{pmatrix}
    -1 & 0\\\\
    -1 & -1 \\\\
    0 & -1
\end{pmatrix}$, the expressions given in later equations can then be used to compute the covariance of the outputs. It would be great to have a comparison between the proposed model and LFMs for system 3, at least i think a comment on this equivalence is necessary. Again I may be mistaken in this point, but I would appreciate a comment from the authors if so. 

**Evaluation**
A minor point, but I'm not sure "ODE error" shown table is interesting to compare with the GP. The metric is useful as a check that the model is working, but it doesn't seem fair to compare with the GP since the information about the differential equation is not available to the GP, and the error should be small for the proposed model almost by definition. 

Table 2 shows results for 10 repeats, with bold text indicating best performance, but does not include the standard deviation or error on these results. This means it is not possible for the reader to assess the significance of the results. Please add sandard errors or deviations to this table

##### References
- [Alvarez et al. 2009] Alvarez, Mauricio, David Luengo, and Neil D. Lawrence. "Latent force models." Artificial Intelligence and Statistics. PMLR, 2009. 
questions: Below are my questions/suggestions for the authors, which are mostly a summary of what is discussed in the weaknesses section. I've numbered the points separately to allow the authors to respond more easily. 

1. Please can the authors comment on the significance of the ability to model noncontrollable systems, and provide some real world motivation? In doing so could the authors expand the discussion of how this work represents a significant contribution over [34]?
2. How does the proposed method compare to existing models of this type from problems that are fully specified?  
3. Can the authors comment on the relationship between this model and LFMs? If they are indeed equivalent, the authors should discuss this in the text, and add a comparison for system 3 if possible?
4. I suggest that the authors add standard errors/deviations to table 2, and indicate if the standard errors are narrow enough to indicate one model being significantly better than the other. 
   
If the above points are addressed I would be more than happy to raise my score. 


Below are some very minor typos/formatting/wording issues:

* Line 23 "allow to"
* Line 193 not sure what "All experiments ran numerically stable and robust using Adam [28]." means, please rephrase.
* Line 233 "The LODE-GP can be constructed using parameters." this is not clear, please rephrase.
* Many of the figures don't have titles/axis labels, so are difficult to read, although mostly the relevant information is in the caption/legend. For figure 2 particularly it is not at all clear which plot is which output, presumably it is f1, f2, u from left to right, but please clarify via titles.
limitations: I do not believe there are any ethical concerns with this work or problems with potential negative societal impact. 

Aside from what I have already mentioned, I believe the authors reasonably fairly represent the limitations of their work, although I do think the systems shown in the evaluation are somewhat similar, and it would be good to get some idea of how these models would perform on real world data, not just synthetic systems.
review_title: null

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: This paper proposes a method, based on Smith normal form, to construct Gaussian processes indexed by R such that they satisfy a specified differential equation. 


strengths_and_weaknesses: 
S1 - The paper is clearly written. 

W1 - While the proof of Lemma 1 is indeed trivial enough to defer to appendix, it may be useful to point out in the main text that (a) the proof is indeed given in appendix, and (b) the covariances are obtained by expanding f in a basis of the null-space of d with i.i.d standard normal coefficients. 

W2 - The entire noise free data is not shown in fig 2. It woudl certainly be interesting to see how well the method extrapolates. 



questions: 

Q1 - I guess \ell_1 should really be \ell_2 in entry (3,2) of A in eq. (1) ? 
limitations: Its fine.
review_title: null

Meta-review: metareview: The reviewers found this paper novel, significant for the community, and well written. All four reviewers recommended accepting the paper. I also appreciated the numerous illustrative examples in the paper. You addressed many of the remaining questions/concerns in your rebuttal and in the author-reviewer discussion. Please, go through the reviews once more for the camera-ready and take these into account.
confidence: Certain
Decision: Accept