Linking Gaussian Process regression with data-driven manifold embeddings for nonlinear data fusion

TL;DR

This paper explores advanced Gaussian Process methods for nonlinear data fusion, leveraging manifold embeddings and derivatives of low-fidelity models to improve high-fidelity function approximation with limited data.

Contribution

It introduces mathematical algorithms that incorporate manifold embeddings and derivatives to enhance multifidelity Gaussian Process regression for complex, non-simple high-fidelity functions.

Findings

Using derivatives and shifts of low-fidelity models improves high-fidelity approximation.

Connections established between Gaussian Processes and embedology techniques.

Enhanced data fusion methods require fewer high-fidelity data points.

Abstract

In statistical modeling with Gaussian Process regression, it has been shown that combining (few) high-fidelity data with (many) low-fidelity data can enhance prediction accuracy, compared to prediction based on the few high-fidelity data only. Such information fusion techniques for multifidelity data commonly approach the high-fidelity model $f_h(t)$ as a function of two variables $(t,y)$, and then using $f_l(t)$ as the $y$ data. More generally, the high-fidelity model can be written as a function of several variables $(t,y_1,y_2....)$; the low-fidelity model $f_l$ and, say, some of its derivatives, can then be substituted for these variables. In this paper, we will explore mathematical algorithms for multifidelity information fusion that use such an approach towards improving the representation of the high-fidelity function with only a few training data points. Given that $f_h$ may not be a simple function -- and sometimes not even a function -- of $f_l$, we demonstrate that using additional functions of $t$, such as derivatives or shifts of $f_l$, can drastically improve the approximation of $f_h$ through Gaussian Processes. We also point out a connection with "embedology" techniques from topology and dynamical systems.

Figures

Equations

$$f_h(t)={\sin }^2(8\pi t),f_l(t)=\sin (8\pi t),$$

$$g(x,y)\sim GP\left(0,k((x,y),(x^{\prime },y^{\prime }))\right),f_h(x)=g(x,f_l(x)).$$

$$\left[\begin{array}{c}{f}_l(x)\\ f_h(x)\end{array}\right]\sim GP\left(\left[\begin{array}{c}0\\ 0\end{array}\right],\left[\begin{array}{cc}{k}_1(x,x^{\prime })& \rho k_1(x,x^{\prime })\\ \rho k_1(x,x^{\prime })& {\rho }^2{k}_1(x,x^{\prime })+k_2(x,x^{\prime })\end{array}\right]\right).$$

$$k(x,x^{\prime };\theta )=\exp \left(-\frac{1}{2}\sum ^d{\theta }_i(x_i-x_i^{\prime }{)}^2\right).$$

$$f_h(t)=t^2+{\sin }^2(8\pi t),f_l(t)=\sin (8\pi t)t\in [0,0.25].$$

$$g(t,y)=t^2+y^2.$$

$$\varphi (p)={\left[h(p),h(f(p,-\Delta t)),\dots ,h(f(p,-(n-1)\Delta t))\right]}^T.$$

$$f_h(x)=g(x,f_l(x),f_l^{(1)}(x),\dots ,f_l^{(K)}(x)),$$

$${\Vert \frac{\partial }{\partial x_l}g(x_1,\dots ,x_{K+1})\Vert }_{L^{\infty }}\le c_g\forall l\in \{1,\dots ,K+2\}.$$

$$f_h(x)=f_h(0)+\frac{\partial }{\partial x}{f}_h(0)x+\mathcal{O}\left(\Vert \frac{{\partial }^2}{\partial x^2}{f}_h\Vert \right).$$

$$\frac{\partial }{\partial x}{f}_h(x)$$

$$=$$

$$f_h(x)$$

$$f_h(t)=t^2+{\sin }^2(8\pi t+\pi /10),f_l(t)=\sin (8\pi t).$$

$$f_h(t)=t^2+(\sin (8\pi t)\cos (\pi /10)+\cos (8\pi t)\sin (\pi /10){)}^2.$$

$$f_h(t)\mapsto g(t,f_l(t),f_l(t-\tau ),f_l(t-2\tau )).$$

$$f_h(t)\mapsto g(f_l(t),f_l(t-\tau ),f_l(t-2\tau )).$$

