A grand challenge with great opportunities facing researchers is to develop a coherent framework that enables scientists to blend conservation laws expressed by differential equations with the vast data sets available in many fields of engineering, science and technology. In particular, this article investigates conservation laws of the form
which model the relationship between two black-box functions u(x) and f(x). Here,
and is a parametric linear operator identified by its parameters
noisy observations
), respectively, the aim is to learn the parameters
and hence the governing equation which best describes the data. Such problems are ubiquitous in science, and in mathematical literature are often referred to as “inverse problems” (see [1, 2]). To provide a unified framework for resolving such problems, this work employs and modifies Gaussian processes (GPs) (see [3, 4]), which is a non-parametric Bayesian machine learning technique. Quoting Diaconis [5], “once we allow that we don’t know f (and u), but do know some things, it becomes natural to take a Bayesian approach”. The Bayesian procedure adopted here, namely Gaussian processes, provides a flexible prior distribution over functions, enjoys analytical tractability, and has a fully probabilistic work-flow that returns robust posterior variance estimates which quantify uncertainty in a natural way. Moreover, Gaussian processes are among a class of methods known as kernel machines (see [6, 7, 8]) and have analogies with regularization approaches (see [9, 10, 11]). However, they are distinguished by their probabilistic viewpoint and their powerful traning procedure. Along exactly the same lines, the methodology outlined in this work is related to and yet fundamentally differentiated from the so-called “meshless” methods (see [12, 13, 14, 15]) for differential equations. Furthermore, differential equations are the cornerstone of a diverse range of applied sciences and engineering fields. However, use within statistics and machine learning, and combination with probabilistic models is less explored. Perhaps the most significant related work in this direction is latent force models [16, 17]. Such models generalize latent variable models [18, 19, 20] using differential equations. In sharp contrast to latent force models, this work bypasses the need for solving differential equations either analytically or numerically by placing the Gaussian process prior on u(x) rather than on f(x). Additionally, equation 1 can be further motivated by considering the familiar cases where
for some kernel . This particular instance corresponds to the well-known convolved Gaussian processes [21, 22] which are suitable for multi-output purposes [23]. Moreover, yet another special and interesting case arises by assuming
being the Kronecker delta, which yields
This results in the so-called recursive co-kriging [24, 25] model v(x). Here, v(x) is a latent function and is included to capture unexplained factors. Recursive co-kriging models can be employed to create a platform for blending multiple information sources of variable fidelity, e.g., experimental data, high-fidelity numerical simulations, expert opinion, etc. The main assumption in multi-fidelity settings is that the data on the high-fidelity function f(x) are scarce since they are generated by an accurate but costly process, while the data on the low fidelity function u(x) are less accurate, cheap to generate, and hence abundant.
The proposed data-driven algorithm for learning general parametric linear equations of the form 1 employs Gaussian process priors that are tailored to the corresponding differential operators.
2.1. Prior
Specifically, the algorithm starts by making the assumption that u(x) is Gaussian process [3] with mean 0 and covariance function
where denotes the hyper-parameters of the kernel
. The key observation to make is that any linear transformation of a Gaussian process such as differentiation and integration is still a Gaussian process. Consequently,
with the following fundamental relationship between the kernels
Moreover, the covariance between ), and similarly the one between
), are given by
), respectively. The main contribution of this work can be best recognized by noticing how the parameters
of the operator
are turned into hyper-paramters of the kernels
Without loss of generality, all Gaussian process priors used in this work are assumed to have a squared exponential covariance function, i.e.,
where is a variance parameter, x is a D-dimensional vector that includes spatial or temporal coordinates, and
. Moreover, anisotropy across input dimensions is handled by Automatic Relevance Determination (ARD) weights
. From a theoretical point of view, each kernel gives rise to a Reproducing Kernel Hilbert Space [26, 27, 28] that defines a class of functions that can be represented by this kernel. In particular, the squared exponential covariance function chosen above, implies smooth approximations. More complex function classes can be accommodated by appropriately choosing kernels.
2.2. Training
The hyper-parameters and more importantly the parameters
linear operator
can be trained by employing a Quasi-Newton optimizer L-BFGS [29] to minimize the negative log marginal likelihood [3]
Here, are included to capture noise in the data. The implicit underlying assumption is that
). Moreover,
are assumed to be independent. It is worth mentioning that the marginal likelihood does not simply favor the models that fit the training data best. In fact, it induces an automatic trade-off between data-fit and model complexity. This effect is called Occam’s razor [30] after William of Occam 1285–1349 who encouraged simplicity in explanations by the principle: “plurality should not be assumed without necessity”. The flexible training procedure outlined above distinguishes Gaussian processes and consequently this work from other kernel-based methods. The most computationally intensive part of training is associated with inverting dense covariance matrices K. This scales cubically with the number of training data in y. Although this scaling is a well-known limitation of Gaussian process regression, it must be emphasize that it has been effectively addressed by the recent works of [31, 32].
2.3. Prediction
Having trained the model, one can predict the values u(x) and f(x) at a new test point x by writing the posterior distributions
with
u(xf(x
where for notational convenience the dependence of kernels on hyper-parameters and parameters is dropped. The posterior variances used as good indicators of how confident one could be about the estimated parameters
of the linear operator
and predictions made based on these parameters. Furthermore, such built-in quantification of uncertainty encoded in the posterior variances is a direct consequence of the Bayesian approach adopted in this work. Although not pursued here, this information is very useful in designing a data acquisition plan, often referred to as active learning [33, 34, 35], that can be used to optimally enhance our knowledge about the parametric linear equation under consideration.
The proposed algorithm provides an entirely agnostic treatment of linear operators, which can be of fundamentally different nature. For example, one can seamlessly learn parametric integro-differential, time-dependent, or even fractional equations. This generality will be demonstrated using three benchmark problems with simulated data. Moreover, the methodology will be applied to a fundamental problem in functional genomics, namely determining the structure and dynamics of genetic networks based on real expression data [36] on early Drosophila melanogaster development.
Figure 1: Integro-differential equation in 1D: (A) Exact left-hand-side function u(x), “noise-free” training data , predictive mean u(x), and two-standard-deviation band around the mean. (B) Exact right-hand-side function f(x), “noise-free” training data
, predictive mean f(x), and two-standard-deviation band around the mean. (C) Exact left-hand-side function u(x), “noisy” training data
, predictive mean u(x), and two-standard-deviation band around the mean. (D) Exact right-hand-side function f(x), “noisy” training data
, predictive mean f(x), and two-standard-deviation band around the mean.
3.1. Integro-differential equation in 1D
Consider the following integro-differential equation,
Note that for (5), the functions
2
) satisfy this equation. In the following, the parameters (
) will be infered from two types of data, namely, noise-free and noisy observations.
Noise-free data
Assume that the noise-free data generated according to
constituting of
spectively. Given these noise-free training data, the algorithm learns the parameters (
) to have values (2.012627, 4.977879). It should be emphasized that the algorithm is able to learn the parameters of the operator using only 7 training data. Moreover, the resulting posterior distributions for u(x) and f(x) are depicted in Figure 1(A, B). The posterior variances could be used as good indicators of how uncertain one should be about the estimated parameters and predictions made based on these parameters.
Noisy data
Consider the case where the noisy data are generated according to
constituting of
interval [0, 1], respectively. Here, the noise
are randomly generated according to the normal distributions
respectively. Given these noisy training data, the algorithm learns the parameters (
) to have values (2.073054, 5.627249). It should be emphasized that for this example the data is deliberately chosen to have a sizable noise. This highlights the ability of the method to handle highly noisy observations without any modifications. The resulting posterior distributions for u(x) and f(x) are depicted in Figure 1(C, D). The posterior variances not only quantify scarcity in observations but also signify noise levels in data.
Figure 2: Heat equation: (A) Exact left-hand-side function u(t, x) and training data Exact right-hand-side function f(t, x) and training data
(C) Absolute point-wise error betrween the predictive mean u(t, x) and the exact function u(t, x). The relative
error for the left-hand-side function is 1
point-wise error betrween the predictive mean f(t, x) and the exact function f(t, x). The relative
error for the right-hand-side function is 4
deviations
, respectively.
3.2. Heat Equation
This example is chosen to highlight the capability of the proposed framework to handle time-dependent problems using only scattered space-time observations. To this end, consider the heat equation
Note that for ) satisfy this equation. Assume that the noise-free data
) are generated according to
) constituting of
spectively. Given these training data, the algorithm learns the parameter
to have value 0.999943. The resulting posterior distributions for u(t, x) and f(t, x) are depicted in Figure 2. A visual inspection of this figure illustrates how closely uncertainty in predictions measured by posterior variances (see Figure 2(E, F)) correlate with prediction errors (see Figure 2(C, D)). Remarkably, the proposed methodology circumvents the need for temporal discretization, and is essentially immune to any restrictions arising due to time-stepping, e.g., the fundamental consistency and stability issues in classical numerical analysis.
3.3. Fractional Equation
Consider the one dimensional fractional equation
where is the fractional order of the operator that is defined in the Riemann-Liouville sense [37]. Fractional operators often arise in modeling anomalous diffusion processes. Their non-local behavior poses serious computational challenges as it involves costly convolution operations for resolving the underlying non-Markovian dynamics [37]. However, the machine leaning approach pursued in this work bypasses the need for numerical discretization, hence, overcomes these computational bottlenecks, and can seamlessly handle all such linear cases without any modifications. The only technicality induced by fractional operators has to do with deriving the kernel
) was obtained by taking the inverse Fourier transform [37] of
Figure 3: Fractional equation in 1D: (A) Exact left-hand-side function u(x), training data , predictive mean u(x), and two-standard-deviation band around the mean. (B) Exact right-hand-side function f(x), training data
, predictive mean f(x), and two-standard-deviation band around the mean.
where ) is the Fourier transform of the kernel
larly, one can obtain
Note that for 2, the functions
and ) satisfy the fractional equation. Assume that the noise-free data
) are generated according to
constituting of
1], respectively. Given these noise-free training data, the algorithm learns the parameter
have value 1.412104. The resulting posterior distributions for u(x) and f(x) are depicted in Figure 3.
3.4. Drosophila melanogaster gap gene dynamics [36, 16]
The gap gene dynamics of protein (see figure 4) can be modeled using a reaction-diffusion partial differential equation
where ) denotes the relative concentration of gap protein a (unitless, ranging from 0 to 255) at space point x (from 35% to 92% of embryo length) and time t (0 min to 68 min after the start of cleavage cycle 13). Here,
Figure 4: Maternal and Gap Gene Expression (see [36]): (A–C) Drosophila embryos at early blastoderm stage (cleavage cycle 13) fluorescently stained for Bcd (A), Cad (B), and Hb (C) protein. (D–H) Drosophila embryos at late blastoderm stage (late cleavage cycle 14A) fluorescently stained for Tll (D), Hb (E), Kr (F), Kni (G), and Gt (H) protein. Anterior is to the left, dorsal is up. Black bars indicate the modeled A-P extent. (I–L) Mean relative gap protein concentration as a function of A-P position (measured in percent embryo length) for Hb (I), Kr (J), Kni (K), and Gt (L). Expression levels are from images and are unitless, ranging from 0 to 255. Images and expression profiles are from the FlyEx database [38]. Embryo IDs: bd3 (A,B), hz30 (C), tb6 (D), kf9 (E), kd17 (F), fq1 (G), nk5 (H). Abbreviations: A-P, anterior-posterior; Bcd, Bicoid; Cad, Caudal; Gt, Giant; Hb, hunchback; Kni, Knirps; Kr, Kr¨uppel; Tll, tailless.
and are decay and diffusion rates of protein a, respectively. Moreover, the right-hand-side is given by
where the term
models the doubling of nuclei and shutdown of transcription during mitosis and
specifies the production rate of protein a. The model combines the processes of transcription and translation into a single production process
Here, is the maximum production rate,
and ranges over the seven genes (see figure 4). The regulatory weights
, encode the effect protein b has on the production rate of protein
0), then gene b is interpreted as being an activator (or a repressor) of gene a.
Table 1: Parameters are assumed to be exogenously given and their values are taken from [36].
This work assumes the maximum production rate , the regulatory weights
, and the bias or offset
to be specified as in table 3.4 and seeks to learn the decay
and diffusion
rates of protein a. In fact, table 3.4 summarizes the values learned by the algorithm for these parameters and figure 5 depicts the corresponding posterior distributions for
). Indeed, figure 5 gives a good indication of how certain one could be about the estimated parameters and the predictions made based on them.
Table 2: Inferred parameter values for the decay and diffusion
rates of protein a.
Figure 5: Predictive expression at seven points in time: (A) Predictive mean expression along with the two-standard-deviations band around the mean for Hb, Kr, Gt, and Kni. The vertical axis represents relative protein concentration corresponding to fluorescence intensity from quantitative gene expression data [38, 36]. (B) Predictive mean along with the two-standard-deviations band around the mean for the right-hand-side function corresponding to Hb, Kr, Gt, and Kni. The horizontal axis in each plot is A-P position, ranging from 35% to 92% of embryo length. No data points are available at time t = 33 min.
In summary, this work introduced a possibly disruptive probabilistic technology for learning general parametric linear equations from noisy data. This generality was demostrated using various bechmark problems with utterly different attributes along with an example application in functional genomics. Furthermore, the current methodology can be applied to inverse problems involving characterization of materials, tomography and electrophysiology, design of effective metamaterials, etc. The methodology can be straightforwardly generalized to address data with multiple levels of fidelity [24, 39] and equations with variable coefficients and complex geometries. Non-Gaussian and input-dependent noise models (e.g., student-t, heteroscedastic, etc.) [3] can also be accommodated. Moreover, systems of linear integro-differential equations can be addressed using multi-output Gaussian process regressions [40, 23, 22]. These scenarios are all feasible because they do not affect the key observation that any linear transformation of a Gaussian process is still a Gaussian process. In its current form, despite its generality regarding linear equations, the proposed framework cannot deal with non-linear equations. However, some specific non-linear operators can be addressed with extensions of the current framework by transforming such equations into systems of linear equations [41, 42] – albeit in high dimensions. In the end, the proposed methodology in this work, being essentially a regression technology, is suitable for resolving such high-dimensional problems.
[1] J. Kaipio, E. Somersalo, Statistical and computational inverse problems, volume 160, Springer Science & Business Media, 2006.
[2] A. M. Stuart, Inverse problems: a bayesian perspective, Acta Numerica 19 (2010) 451–559.
[3] C. K. Williams, C. E. Rasmussen, Gaussian processes for machine learn- ing, the MIT Press 2 (2006) 4.
[4] K. P. Murphy, Machine learning: a probabilistic perspective, MIT press, 2012.
[5] P. Diaconis, Bayesian numerical analysis, Statistical decision theory and related topics IV 1 (1988) 163–175.
[6] V. Vapnik, The nature of statistical learning theory, Springer Science & Business Media, 2013.
[7] B. Sch¨olkopf, A. J. Smola, Learning with kernels: support vector ma- chines, regularization, optimization, and beyond, MIT press, 2002.
[8] M. E. Tipping, Sparse bayesian learning and the relevance vector ma- chine, The journal of machine learning research 1 (2001) 211–244.
[9] A. Tikhonov, Solution of incorrectly formulated problems and the reg- ularization method, in: Soviet Math. Dokl., volume 5, pp. 1035–1038.
[10] A. N. Tikhonov, V. Y. Arsenin, Solutions of Ill-posed problems, W.H. Winston, 1977.
[11] T. Poggio, F. Girosi, Networks for approximation and learning, Pro- ceedings of the IEEE 78 (1990) 1481–1497.
[12] C. Franke, R. Schaback, Solving partial differential equations by collo- cation using radial basis functions, Applied Mathematics and Computation 93 (1998) 73–82.
[13] G. E. Fasshauer, Q. Ye, A kernel-based collocation method for elliptic partial differential equations with random coefficients, in: Monte Carlo and Quasi-Monte Carlo Methods 2012, Springer, 2013, pp. 331–347.
[14] H. Owhadi, Bayesian numerical homogenization, Multiscale Modeling & Simulation 13 (2015) 812–828.
[15] J. Cockayne, C. Oates, T. Sullivan, M. Girolami, Probabilistic meshless methods for partial differential equations and bayesian inverse problems, arXiv preprint arXiv:1605.07811 (2016).
[16] M. A. ´Alvarez, D. Luengo, N. D. Lawrence, Linear latent force models using gaussian processes, IEEE transactions on pattern analysis and machine intelligence 35 (2013) 2693–2705.
[17] M. A. Alvarez, D. Luengo, N. D. Lawrence, Latent force models., in: AISTATS, volume 12, pp. 9–16.
[18] N. Lawrence, Probabilistic non-linear principal component analysis with gaussian process latent variable models, Journal of Machine Learning Research 6 (2005) 1783–1816.
[19] N. D. Lawrence, Gaussian process latent variable models for visualisa- tion of high dimensional data, Advances in neural information processing systems 16 (2004) 329–336.
[20] M. K. Titsias, N. D. Lawrence, Bayesian gaussian process latent variable model., in: AISTATS, pp. 844–851.
[21] D. Higdon, Space and space-time modeling using process convolutions, in: Quantitative methods for current environmental issues, Springer, 2002, pp. 37–56.
[22] P. Boyle, M. Frean, Dependent gaussian processes, in: Advances in neural information processing systems, pp. 217–224.
[23] M. Alvarez, N. D. Lawrence, Sparse convolved gaussian processes for multi-output regression, in: Advances in neural information processing systems, pp. 57–64.
[24] M. C. Kennedy, A. O’Hagan, Predicting the output from a complex computer code when fast approximations are available, Biometrika 87 (2000) 1–13.
[25] L. Le Gratiet, J. Garnier, Recursive co-kriging model for design of com- puter experiments with multiple levels of fidelity, International Journal for Uncertainty Quantification 4 (2014).
[26] N. Aronszajn, Theory of reproducing kernels, Transactions of the Amer- ican mathematical society 68 (1950) 337–404.
[27] S. Saitoh, Theory of reproducing kernels and its applications, volume 189, Longman, 1988.
[28] A. Berlinet, C. Thomas-Agnan, Reproducing kernel Hilbert spaces in probability and statistics, Springer Science & Business Media, 2011.
[29] D. C. Liu, J. Nocedal, On the limited memory bfgs method for large scale optimization, Mathematical programming 45 (1989) 503–528.
[30] C. E. Rasmussen, Z. Ghahramani, Occam’s razor, Advances in neural information processing systems (2001) 294–300.
[31] E. Snelson, Z. Ghahramani, Sparse gaussian processes using pseudo- inputs, in: Advances in neural information processing systems, pp. 1257–1264.
[32] J. Hensman, N. Fusi, N. D. Lawrence, Gaussian processes for big data, arXiv preprint arXiv:1309.6835 (2013).
[33] D. A. Cohn, Z. Ghahramani, M. I. Jordan, Active learning with statis- tical models, Journal of artificial intelligence research (1996).
[34] A. Krause, C. Guestrin, Nonmyopic active learning of gaussian processes: an exploration-exploitation approach, in: Proceedings of the 24th international conference on Machine learning, ACM, pp. 449–456.
[35] D. J. MacKay, Information-based objective functions for active data selection, Neural computation 4 (1992) 590–604.
[36] T. J. Perkins, J. Jaeger, J. Reinitz, L. Glass, Reverse engineering the gap gene network of drosophila melanogaster, PLoS Comput Biol 2 (2006) e51.
[37] I. Podlubny, Fractional differential equations: an introduction to frac- tional derivatives, fractional differential equations, to methods of their solution and some of their applications, volume 198, Academic press, 1998.
[38] E. Poustelnikova, A. Pisarev, M. Blagov, M. Samsonova, J. Reinitz, A database for management of gene expression data in situ, Bioinformatics 20 (2004) 2212–2221.
[39] L. Le Gratiet, Multi-fidelity Gaussian process regression for computer experiments, Ph.D. thesis, Universit´e Paris-Diderot-Paris VII, 2013.
[40] M. A. Osborne, S. J. Roberts, A. Rogers, S. D. Ramchurn, N. R. Jen- nings, Towards real-time information processing of sensor network data using computationally efficient multi-output gaussian processes, in: Proceedings of the 7th international conference on Information processing in sensor networks, IEEE Computer Society, pp. 109–120.
[41] R. Zwanzig, Ensemble method in the theory of irreversibility, The Journal of Chemical Physics 33 (1960) 1338–1341.
[42] A. J. Chorin, O. H. Hald, R. Kupferman, Optimal prediction and the mori–zwanzig representation of irreversible processes, Proceedings of the National Academy of Sciences 97 (2000) 2968–2973.