arXiv:1611.05843v1 [math.ST] 17 Nov 2016

Bayesian Analysis (2016)

TBA, Number TBA, pp. 1–2

Contributed Discussion to Bayesian Solution
Uncertainty Quantification for Differential
Equations by Oksana A. Chkrebtii, David A.
Campbell, Ben Calderhead, and Mark A.
Girolami
William Weimin Yoo∗
Abstract. We begin by introducing the main ideas of the paper, and we give a
brief description of the method proposed. Next, we discuss an alternative approach
based on B-spline expansion, and lastly we make some comments on the method’s
convergence rate.
Keywords: Differential equation, discretization uncertainty, B-splines, tensor
product B-splines, convergence rate.

I would like to congratulate the authors for such an interesting research. The Bayesian
method with the probabilistic solver introduced is highly innovative and practical. The
various examples presented in the paper show the wide applicability of the proposed
method. However, I do find the title a bit of a misnomer, as I initially thought that
the authors are constructing credible sets for the fixed but unknown solution u∗ of the
differential equation.
The inverse problem that the authors are trying to solve, in its most basic form is
this: Suppose you have observations Y = Au + ε, where ε is some normal errors and u
follows ut = f (t, u, θ). Here, A is a known transformation from the state space u to the
observation space Y , ut is the first order derivative with respect to its argument t, f is
the known form of the differential equation, and θ’s are the equation’s parameters. The
method proposed consists of two steps, with one nested within the other. First, solve
for u probabilistically to obtain a discretized solution at some grid points. Then we
embed these discretized version of u in a Bayesian hierarchical framework to estimate
θ. To model discretization uncertainty associated with using only u evaluated at grid
points, the authors endow priors based on Gaussian process jointly on u and ut , where
the covariance function is constructed by convolving kernels.
There is an alternative and perhaps a conceptually easier way P
to achieve the same
J
result. We can first represent u by a B-spline series, i.e., u(t) = j=1 ϑj Bj,q (t) with
Bj,q (·) denoting the jth B-spline of order q, and we endow the coefficients ϑj ’s with
normal priors. Here, the number of basis J plays the role of 1/λ, where λ is the lengthscale parameter defined in the paper. It turns out that the first derivative of this u is
PJ−1 (1)
(1)
is some weighted first
another B-spline series ut (t) =
j=1 ϑj Bj,q−1 (t) where ϑj
∗ Mathematical

Institute, Leiden University, The Netherlands yooweimin0203@gmail.com

c 2016 International Society for Bayesian Analysis

DOI: 0000

imsart-ba ver. 2014/10/16 file: BA1040.tex date: March 12, 2018

2
order finite difference of the ϑj ’s ((4.23) of Schumaker (2007)). Therefore, u and ut
are jointly normal and their associated covariance matrices are banded due the support
separation property of B-splines. To enforce the given initial condition, we can condition
the joint prior (u, ut ) on u∗ (0).
Moreover, this approach can be generalized to the partial differential equation case,
where we take tensor product
P ofPB-splines to model both the spatial and temporal
components, i.e., u(x, t) = Jj1 =1 Jj2 =1 θj1 ,j2 Bj1 ,q (x)Bj2 ,q (t). As in the univariate case,
partial derivatives of tensor product B-splines will be another tensor-product B-splines
((3.2) of Yoo and Ghosal (2016)). Hence we will obtain the same Gaussian process prior
for u and all its mixed partial derivatives if we endow normal priors on the coefficients.
As before, we enjoy some simplification in computing the covariance matrices because
they are banded.
In addition, I would like to comment on the effect of grid point distribution on
the convergence rate of the proposed algorithm. Intuitively, one would expect that the
grid points should be chosen roughly uniformly across the domain [0, L]. Suppose we
choose grid points {t1 , t2 , . . . , tN } and we further assume that they are quasi-uniform,
i.e., h/ mini (ti − ti−1 ) ≤ C for some constant C > 0 with h = maxi (ti − ti−1 ). In other
words, the max grid increment is of the same order as the min grid increment. Then
it follows that h is of the order of 1/N and by Theorem 1 of Chkrebtii et al. (2016),
the rate of convergence is O(N −1 ). Therefore for quasi-uniform grids (which includes
uniform discrete grids), increasing the number of grid points will result in more accurate
solution.
The paper under discussion Chkrebtii et al. (2016) makes significant contribution to
the new field of probabilistic numerics. I have learnt a great deal by reading this paper,
and I hope that there will be more papers in uncertainty quantification for differential
equation models in the future.

References
Chkrebtii, O. A., Campbell, D. A., Calderhead, B., and Girolami, M. A. (2016).
“Bayesian solution uncertainty quantification for differential equations.” Bayesian
Anal., 29 pages. Advance publication.
Schumaker, L. (2007). Spline Functions: Basic Theory. Cambridge University Press,
New York, third edition.
Yoo, W. W. and Ghosal, S. (2016). “Supremum norm posterior contraction and credible
sets for nonparametric multivariate regression.” Ann. Statist., 44(3): 1069–1102.

imsart-ba ver. 2014/10/16 file: BA1040.tex date: March 12, 2018