$$f_h(t)=\sin (8\pi t+\pi /10),f_l(t)=\sin (6\sqrt{2}\pi t).$$

$$f_h(t)=\sin (8\pi t)\cos (\pi /10)+\cos (8\pi t)\sin (\pi /10).$$

$$\sin (at-bt)=\cos (bt)\sin (at)-\cos (at)\sin (bt),$$

$$\sin (8\pi t)=\cos (bt)f_l(t)-\sin (bt)f_l^{(1)}(t).$$

$$f_l(t)=0.5(6t-2{)}^2\sin (12t-4)+10(t-0.5)-5,$$

$$f_l(t)=0.5(6t-2{)}^2\sin (12t-4)+10(t-0.5),$$

$$f_h(t)=2f_l(t)-20t+20=g(t,f_l(t)).$$

$$C_m\frac{d{V}_m}{dt}+I_{ion}=I_{ext},$$

$$I_{ion}={\bar{g}}_{Na}{m}^{3}h(V_m-E_{Na})-{\bar{g}}_K{n}^4(V_m-E_K)-{\bar{g}}_L(V_m-E_L),$$

$$\frac{dn}{dt}={\alpha }_n(V)(1-n)-{\beta }_n(V)n,$$

$$\frac{dm}{dt}={\alpha }_m(V)(1-m)-{\beta }_m(V)m,$$

$$\frac{dh}{dt}={\alpha }_h(V)(1-h)-{\beta }_h(V)h.$$

Full text

††preprint: APS/123-QED

I Introduction

Recent advances in both algorithms and hardware are increasingly making machine learning an important component of mathematical modeling for physicochemical, engineering, as well as biological systems (e.g.  Tarca07 ; Holzinger14 ; Libbrecht15 ; Villoutreix17 ). Part of these developments focus on multiresolution and multifidelity data fusion Lanckriet04 ; Willett13 . Fusing information from models constructed at different levels of resolution/fidelity has been shown to enhance prediction accuracy in data-driven scientific computing (e.g. lee17R ; lee17D ).

In addition, if high-fidelity data are costly to obtain (experimentally or computationally) while low-fidelity data are relatively cheap, a combination of a few high-fidelity data and many low-fidelity data can also lead to overall computational efficiency. Recently, Gaussian Process (GP) regression has been widely used to effectively combine multiple fidelity data le14 ; perdikaris15 ; perd16 ; parussini17 ; perdikaris17 . “Classical” information fusion by GP had focused only on linear correlations between high- and low-fidelity data via auto-regressive schemes such as the coKriging approach of Kennedy and O’Hagan kennedy00 . If two (or more) data sets have nonlinear correlations, such linear-type approaches will lose their effectiveness.

When a high-fidelity model $f_{h}(t)$ is nonlinearly correlated with a low-fidelity model $f_{l}(t)$, Nonlinear Auto-Regressive Gaussian Process (NARGP) perdikaris17 achieves highly accurate prediction results by introducing an additional dimension: $f_{h}(t)$ is approximated as a curve on a two-dimensional manifold parametrized by $t$ and a second variable, $y$. The data points on this manifold are of the form $f_{h}(t),t,f_{l}(t)$. More generally, NARGP finds a strong correlation between $d$-dimensional high- and low-fidelity models in a $(d+1)$-dimensional space. In this framework, the low-fidelity model provides the additional “latent” variable of the high-fidelity model.

Deep multifidelity GPs were introduced raissi16 as an improvement, especially in the context of discontinuities in $f_{h}$ with respect to $t$ (and $f_{l}$). This approach focuses on constructing a useful transformation $T(t)$ of the input variables $t$ through a (deep) neural network. Then, the high-fidelity function $f_{h}(t)$ is approximated as a GP of (just) $f_{l}(T(t))$. One must now, of course, perform optimization for the additional network hyperparameters.

In this paper we discuss a connection of NARGP with data-driven embeddings in topology/dynamical systems, and extend it (in the spirit of such data driven embeddings) in an attempt to improve the numerical approximation of $f_{h}$ in the “sparse $f_{h}$, rich $f_{l}$” data setting. In what follows we will (rather arbitrarily) ascribe the characterization “high fidelity” or “low fidelity” to different functions used in our illustrations; this characterization is solely based on the number of available data points for each.

For our first example, the “high-fidelity” function $f_{h}$, for which we only have a few data points, is a function of $t$; but it actually also happens that we can describe it as a function of $f_{l}$,

[TABLE]

see FIG.1.

When we choose $t$ as the coordinate parametrizing $f_{h}(t)$ (the green curve in FIG.1(a)), the GP regression fails to represent the high-frequency sine function with just a few training data points as shown in FIG.1(b). However, as FIG.1(c) shows, if we choose $f_{l}$ as the coordinate of $f_{h}(f_{l})$ (colored by red in FIG.1(a)), the GP regression can represent the simple quadratic function quite effectively. If we still need to know the parametrization of $f_{h}$ by $t$, we can obtain it through the “data rich” $f_{l}$: $f_{h}(t)\equiv f_{h}(f_{l}(t))$.

If $f_{h}$ is not a function of $f_{l}$, however (as has been observed in the literature perdikaris17 and as can be seen in Figure 2(a)) more variables are necessary to create a useful parametrization of $f_{h}$ - more precisely, a domain over which $f_{h}$ is a function.

In the NARGP framework, the variable used in addition to $f_{l}$ is $t$ itself. In this paper, we will also advocate the use of delays or derivatives of $f_{l}$ as additional variables. This approach can also help remove the explicit dependence of $f_{h}$ on $t$, since embedding theories in dynamical systems Whitney36 ; Nash66 ; Takens81 ; kostelich90 ; sauer91 ; kennel92 guarantee that we can choose any generic observation of $t$, or derivatives and/or delays of this observation, as a replacement for $t$, see section II.2 for more details.

In section III.4, we apply the proposed framework to the Hodgkin-Huxley model, describing the behavior of action potentials in a neuron. Here, the high- and the low-fidelity functions are action potentials at two different values of the external current. This is a case where $f_{h}$ is a complicated function of $t$ and does not only depend on $f_{l}$; yet as we will see, delays of $f_{l}$ will help us construct an accurate approximation of $f_{h}$.

The paper is organized as follows: In section II, we review the NARGP framework and concepts of “embedology”. Also, we illustrate why and when this framework is successful. In section III, we demonstrate the effectiveness of our proposed framework via some pedagogical examples and the biologically motivated application to the Hodgkin-Huxley model. In section IV, we summarize our results and discuss open issues for further development of general multifidelity information fusion in modeling and simulation practice across scientific domains.

II Methods

II.1 Nonlinear information fusion algorithms

“Classical” multifidelity data fusion algorithms require a linear (or almost linear) correlation between different fidelity models. Under this constraint, we can merge two or more data sets by using scaling and shifting parameters such as the Kennedy and O’Hagan coKriging approach kennedy00 . However, more generally, the data sets are nonlinearly correlated, which typically degrades the quality of results of linear information fusion algorithms. In order to resolve correlation nonlinearity between data sets, the use of a space-dependent scaling factor $\rho(x)$ le13 or, alternatively, deep multifidelity Gaussian Processes raissi16 have been introduced; clearly, the improvement they bring requires additional hyperparameter optimization.

When the high-fidelity model $f_{h}$ is nonlinearly correlated with the low-fidelity model $f_{l}$, but can be written as a simple function of $x$ and $f_{l}$, the Nonlinear Auto-Regressive GP (NARGP) perdikaris17 is an appropriate choice. In this framework, a one dimensional high-fidelity function $f_{h}$ is assumed to be a “simple” function $g$ of two variables $(x,y)$, i.e. it is a curve that lies in the two-dimensional manifold described by $g$. Then, GP regression in the two dimensional space is performed, where the data for $y$ are the $f_{l}(x)$ data,

[TABLE]

In  perdikaris17 the gain in accuracy of NARGP, compared to an auto-regressive scheme with a constant scaling factor, as well as a scaling factor that was modeled as a (space-dependent) Gaussian process, was documented.

Algorithmically, classical autoregressive GPs employ an explicit method such as a scaling constant ($\rho$) between two covariance kernels, $k_{1}$ and $k_{2}$ as

[TABLE]

The NARGP framework, on the other hand, employs an implicit approach by the automatic relevance determination (ARD) weight rasmussen06 in the extended space parametrized by $x$ and $y$: a different scaling hyperparameter for each of the two dimensions in the kernel. In many applications, the radial basis function (see the equation (4)) has been used as the covariance kernel (where ARD implies a different scaling hyperparameter $\theta_{i}$ for each dimension):

[TABLE]

FIG.2 showcases an example where $f_{h}$ cannot be written as a function of $f_{l}$:

[TABLE]

Following the NARGP framework, we choose the low-fidelity data as an additional variable $y=f_{l}(t)$; we then approximate the two-dimensional function

[TABLE]

Approximating $g$ only requires a few training data point pairs for the GP regression. Then, $f_{h}$ can be written as $f_{h}(t)=g(t,f_{l}(t))$, see FIG.2(b). FIG.2(c) demonstrates that we can, alternatively, use delays of $f_{l}$ instead of $t$ as an additional variable. A rationalization of this follows in the next section.

II.2 Data-driven Higher-dimensional Embeddings

The theorem of Whitney Whitney36 states that any sufficiently smooth manifold of dimension $d\in\mathbb{N}$ can be embedded in Euclidean space $\mathbb{R}^{n}$, with the tight bound $n\geq 2d+1$. Nash Nash66 showed that this embedding can even be isometric if the manifold is compact and Riemannian, even though the bound on $n$ is higher. Many results on the reconstruction of invariant sets in the state spaces of dynamical systems are based on these two theorems Takens81 ; kostelich90 ; sauer91 ; kennel92 . Here, the $n$ embedding dimensions are usually formed by $n$ scalar observations of the system state variables. Instead of $n$ different observations, recording $n$ time delays of a single scalar observable is also possible, as originally formulated by Takens Takens81 (see also Packard80 .)

Given a smooth, $d$-dimensional manifold $M$, a vector field $X$ on $M$ with associated flow map $f:M\times\mathbb{R}\to M$, as well as an observable $h:M\to\mathbb{R}$, it is possible to construct an embedding $\phi:M\to\mathbb{R}^{n}$ through

[TABLE]

Sauer et al. sauer91 further specified the conditions on the observable, with respect to the underlying vector field of the dynamical system, that have to be satisfied such that the state space can be embedded successfully. They also extended the results on embeddings of invariant sets with fractal dimension. The observable and the vector field together build a tuple that has to be generic (Takens Takens81 ) or prevalent (Sauer et al. sauer91 ). The two papers show that “almost all” tuples are admissible, where the notions of genericity and prevalence are defining that probabilistic statement in the respective function spaces. These results are crucial for the numerical exploitation of the embedding theorems, since they show that “almost all” observables of a dynamical system can be chosen for a delay embedding.

In many applications, the intrinsic dimension $d$ of $M$ is unknown, and $n$ is preemptively chosen large (larger than necessary). Manifold learning techniques are then often capable of reducing the embedding dimension, bringing it closer to the minimal embedding dimension necessary for the given manifold.

We also note that in our framework, if we can obtain low-fidelity data from a (dynamic) process instead of just single measurements, the same embedding and manifold learning techniques can be used even if the “independent variable” $t$ is not known Dietrich18 .

Let us consider the equation (5) again. NARGP performs GP regression with two observations $\{t,f_{l}(t)\}$. Using embedding theory, we can rationalize (a) why the two observations were necessary; and (b) why performing GP regression with an additional delay of the scalar observable, $\{f_{l}(t),f_{l}(t-\tau)\}$, is equally appropriate for a relatively small time horizon, see FIG.2(c). Notice that delay coordinates $f_{l}(t)$ and $f_{l}(t-\tau)$ lie on an ellipse with period of 0.25. Hence, if the data are collected over times longer than $0.25$, using only delays will fail to represent the high-fidelity function due to multivaluedness (see FIG.3(e) in section III.1); $t$ itself used as an observable will resolve this.

II.3 Extending the formulation through delays

Now we provide a mathematical formulation of the approach. We assume that $f_{h}:\mathbb{R}\to\mathbb{R}$ and $f_{l}:\mathbb{R}\to\mathbb{R}$ are $C^{K}$-smooth functions with $K\geq 2$, and we want to construct an interpolant of $f_{h}$ with a numerical scheme (here, a Gaussian Process). We assume that

1. (1)

only a small number of data points $\{x,f_{h}(x)\}$ as well as the function $f_{l}$ and its $K$ derivatives are available, 2. (2)

$f_{h}$ can be written in the form

[TABLE]

where $f_{l}^{(i)}(x)$ denotes the $i$-th derivative, $i\in\{1,\dots,K\}$, and 3. (3)

$g:\mathbb{R}^{K+2}\to\mathbb{R}$ is a $C^{K}$ function with derivatives bounded in the $L^{\infty}$ norm by a small constant $c_{g}>0$,

[TABLE]

The Taylor series of $f_{h}$ reveals why assumptions (2) and (3) are required: Assume we want to evaluate $f_{h}$ in a small neighborhood of a point $x_{0}\in\mathbb{R}$ (here, w.l.o.g. we set $x_{0}=0$). Then, assuming $f_{h}$ is smooth enough, we can write

[TABLE]

Since we assume the form (8) for $f_{h}$, we can write

[TABLE]

From this, we can see that $\frac{\partial}{\partial x}f_{h}(x)$ can be large (because the derivatives $f_{l}^{(i)}(x)$ can be large), but if we know all $f_{l}^{(i)}(x)$, we only have to estimate $g$ from data $f_{h}$ and $f_{l}^{(i)}$. Crucially, we do not approximate the function $f_{h}(x)$ and its derivatives. The derivatives of $g$ are bounded by $c_{g}$ through assumption (8), and $g$ is a $C^{K}$ function, so only a few data points are necessary for a good fit. If we have access to function values of $f_{l}$ over “delays” (at discrete shifts, say in the form of a finite difference stencil) in space, rather than its derivatives, we can use the Newton series approximation of $f_{h}$ instead of equation (10) for an analogous argumentation: For functions $f_{h}$ that are analytic around the expansion point (here, $x_{0}=0$), the function $f_{h}$ can be evaluated at $x$ close to it by

[TABLE]

where $\Delta_{m}f_{h}$ is the $m$-th finite difference approximation of $f_{h}$, with a small step size $\Delta x$. By equation (8), these differences can be expressed through $g$ and delays of $f_{l}$ (instead of delays of $f_{h}$), analogously to equations (11–13). Using delays in space compared to derivatives has numerical advantages, especially in cases where $f_{h}$ or $f_{l}$ are not differentiable (or even have discontinuities). It also enables us to estimate the derivatives of $f_{l}$ implicitly, in case only the function $f_{l}$ is available (and not its derivatives).

II.4 Outline of the numerical approach

In order to employ the delay coordinates of the low-fidelity function, it is required to know shifts of it. A necessary condition of the proposed framework is that the low-fidelity function is given explicitly or can be well-learned by given data such that low-fidelity function values can be accurately approximated (interpolated) at arbitrary points. If the state variable $t$ is not available, the low fidelity model should be a generic observation of $t$ to be useful in employing Takens’ embedding theorem sauer91 ; Takens81 . Under these conditions, we now present a summary of the workflow.

If the low-fidelity model is given in the form of (rich) data, we train a GP regression model for it from these data $\{\mathbb{t_{l}},\mathbb{y_{l}(t_{l})}\}$ via minimizing a negative log marginal likelihood estimation. This data driven process can be circumvented if the low-fidelity model is explicitly given, as in the above examples. After that, we compute predictive posterior means of the low-fidelity model at the points $t_{h}$ where the high-fidelity data is available. We also compute a number of shifts of the low-fidelity function at the points $t_{h}-k\tau$ and at the test points $t^{*}$. Next, we train another GP regression model for high-fidelity data in the higher-dimensional space, $\{\{\mathbb{t_{h}},\mathbb{f_{l}(t_{h})},\mathbb{f_{l}(t_{h}-\tau)},\dots,\mathbb{f_{l}(t_{h}-n\tau)}\},\mathbb{f_{h}(t_{h})}\}$. Finally, we compute the predictive posterior mean and variance at the test points ($t^{*}$) in the higher-dimensional space. Each new delay burdens the optimization by a single additional hyperparameter. For more details, refer to rasmussen06 ; perdikaris17 . In this paper, all GP computations are performed by the open source python package GPy gpy14 .

III Results

We introduce three pedagogical examples to demonstrate our approach. First, we explore the case where $f_{h}$ is a phase shifted version of $f_{l}$ (section III.1). Then, we show that oscillations with different periods (leading to different recurrences) present a more challenging scenario, which however can still be resolved by using shifts of $f_{l}$. The third example involves discontinuities in $f_{h}$ and $f_{l}$. After these three examples, in section III.4, we demonstrate the approach in the context of the Hodgkin-Huxley model. We investigate the effectiveness of the proposed framework by comparing it to three established frameworks: (1) single Gaussian Process regression with high-fidelity data only (GP or Kriging), (2) Auto-regressive method with a constant scaling parameter $\rho$ (AR1 or coKriging), and (3) nonlinear auto-regressive GP (NARGP) in the same computational environment.

III.1 Phase Shifted Oscillations

Using models at different levels of resolution (e.g. for biological oscillators) will often give oscillations that have very comparable periods but are phase-shifted. Let us start by considering two functions with different phases on $t\in[0,1]$,

[TABLE]

The high-fidelity function can be rewritten by a trigonometric addition formula:

[TABLE]

Now we can explicitly see how the high-fidelity function can be thought of as a combination of three variables: $t^{2}$, the low-fidelity function $f_{l}(t)=\sin(8\pi t)$ and its first derivative $f_{l}^{(1)}(t)=\cos(8\pi t)$.

Using delays of $f_{l}$ with a small stepsize $\tau$ contains enough information to numerically estimate its derivatives, hence we can also write $f_{h}$ as

[TABLE]

The GP regression model for $g$ is trained on this 4-dimensional data. In addition, we perform GP regression in a three dimensional space constructed from only three delays:

[TABLE]

As shown in FIG.3(b), the single GP regression model provides inaccurate predictive posterior means due to lack of high-fidelity data. While the linear auto-regressive model (AR1) also fails to predict the high fidelity values, the NARGP (with 10 high-fidelity data points and 100 low-fidelity data points) catches the trend of the high fidelity data, yet still yields inaccurate results: NARGP is informed only by $t$ and $f_{l}(t)$, and not also by $f_{l}^{(1)}(t)$. Similarly, the GP regression with only delays (no information about $t$) in FIG.3(e) fails to represent the high-fidelity function for these long observation windows. Beyond $0.25$, $t$ cannot be recovered from the shifts of $f_{l}$ because $f_{l}$ is only a generic observer of $t\in[0,0.25]$.

As shown in FIG.3(f), the GP using $t$ and three delays of $f_{l}$ provides an excellent prediction with only 10 high-fidelity data points (and 100 low-fidelity data points). This means that, in the 4-dimensional space, $g$ (see equation 18) has small derivatives, which then helps to employ GP regression successfully.

Next, we investigate the sensitivity and scalability of the proposed framework on the number of high-fidelity data points (training data points). We train all GP regression models with 10, 15, 20, and 25 randomly-chosen high-fidelity data points and 100 uniformly distributed low-fidelity data points. The error is obtained by averaging 10 trials of random data selections. A log $L_{2}$ error with respect to the number of high-fidelity data points is presented in FIG.4(a).

Two established approaches (AR1 and NARGP) and the GP with only delays have no significant accuracy enhancements as the number of training points increases. The reason for the consistently large errors is the lack of additional information provided by the derivatives. The GP in the higher-dimensional space that includes $t$, on the other hand, shows a strong correlation between accuracy and the number of training points – more high fidelity points visibly improve the approximation.

III.2 Different Periodicity

In this example, the high- and the low-fidelity model oscillations are not just phase shifted, but they also are characterized by different periods. In applications, this could arise if we tried to match observations of oscillations of the same model at two different parameter values. Different (possibly irrationally related) oscillation periods dramatically complicate the correlation across the two data sets.

We consider two different period and phase shifted data,

[TABLE]

The high-fidelity function can be rewritten by a trigonometric addition formula:

[TABLE]

In addition, the first term $\sin(8\pi t)$ can be rewritten again by a trigonometric subtraction formula:

[TABLE]

where $a=6\sqrt{2}\pi$ and $b=6\sqrt{2}\pi-8\pi$. Then,

[TABLE]

The second term $\cos(8\pi t)$ can be rewritten in the same way. This shows that the high-fidelity function can be written in terms of $\sin(bt)$, $\cos(bt)$, $f_{l}(t)$, and $f_{l}^{(1)}(t)$. Since $\sin(bt)$ and $\cos(bt)$ have lower frequency compared to the original frequency $8$, the bound $c_{g}$ for the derivatives of $g$ (see section II.3) is smaller. It is then reasonable that we can approximate the high-fidelity function in the higher-dimensional space with a few training data points.

We perform the GP in the two different extended spaces: (1) 3 additional delays, totaling 4-dimensional space (GP+E), and (2) 5 additional delays, totaling 6-dimensional space (GP+E(2)), and compare them to a single GP, AR1, and NARGP. Examples of regression results with 15 high-fidelity data points and 200 uniformly distributed low-fidelity data are shown in FIG.5. The GP in the 4-dimensional space provides better regression result than other established methods, and the GP in the 6-dimensional space presents the best results.

Moreover, as shown in FIG.5(b), the phase discrepancy between the high- and low-fidelity functions increases as time increases, resulting in larger error for larger values of $t$, see FIG.5(b)–(d). However, the GPs in the higher-dimensional spaces provide accurate prediction results over this time observation window.

The sensitivity of the number of high fidelity data is shown in FIG.4(b). The GPs in the 4- and 6-dimensional space show significant computational accuracy gain compared to all other methods. These results demonstrate the capability of the proposed framework for period-shifted and phase-shifted data fusion.

III.3 Models with Discontinuities

In general, a smooth stationary covariance kernel cannot capture discontinuous model data. In order to resolve this problem, a non-stationary kernel has been introduced schmidt03 ; paciorek04 , with space-dependent hyperparameters. Moreover, nested GPs were also used successfully to alleviate this problem hensman14 . Both approaches introduce, of course, additional hyperparameters to optimize.

In this example, we introduce a discontinuous function $f_{l}$ on $t\in[0,0.5)$,

[TABLE]

and on $t\in[0.5,1]$ as

[TABLE]

and the high-fidelity function $f_{h}$

[TABLE]

In the scenario we describe here, the high-fidelity function $f_{h}$ is discontinuous, but can be written in terms of a linear function of $g$ in two variables.

Examples of regression results with 10 high-fidelity data points and 200 uniformly distributed low-fidelity data points are shown in FIG.6. Since $g$ is a linear function of $t$ and $f_{l}(t)$, the NARGP, as well as our GPs in the higher-dimensional spaces, provide highly accurate prediction results with just a few high-fidelity data.

In the sensitivity analysis of the number of high fidelity data points, see FIG.4(c), there is no significant accuracy gain after 10 high-fidelity data points. That is because 10 training data points are enough to represent a linear function accurately. It is worth noting that, here, the NARGP provides better prediction results with 20 training data points compared to the GPs in the higher-dimensional space, possibly due to overfitting.

III.4 The Hodgkin-Huxley Model

Based on the results of the pedagogical examples, we apply the proposed framework to a famous model of an intracellular process, a version of the Hodgkin-Huxley equations Hodgkin52 . In 1952, Hodgkin and Huxley introduced a mathematical model which can describe the initiation and propagation of action potentials in a neuron. Specifically, they invented electrical equivalent circuits to mimic the ion channel, where ions traffic through the cell membrane. The model for intracellular action potentials ($V_{m}$) can be written as a simple ODE:

[TABLE]

where $C_{m}$ is the membrane capacitance and $I_{ion}$ and $I_{ext}$ represent the total ionic current and the external current, respectively.

The total ionic current $I_{ion}=I_{Na}+I_{K}+I_{L}$ is the sum of the three individual currents as a sodium current ($I_{Na}$), a potassium current ($I_{K}$), and a leakage current ($I_{L}$). In order to calculate the three individual currents in time, the Hodgkin-Huxley model introduced gates which regulate the flow of ions through the channel. Specifically, the three ionic currents are affected by the three different gates $n$, $m$, and $h$. Based on these gates, the total ionic currents can be calculated by

[TABLE]

where $\bar{g}_{*}$ represents a normalized constant and $E_{*}$ represents the equilibrium potential for a sodium, $(*\equiv Na)$, a potassium, $(*\equiv K)$, and a leakage, $(*\equiv L)$, current. The three gates $n$, $m$, and $h$ can be modeled by the following ODEs:

[TABLE]

[TABLE]

[TABLE]

In this paper, we set the model parameter values to $\bar{g}_{Na}=1.2$, $\bar{g}_{K}=0.36$, $\bar{g}_{L}=0.003$, $E_{Na}=55.17$, $E_{K}=-72.14$, and $E_{L}=-49.42$ Siciliano12 . We assume that different fidelity data come from simulations at different external currents $I_{ext}$. We set $I_{ext}=1.0$ for the high fidelity model and $I_{ext}=1.05$ for the low fidelity model, resulting in different oscillation period (and a phase shift when we start at the same initial conditions). The action potentials $V_{m}$ of the two different fidelity models are shown in FIG.7(a,b).

Examples of regression results for 5 different methods with 20 high-fidelity data points and 300 uniformly distributed low-fidelity data points are shown in FIG.7(b–f). Since the two data sets are phase-shifted, the single GP, AR1, and NARGP fail to represent a proper function of the high fidelity model effectively. However, GPs in the higher-dimensional spaces provide reasonable prediction results. The GP in the 6-dimensional space (GP+E(2)) shows significant improvement in the form of large uncertainty reduction as well as high prediction accuracy.

The sensitivity to the number of high fidelity data is shown in FIG.4(d). The proposed framework shows computational accuracy gains compared to all other methods, as well as marked improvement when new high fidelity points are added.

IV Conclusion

In this paper we explored mathematical algorithms for multifidelity information fusion and its links with “embedology”, motivated by the NARGP approach. These modifications/extensions of kriging show promise in improving the representation of data-poor “high-fidelity” datasets exploiting data-rich “low-fidelity” datasets. Given that $f_{h}$ may not be a simple function -and sometimes not even a function- of $f_{l}$, we demonstrated that using additional functions of $t$, such as derivatives or shifts of $f_{l}$, can drastically improve the approximation of $f_{h}$ through Gaussian Processes.

The limitations of the proposed framework arise in the form of the curse of dimensionality and of overfitting. As the number of hyperpameters in the GP framework grows in an increasingly higher dimensional input space, the optimization cost grows (and there is always the possibility of converging to local, unsatisfactory minima). Adaptively testing for the “best” number of delays is possible, and will be pursued in future work. The natural option of using multiple low-fidelity models (instead of delays of just one of them) is also being explored. Techniques that systematically find all the local hyperparameter minima (in the spirit of the reduced gradient method quapp03 ) may also be useful in this effort. Another promising research direction involves the construction of data-informed kernels (e.g. through “Neural-net-induced Gaussian process”, NNGP pang18 ) for more realistic and unbiased predictions. Alternatively, it is interesting to consider transformations of the input space using manifold learning techniques and the so-called Mahalanobis distance singer08 ; singer09 , which has been demonstrated to successfully match different (yet conjugate) models dsilva16 ; kemeth18 .

What we believe is a most promising direction for the use of these techniques is the reconciliation of different granularity multiscale models - having, say, an atomistic “high-fidelity” simulation enhanced by a continuum “low-fidelity” approximate closure. Thus, “heterogeneous data” fusion becomes a version of multifidelity data fusion lee17D . In this paper the fusion tools simply “filled in the gaps” in a single manifestation of the high fidelity data. In a time-dependent context, “full space, full time” low fidelity simulations can help complete and thus accelerate “small space, small time” high fidelity simulations - in a form reminiscent of the patch-dynamics approach in equation-free computation kevrekidis03 , see also lee17D ; lee17R . Using a qualitatively correct (even though quantitatively inaccurate) low-fidelity model –as opposed to just the local Taylor series that play the role of low-fidelity modeling in patch dynamics– may very much improve the computational savings of such multiscale computation schemes.

V Acknowledgments

S. Lee, F. Dietrich and I. G. Kevrekidis gratefully acknowledge the partially support of NSF, NIH, and DARPA. G. E. Karniadakis gratefully acknowledge the financial support of DARPA (grant no. N66001-534 15-2-4055) also the MIT ESRDC grant.