claimClaim remarkRemark hypothesisHypothesis

    Kriging with Partial Differential EquationsJialei Chen, Zhehui Chen, Chuck Zhang and C. F. Jeff Wu

    # APIK: Active Physics-Informed Kriging Model with Partial Differential Equations

    Jialei Chen
    The H. Milton Stewart School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta,
    GA 30332 (, , , ).
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

    Zhehui Chen11footnotemark: 1

    Chuck Zhang11footnotemark: 1
    Georgia Tech Manufacturing Institute.

    C. F. Jeff Wu11footnotemark: 1

    ###### Abstract

    Kriging (or Gaussian process regression) is a popular machine learning method for its flexibility and closed-form prediction expressions.
    However, one of the key challenges in applying kriging to engineering systems is that the available measurement data is scarce due to the measurement limitations and high sensing costs.
    On the other hand, physical knowledge of the engineering system is often available and represented in the form of partial differential equations (PDEs).
    We present in this work a PDE Informed Kriging model (PIK), which introduces PDE information via a set of PDE points and conducts posterior prediction similar to the standard kriging method.
    The proposed PIK model can incorporate physical knowledge from both linear and nonlinear PDEs.
    To further improve learning performance, we propose an Active PIK framework (APIK) that designs PDE points to leverage the PDE information based on the PIK model and measurement data.
    The selected PDE points not only explore the whole input space but also exploit the locations where the PDE information is critical in reducing predictive uncertainty. Finally, an expectation-maximization algorithm is developed for parameter estimation. We demonstrate the effectiveness of APIK in two synthetic examples, a shock wave case study, and a laser heating case study.

    ###### keywords:

     Derivative process; Expectation-maximization algorithm; Gaussian process regression; Pseudo points.

    60G15, 62F15, 62G08, 62K20, 62P25, 62P30

    ##
    1 Introduction

    Recent advances in machine learning have attracted increasing attention from the engineering community.
    Among these learning methods, kriging (or Gaussian process regression) [46, 37] is particularly popular for its flexible structure and the fact that both the prediction and its uncertainty quantification enjoy simple and closed-form expressions.
    Therefore, kriging has been adopted in a variety of applications in different domains, including geostatistics [29, 10], financial engineering [40, 20], and computer experiment [36, 7].
    Despite the popularity, one of the key challenges in real-world engineering applications is that the available amount of data is often extremely limited due to the constraints of measurement devices and human labor.
    For example, in a laser heating process in wafer manufacturing [3, 6], the goal is to understand the temperature profile of a rotating wafer heated by a laser beam (more details in Section 5.4).
    The temperature measurements rely on a handheld temperature gun.
    Only a handful of readings during the heating process is available since (i) the measurement procedure is labor-intensive, and (ii) the response time of the temperature gun is around one second. To address this, we present in this work a new kriging modeling framework, which can improve predictive performance with limited data.

    One way to improve learning performance (we adopt) is to incorporate physical knowledge, which can also be viewed as enlarging the existing data set [45, 16].
    Recently, researchers have also proposed several physics-informed kriging methods to incorporate physical knowledge,
    such as bound constraints [8], invariance equations [4], monotonicity information [34], and output’s correlation structure [28]. However, the physical knowledge considered in those works only exists on a case-by-case basis.
    In this work, we develop a kriging framework with partial differential equation (PDE) shown in Fig. 1, which is a rigorous representation of the first principle and exists for almost all engineering systems.
    Take the previous wafer manufacturing application as an example. Besides the limited data, we know that the temperature profile of the wafer is governed by the Fourier equation [12], a PDE describing macroscopic transport properties of energy.
    Our goal is to incorporate the Fourier equation with measurement data, thereby better predicting the temperature profile.

    Figure 1:  An illustration of the proposed PIK and APIK methods. The PIK method can incorporate the PDE information in the kriging framework, and the APIK model can actively select PDE points to leverage the PDE information according to the measurement data.

    There are several related works in the kriging literature. Wang and Berger [43] developed a kriging models to incorporate derivative constraints, i.e., the simplest PDE information, with follow-up works in [42, 17, 19].
    Wheeler et al. [44] proposed a hierarchical Gaussian process model that favors curves consistent with the considered linear PDEs.
    After that, Jidling et al. [21] developed a kriging model with modification of the correlation function to account for linear PDEs.
    Raissi et al. [32] proposed a kriging-based method to evaluate the unknown coefficients in linear PDEs, with further development in [15].
    Besides, Raissi et al. [33] presented a kriging model for the evolution of dynamic systems with a correlation structure inherited from PDEs.
    There are also works on incorporating boundary condition constraints [39, 11].
    However, to the best of our knowledge, there is no systematic modeling framework for incorporating nonlinear PDEs into kriging methods.

    In this work, we propose a new kriging model to fuse general PDEs and measurement data. We develop first a PDE Informed Kriging model (PIK) to utilize a set of pseudo points, called PDE points, to incorporate physical knowledge from linear PDEs and nonlinear PDEs. Specifically, for linear PDEs, we extend the learning method of incorporating gradient information in [43]. For nonlinear PDEs, we introduce a set of latent variables, which converts the nonlinear PDEs to linear PDEs by the conditional process.
    To further improve the learning performance, we propose an Active PIK framework (APIK) that leverages the PDE information by actively designing the PDE points based on the PIK model and measurement points. The selected PDE points not only explore the whole input space but also exploit the locations where the PDE information is the most important in reducing predictive uncertainty. We then develop an expectation-maximization method for APIK, which efficiently estimates model parameters and latent variables, as well as actively selects the locations for PDE points.
    Further discussions on PDE data size, computational simplification, and the connection to numerical PDE solvers are included.
    Finally, we demonstrate the improvements of the proposed APIK method in both synthetic examples and real case studies.

    The remaining part of the article is organized as follows. In Section 2, we present the PIK model with linear and nonlinear PDEs. In Section 3, we develop the APIK model based on PIK. We discuss a parameter estimation algorithm in Section 4. PIK and APIK are then applied to four applications in Section 5.
    Finally, Section 6 concludes this work.

    ##
    2 PIK model

    In this section, we first review the kriging method. We then present the PIK model with linear PDEs. Finally, we extend the PIK model to nonlinear PDEs.

    ###
    2.1 Kriging model

    Let $\mathbf{x}_{i}\in\mathcal{X}\subset\mathbb{R}^{d}$ be a vector of $d$ input variables in a bounded space $\mathcal{X}$, and let $y_{i}\in\mathbb{R}$ be the corresponding measurement of the physical quantity. We consider the following model for the physical quantity of interest:

    (1)

    $\displaystyle y_{i}=y(\mathbf{x}_{i})+\epsilon_{i},\quad i=1,2,\cdots,n,$

    where $y(\mathbf{x}_{i})$ is the underlying output physical quantity at input $\mathbf{x}_{i}$, $\epsilon_{i}$ is the corresponding measurement error, and $n$ is the size of measurement data. The experimental noise $\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2}_{\epsilon})$ is assumed to be i.i.d. normally distributed and independent to $y(\cdot)$. Kriging model further supposes that the input-output relationship $y(\cdot):\mathcal{X}\mapsto\mathbb{R}$ follows a Gaussian process

    (2)

    $\displaystyle y(\mathbf{x})\sim\mathcal{GP}(\mu(\mathbf{x}),\sigma^{2}R_{\boldsymbol{\theta}}(\cdot,\cdot)),\quad\text{with}\quad\mu(\mathbf{x})=\mathbf{p}^{\top}(\mathbf{x})\boldsymbol{\beta}.$

    Here, $\mathbf{p}(\mathbf{x})=[p_{1}(\mathbf{x}),\cdots,p_{q}(\mathbf{x})]^{\top}$ consists of $q$ basis functions for the mean function $\mu(\mathbf{x})$, $\boldsymbol{\beta}\in\mathbb{R}^{q}$ denotes the corresponding coefficients, $\sigma^{2}$ is the process variance, and $R_{\boldsymbol{\theta}}(\cdot,\cdot)$ is the correlation function with parameter $\boldsymbol{\theta}$.
    Denote the response vector of measurement data $\mathbf{y}_{1:n}=[y_{1},\cdots,y_{n}]^{\top}$. We have the following proposition for the posterior prediction.

    ######
    Proposition 2.1 (Page 17 in [46]).

    Suppose that the Gaussian process model (1)-(2) holds with known model parameters $\{\boldsymbol{\beta},\sigma^{2},\boldsymbol{\theta}\}$.
    Then the posterior prediction $y(\mathbf{x}_{\rm new})|\mathbf{y}_{1:n}$ at a new input location $\mathbf{x}_{\rm new}$ is normally distributed with mean and variance

    (3)

    $\displaystyle\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\mathbf{y}_{1:n}\right]$
    $\displaystyle=\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\mathbf{r}^{\top}(\mathbf{x}_{\rm new})\mathbf{R}_{\boldsymbol{\theta}}^{-1}(\mathbf{y}_{1:n}-\mathbf{P}\boldsymbol{\beta}),$

    (4)

    $\displaystyle\text{\rm Var}\left[y(\mathbf{x}_{\rm new})|\mathbf{y}_{1:n}\right]$
    $\displaystyle=\sigma^{2}(1-\mathbf{r}^{\top}(\mathbf{x}_{\rm new})\mathbf{R}_{\boldsymbol{\theta}}^{-1}\mathbf{r}(\mathbf{x}_{\rm new})).$

    where, $\mathbf{P}=[\mathbf{p}(\mathbf{x}_{1}),\cdots,\mathbf{p}(\mathbf{x}_{n})]^{\top}$, $\mathbf{r}(\mathbf{x}_{\rm new})=[R_{\boldsymbol{\theta}}(\mathbf{x}_{\rm new},\mathbf{x}_{1}),\cdots,R_{\boldsymbol{\theta}}(\mathbf{x}_{\rm new},\mathbf{x}_{n})]^{\top}$, and $\mathbf{R}_{\boldsymbol{\theta}}={[R_{\boldsymbol{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})]_{i=1}^{n}}_{j=1}^{n}+\sigma^{2}_{\epsilon}/\sigma^{2}\mathbf{I}_{n}$ with $\mathbf{I}_{n}$ the $n\times n$ identity matrix.

    Here, $\mathbf{P}$ is the $n\times q$ model matrix at $n$ measurement points, $\mathbf{r}(\mathbf{x}_{\rm new})$ is the correlation vector between the $n$ points and the predictive location, and $\mathbf{R}_{\boldsymbol{\theta}}$ is the corresponding $n\times n$ correlation matrix. The conditional mean in (3) can be used for prediction, and the conditional variance in (4) can be used for constructing point-wise predictive intervals.

    ###
    2.2 PIK with linear PDEs

    As discussed, PDE information on the physical quantity $y(\mathbf{x})$ is available. Consider for now the PDEs are linear; Nonlinear PDEs will be discussed later in Section 2.3. We have the following definition for linear PDEs.

    ######
    Definition 2.2 (Linear PDE).

    PDE $\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})=b(\mathbf{x})$ with a non-homogeneous term $b(\mathbf{x})$ is called a linear PDE if $\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})$ is a linear differential operator:

    (5)

    $\displaystyle\mathcal{F}_{\mathbf{x}}[y](\mathbf{x}):=\sum_{i=1}^{L}c_{i}(\mathbf{x})\nabla_{\boldsymbol{\alpha}_{i}}y(\mathbf{x}),\quad\text{with}\quad\nabla_{\boldsymbol{\alpha}_{i}}y:=\frac{\partial^{\alpha_{i,1}}}{\partial{x}_{1}^{\alpha_{i,1}}}\cdots\frac{\partial^{\alpha_{i,d}}}{\partial{x}_{d}^{\alpha_{i,d}}}y,$

    $\boldsymbol{\alpha}_{i}=\left[\alpha_{i,1},\cdots,\alpha_{i,d}\right]$, and $\mathbf{x}=[x_{1},\cdots,x_{d}]$.

    Here, $L$ is the number of derivatives, $\nabla_{\boldsymbol{\alpha}_{i}}y,\;i=1,\cdots,L$ are the $L$ derivatives in the differential operator $\mathcal{F}_{\mathbf{x}}[y]$, with each a $d$-dimensional multiple index $\boldsymbol{\alpha}_{i}\in\mathbb{N}^{d}$ indicating the order of derivative for each input dimension, and $c_{i}(\mathbf{x})$ is the corresponding coefficient at $\mathbf{x}\in\mathcal{X}$.

    We refer to the corresponding derivative process $y^{\mathcal{F}}(\mathbf{x}):=\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})$ as the PDE process, and our goal is to find the posterior distribution of $y(\mathbf{x}_{\rm new})$ given both measurement data $\mathbf{y}_{1:n}$ and the PDE process $y^{\mathcal{F}}(\mathbf{x})$. To this end, we introduce a set of pseudo points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$, named PDE points, and incorporate the PDE information only at those $m$ points.
    A design method for $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ will be discussed later in Section 3.
    Suppose $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ are selected, and denote the response vector of the PDE process $\mathbf{y}^{\mathcal{F}}_{1:m}=[b(\mathbf{x}^{\mathcal{F}}_{1}),\cdots,b(\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$.
    We further let $\mathcal{I}[\cdot]$ be the identity operator with $\mathcal{I}_{\mathbf{x}}[y](\mathbf{x})=y(\mathbf{x})$, and then the derivatives of the correlation function and mean function are denoted as

    (6)

    $\displaystyle R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\prime})$
    $\displaystyle=(\mathcal{F}_{\mathbf{x}}\times\mathcal{F}_{\mathbf{x}^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime}):=\mathcal{F}_{\mathbf{x}}[\mathcal{F}_{\mathbf{x}^{\prime}}[R_{\boldsymbol{\theta}}]](\mathbf{x},\mathbf{x}^{\prime}),$

    (7)

    $\displaystyle R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\prime})$
    $\displaystyle=(\mathcal{I}_{\mathbf{x}}\times\mathcal{F}_{\mathbf{x}^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime}):=\mathcal{I}_{\mathbf{x}}[\mathcal{F}_{\mathbf{x}^{\prime}}[R_{\boldsymbol{\theta}}]](\mathbf{x},\mathbf{x}^{\prime}),$

    (8)

    $\displaystyle\text{and}\quad\mu^{\mathcal{F}}(\mathbf{x})$
    $\displaystyle=\mathbf{p}^{\top}_{\mathcal{F}}(\mathbf{x})\boldsymbol{\beta}.$

    Here, $\mathbf{p}_{\mathcal{F}}(\mathbf{x})=[\mathcal{F}_{\mathbf{x}}[p_{1}](\mathbf{x}),\cdots,\mathcal{F}_{\mathbf{x}}[p_{q}](\mathbf{x})]^{\top}$ contains the $q$ derivatives of the basis functions, and notation $\mathcal{F}_{\mathbf{x}}\times\mathcal{F}_{\mathbf{x}^{\prime}}$ stands for taking two differential operators $\mathcal{F}_{\mathbf{x}}[\cdot]$ and $\mathcal{F}_{\mathbf{x}^{\prime}}[\cdot]$ on the two inputs of the correlation function $R_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$, respectively. Note that we select the basis functions $\{p_{i}(\cdot)\}_{i=1}^{q}$ and the correlation function $R_{\boldsymbol{\theta}}(\cdot,\cdot)$ such that the above derivatives exist. Since both covariance operator and differential operator are linear, we have the following proposition for the response vector of measurement points $\mathbf{y}_{1:n}$ and that of PDE points $\mathbf{y}^{\mathcal{F}}_{1:m}$.

    ######
    Proposition 2.3 (Theorem 2.2.2 in [1]).

    Suppose that the Gaussian process model (1)-(2) holds with known model parameters $\{\boldsymbol{\beta},\sigma^{2},\boldsymbol{\theta}\}$, and the differential operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ has a linear form (5). Then, response vectors $\mathbf{y}_{1:n}$ and $\mathbf{y}^{\mathcal{F}}_{1:m}$ are jointly normal distributed:

    (9)

    $\displaystyle\begin{bmatrix}\mathbf{y}_{1:n}\vspace{0.1cm}\\
    \mathbf{y}^{\mathcal{F}}_{1:m}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{P}\vspace{0.1cm}\\
    \mathbf{P}_{\mathcal{F}}\end{bmatrix}\boldsymbol{\beta},\;\sigma^{2}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}\right),\quad\text{with}\quad\boldsymbol{\Gamma}_{\boldsymbol{\theta}}=\begin{bmatrix}\mathbf{R}_{\mathcal{I}\mathcal{I}}&\mathbf{R}_{\mathcal{I}\mathcal{F}}\vspace{0.1cm}\\
    \mathbf{R}_{\mathcal{I}\mathcal{F}}^{\top}&\mathbf{R}_{\mathcal{F}\mathcal{F}}\end{bmatrix}.$

    Here, $\mathbf{P}_{\mathcal{F}}=[\mathbf{p}_{\mathcal{F}}(\mathbf{x}^{\mathcal{F}}_{1}),\cdots,\mathbf{p}_{\mathcal{F}}(\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$,
    $\mathbf{R}_{\mathcal{I}\mathcal{I}}=\mathbf{R}_{\boldsymbol{\theta}}$, $\mathbf{R}_{\mathcal{I}\mathcal{F}}={[R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x}_{i},\mathbf{x}^{\mathcal{F}}_{j})]_{i=1}^{n}}_{j=1}^{m}$, and $\mathbf{R}_{\mathcal{F}\mathcal{F}}={[R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\mathbf{x}_{i}^{\mathcal{F}},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}$.

    In words, $\mathbf{P}_{\mathcal{F}}$ is the $m$ by $q$ model matrix for PDE points with the derivatives of basis functions from (8), and the $(n+m)$ by $(n+m)$ correlation matrix
    $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ can be written in a 2-by-2 block matrix form, with
    $\mathbf{R}_{\mathcal{I}\mathcal{I}}$ the standard $n$ by $n$ correlation matrix among observations at $\{\mathbf{x}_{i}\}_{i=1}^{n}$ in (3), $\mathbf{R}_{\mathcal{F}\mathcal{F}}$ the $m$ by $m$ correlation matrix among the PDE points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ with the correlation function $R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\cdot,\cdot)$, and $\mathbf{R}_{\mathcal{I}\mathcal{F}}$ the $n$ by $m$ correlation matrix between $\{\mathbf{x}_{i}\}_{i=1}^{n}$ and $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ with the correlation function $R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\cdot,\cdot)$.

    Now, let $\tilde{\mathbf{y}}_{n+m}=[\mathbf{y}_{1:n}^{\top},\mathbf{y}^{\mathcal{F~{}\top}}_{1:m}]^{\top}$ be an $(n+m)$-vector concatenating the response vectors of measurement points and PDE points. Combining
    Propositions 2.1 and 2.3, we have the following theorem for the posterior distribution of the PIK model with linear PDEs.

    Suppose that the Gaussian process model (1)-(2) holds with known model parameters $\{\boldsymbol{\beta},\sigma^{2},\boldsymbol{\theta}\}$, and the differential operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ has a linear form (5).
    Then, given the response vector $\tilde{\mathbf{y}}_{n+m}$, the posterior distribution $y(\mathbf{x}_{\rm new})$ at a new input location $\mathbf{x}_{\rm new}$ is a normal distribution with mean and variance

    (10)

    $\displaystyle\hat{y}(\mathbf{x}_{\rm new})$
    $\displaystyle:=\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]=\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}(\tilde{\mathbf{y}}_{n+m}-\tilde{\mathbf{P}}\boldsymbol{\beta}),$

    (11)

    $\displaystyle s^{2}(\mathbf{x}_{\rm new})$
    $\displaystyle:=\text{\rm Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]=\sigma^{2}(1-\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})),$

    with $\tilde{\mathbf{P}}=[\mathbf{P}^{\top},\mathbf{P}_{\mathcal{F}}^{\top}]^{\top}$ and $\boldsymbol{\gamma}(\mathbf{x})=[\mathbf{r}^{\top}(\mathbf{x}),R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$.

    Here, $\tilde{\mathbf{P}}$ is the $(n+m)$ by $q$ model matrix for measurement points and PDE points, $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is the correlation matrix defined in (9), and $\boldsymbol{\gamma}(\mathbf{x})$ is the corresponding correlation vector. The proof of Proposition 2.3 can be found in Appendix A, which can be derived from the conditional normal formulation.

    With linear PDEs, Proposition 2.3 indicates that the proposed PIK model can elegantly incorporate the PDE information into the Gaussian process probabilistic modeling framework, via a set of PDE points similar to the measurement data.
    Despite more cumbersome notation, both conditional mean (10) and conditional variance (11) enjoy the closed-form expressions, which can be efficiently computed via the measurement data and the PDEs at hand. Similar to the standard kriging method, the conditional mean can be used for prediction, and the conditional variance can be used for constructing point-wise predictive intervals.

    We are also interested in the predictive distribution of the derivative of the physical quantity of interest, i.e., $y^{\mathcal{G}}(\mathbf{x}_{\rm new}):=\mathcal{G}_{\mathbf{x}}[y](\mathbf{x}_{\rm new})$, where $\mathcal{G}_{\mathbf{x}}[\cdot]$ is a linear differential operator (5). We have the following corollary for the posterior prediction.

    ######
    Corollary 2.4.

    With the same model in Proposition 2.3, the posterior distribution of a derivative process $y^{\mathcal{G}}(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}$ at a new input location $\mathbf{x}_{\rm new}$ is a normal distribution with

    (12)

    $\displaystyle\mathbb{E}\left[y^{\mathcal{G}}(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$
    $\displaystyle=\mathbf{p}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}(\tilde{\mathbf{y}}_{n+m}-\tilde{\mathbf{P}}\boldsymbol{\beta}),$

    (13)

    $\displaystyle\text{\rm Var}\left[y^{\mathcal{G}}(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$
    $\displaystyle=\sigma^{2}(r_{\mathcal{G}}(\mathbf{x})-\boldsymbol{\gamma}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}_{\mathcal{G}}(\mathbf{x}_{\rm new})),$

    where $\mathbf{p}_{\mathcal{G}}(\mathbf{x})=[\mathcal{G}_{\mathbf{x}}[p_{1}](\mathbf{x}),\cdots,\mathcal{G}_{\mathbf{x}}[p_{q}](\mathbf{x})]^{\top}$, $\boldsymbol{\gamma}_{\mathcal{G}}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{I}}(\mathbf{x},\mathbf{x}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{I}}(\mathbf{x},\mathbf{x}_{n}),\break R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$, and $r_{\mathcal{G}}(\mathbf{x})=R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{G}}(\mathbf{x},\mathbf{x})$.

    Here, $\mathbf{p}_{\mathcal{G}}(\mathbf{x})$ consists of the derivatives of the $q$ basis functions via $\mathcal{G}_{\mathbf{x}}[\cdot]$, and $\boldsymbol{\gamma}_{\mathcal{G}}(\mathbf{x})$ is the correlation vector of the considered derivative process $\mathcal{G}_{\mathbf{x}}[y]$ and data $\tilde{\mathbf{y}}_{n+m}$; The notation of the derivatives follow from (8).
    The proof of Corollary 2.4 can be found in Appendix B. Corollary 2.4 shows that, given the response vectors, the posterior distribution of a derivative process at a new point also follows a normal distribution. We will use both conditional mean (12) and variance (13) for in model estimation with nonlinear PDEs (see Section 4).

    ###
    2.3 PIK with nonlinear PDEs

    Now we extend the above PIK model for nonlinear PDEs. Here, we consider the following class of nonlinear PDEs.

    ######
    Definition 2.5 (Nonlinear PDE).

    $\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})=b(\mathbf{x})$ is called nonlinear PDE if the differential operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ is a nonlinear differential operator:

    (14)

    $\displaystyle\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})=\sum_{i=1}^{L}c_{i}(\mathbf{x})\prod_{j=1}^{l_{i}}\nabla_{\boldsymbol{\alpha}_{ij}}y(\mathbf{x}),~{}~{}\text{where}~{}~{}\boldsymbol{\alpha}_{ij}=\left[\alpha_{ij,1},\cdots,\alpha_{ij,d}\right]^{\top}~{}\textrm{and}~{}~{}\prod_{i=1}^{L}l_{i}>1.$

    The nonlinear operator in (14) consists of $L$ terms with each containing a product of $l_{i}\geq 1$ derivatives, and at least one of $l_{i}$’s is strictly greater than one, making the operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ nonlinear (see an example below).
    Our goal is to find the posterior distribution of $y(\mathbf{x}_{\rm new})$ given both the measurement data and the derivative process outputs at PDE points $\tilde{\mathbf{y}}_{n+m}$.
    Due to the existence of the nonlinear terms in (14), the PDE process $y^{\mathcal{F}}(\mathbf{x})=\mathcal{F}_{\mathbf{x}}[y](\mathbf{x})$ is not a Gaussian process. This compromises the elegant kriging framework in Proposition 2.3.
    However, by assuming some derivatives are known, the “conditional” differential operator would be still linear.

    One example for nonlinear PDEs is the inviscid Burger’s equation [2], which we will look into later in Section 5.3; It contains the following nonlinear differential operator

    (15)

    $\displaystyle\mathcal{F}_{\mathbf{x}}[y]:=\frac{\partial y}{\partial t}+y\frac{\partial y}{\partial z}.$

    Here, $\mathbf{x}=[t,z]^{\top}$ is the input with dimension $d=2$, and the coefficients are $c_{1}=c_{2}=1$. The $L=3$ considered derivatives are $\frac{\partial y}{\partial t}$, $y$, and $\frac{\partial y}{\partial z}$ with the corresponding orders $\boldsymbol{\alpha}_{11}=[1,0]$, $\boldsymbol{\alpha}_{21}=[0,0]$, and $\boldsymbol{\alpha}_{22}=[0,1]$, respectively. The nonlinearity is due to the interaction term of $y$ and $\frac{\partial y}{\partial z}$. However, assuming the derivative $\frac{\partial y}{\partial z}=f$ is known, the conditional differential operator $\mathcal{F}_{\mathbf{x}}[y]|_{\partial y/\partial z=f}=\frac{\partial y}{\partial t}+fy$ becomes a linear differential operator. We can then obtain the posterior mean and variance similar to the PIK method with linear PDEs as in (10) and (11).

    With the above intuition in mind, we first select a subset of the derivatives $\{\bar{\boldsymbol{\alpha}}_{k}\}_{k=1}^{K}\subset\{\boldsymbol{\alpha}_{ij}\}$ with the minimal cardinality $K<\sum_{i}l_{i}$, such that one can linearize the differential operator by constructing a new conditional operator $\bar{\mathcal{F}}[\cdot]=\mathcal{F}[\cdot]|\partial[\cdot]$ with $\partial=\{\nabla_{\bar{\boldsymbol{\alpha}}_{1}},\cdots,\nabla_{\bar{\boldsymbol{\alpha}}_{K}}\}$. Such a linearization is always possible with the considered nonlinear PDEs (14). Furthermore, denote the values of those derivatives $\bar{\boldsymbol{\alpha}}_{k}$ ($K$ in total) at the PDE points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ as $\mathbf{z}={\{\partial[y](\mathbf{x}_{j}^{\mathcal{F}})\}_{j=1}^{m}}$; We then treat $\mathbf{z}\in\mathbb{R}^{mK}$ as latent variables. The posterior distribution of $\mathbf{z}$ would be learned together with the model parameters via an expectation-maximization method (see Section 4).
    We have the following theorem for the posterior prediction of the PIK model with nonlinear PDEs.

    Suppose that the Gaussian process model (1)-(2) holds with known model parameters $\{\boldsymbol{\beta},\sigma^{2},\boldsymbol{\theta}\}$, and the nonlinear differential operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ has the form (14). Further assume that $\bar{\mathbf{z}}:=\mathbb{E}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}]$ and $\text{\rm Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}]$ are known.
    Then, at a new input location $\mathbf{x_{\rm new}}$, we have the following posterior mean $\hat{y}(\mathbf{x}_{\rm new})$ and variance $s^{2}(\mathbf{x}_{\rm new})$:

    (16)

    $\displaystyle\hat{y}(\mathbf{x}_{\rm new})$
    $\displaystyle=\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{-1}([\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}),$

    (17)

    $\displaystyle s^{2}(\mathbf{x}_{\rm new})$
    $\displaystyle=\sigma^{2}-\sigma^{2}\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\Sigma}_{\mathbf{z}}\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new}).$

    Here, $\tilde{\mathbf{P}}=[\mathbf{P}^{\top},\mathbf{P}_{\bar{\mathcal{F}}}^{\top},\mathbf{P}_{\partial}^{\top}]^{\top}$, $\boldsymbol{\gamma}(\mathbf{x})=[\mathbf{r}^{\top}(\mathbf{x}),R_{\boldsymbol{\theta}}^{\mathcal{I}\bar{\mathcal{F}}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1:m}),R_{\boldsymbol{\theta}}^{\mathcal{I}{\partial}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1:m})]^{\top}$, the 3-by-3 block matrix $\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}=[\mathbf{R}_{\mathcal{I}\mathcal{I}},\mathbf{R}_{\mathcal{I}\bar{\mathcal{F}}},\mathbf{R}_{\mathcal{I}{\partial}};\mathbf{R}_{\mathcal{I}\bar{\mathcal{F}}}^{\top},\mathbf{R}_{\bar{\mathcal{F}}\bar{\mathcal{F}}},\mathbf{R}_{\bar{\mathcal{F}}\partial};\mathbf{R}^{\top}_{\mathcal{I}{\partial}},\mathbf{R}^{\top}_{\bar{\mathcal{F}}\partial},\mathbf{R}_{\partial\partial}]$ with its explicit expression in Appendix C, and $\boldsymbol{\Sigma}_{\mathbf{z}}=\text{\rm diag}(\mathbf{0}_{n+m},\text{\rm Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}])$.

    In words, $\tilde{\mathbf{P}}$ is the $(n+m+mK)\times q$ model matrix, concatenating basis functions for the measurement points, the corresponding basis functions under linearized differential operator $\bar{\mathcal{F}}$ at PDE points, and those under latent derivatives $\partial$ also at PDE points, vector $\boldsymbol{\gamma}(\mathbf{x})$ is the correlation vector combining the same three parts, and $\bar{\boldsymbol{\Gamma}}_{\boldsymbol{\theta}}$ is the $(n+m+mK)\times(n+m+mK)$ corresponding correlation matrix following from (9). The proof of Definition 2.5 can be found in Appendix C; Note that the law of total expectation and the law of total variance are used for the posterior mean and variance, respectively. We have overloaded some of the notation from (10) and (11); The difference should be clear from the content.

    According to Definition 2.5, the computation of the posterior mean with nonlinear PDEs (16) can be viewed as first evaluating the expected value of the latent variable $\bar{\mathbf{z}}$ and then plugging it into the prediction formulation with linear PDEs (10).
    For the posterior variance (17), we notice an additional variance term $\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\Sigma}_{\mathbf{z}}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})$ with the posterior variance of the latent variables $\boldsymbol{\Sigma}_{\mathbf{z}}=\text{diag}(\mathbf{0}_{n+m},\text{Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}])$, compared to that with linear PDEs (11). This can be interpreted as accounting for the additional uncertainty due to the variability of the latent $\mathbf{z}$. Similarly, the conditional mean in (16) can be used for prediction, and the conditional variance in (17) can be used for quantifying predictive uncertainty.

    ###
    2.4 Difference to numerically solving PDEs

    There are two major differences between the proposed PIK method and the numerical methods for solving PDEs, e.g., finite difference methods [41] and finite element methods [22]. First,
    those numerical methods require boundary conditions and initial conditions along with the governing PDEs. In many real applications, however, the boundary conditions are difficult to obtain. Take the previous wafer heating application as an example. The boundary condition of temperature is difficult to specify since the wafer is placed on a rotating platform with unknown temperature and complex heat flux.
    Another example is to understand the flood flow in healthcare applications [31, 5]. While the governing PDEs are known, i.e., the Naiver-Stocks equation, the boundary conditions are extremely difficult to obtain, considering the patient-specific blood vessel geometry and the interaction between blood flow and soft biological tissues.
    In contrast, the proposed PIK method can utilize measurement points together with PDEs, and provide a prediction of the physical field without boundary conditions.

    Second, the PIK method provides a natural way to quantify the uncertainty via point-wise predictive intervals, similar to the standard kriging method.
    This is because PIK assigns a probabilistic model on the target physical field and then finds the posterior prediction conditional on both measurement data and PDEs. In contrast, standard numerical methods typically adopt a finite discretization. Therefore, the solution is deterministically computed from the boundary to the field without underestimating the associated uncertainty. This is particularly true, when only part of the governing PDE system is known. In such cases, standard numerical methods may not provide any results since the physical field is under-determined. A simple example is that one wants to understand a flow field with velocities of two directions while only having one PDE. Nevertheless, our PIK method can still provide a prediction with uncertainty quantification of the physical field via posterior distributions.

    ##
    3 APIK model

    This section presents the APIK model, which actively designs PDE points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ based on the PIK model and the existing measurements $\{\mathbf{x}_{i}\}_{i=1}^{n}$ (see Fig. 1). We first discuss the design criterion and then present a heuristic method for selecting the PDE data size $m$.

    ###
    3.1 Design PDE points

    We propose to select the input locations of the PDE points by minimizing the following integral mean-squared error (IMSE) criterion [36, 37]:

    (18)

    $\displaystyle\{\hat{\mathbf{x}}^{\mathcal{F}}_{j}\}_{j=1}^{m}$
    $\displaystyle=\operatorname*{\arg\!\min}_{\mathbf{x}^{\mathcal{F}}_{1},\cdots,\mathbf{x}^{\mathcal{F}}_{m}}\int_{\mathcal{X}}\text{Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]\;d\mathbf{x}_{\rm new}.$

    IMSE design criterion (18) can be interpreted as follows. The term $\text{Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$ quantifies the posterior variance (i.e., mean-squared error, MSE) of the physical quantity at an untested input $\mathbf{x}_{\rm new}$, given both measurement data and the potential PDE points.
    An integral is then taken to find the average predictive uncertainty over the whole input space.
    Finally, we select the PDE points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ by minimizing IMSE (18) to ensure the designed PDE points can yield low predictive uncertainty and better predictive accuracy.

    The selected PDE points by (18) can leverage the PDE information in two aspects. First, we explicitly minimize the predictive uncertainty for the potential PIK models. To this end, the selected PDE points not only improve the predictive accuracy of the learned PIK model but also provide a better quantification of uncertainty with narrow predictive intervals. Second, IMSE criterion (18) considers the predictive uncertainty conditional on both measurement data and PDE data, to explore their correlation.
    This ensures that the selected PDE locations minimize the predictive uncertainty according to the measurement data.
    As a result, our APIK model leverages the PDE information via a carefully designed set of PDE points, and therefore, can be shown to achieve better learning performance.

    ###
    3.2 PDE data size

    We present a heuristic method for selecting the size of the PDE points. For simplicity, we consider only linear PDEs in this subsection. As for the nonlinear PDEs, we suggest using the corresponding linearized PDEs in Section 2.3.
    We first define the variance reduction ratio (VRR) of the process $y(\cdot)$, with input locations $\{\mathbf{x}_{i}\}_{i=1}^{n}$ and the corresponding responses $\mathbf{y}_{1:n}$:

    (19)

    $\displaystyle\text{\rm VRR}(\{\mathbf{x}_{i}\}_{i=1}^{n};y(\cdot)):=\frac{\int_{\mathcal{X}}\text{\rm Var}\left[y(\mathbf{x}_{\rm new})\right]d\mathbf{x}_{\rm new}-\int_{\mathcal{X}}\text{\rm Var}\left[y(\mathbf{x}_{\rm new})|\mathbf{y}_{1:n}\right]d\mathbf{x}_{\rm new}}{\int_{\mathcal{X}}\text{\rm Var}\left[y(\mathbf{x}_{\rm new})\right]d\mathbf{x}_{\rm new}}.$

    Given a correlation function $R_{\boldsymbol{\theta}}(\cdot,\cdot)$ for $y(\cdot)$, we have the following proposition for VRR.

    ######
    Proposition 3.1 (Theorem 1 in [18]).

    Given budget $n$ (i.e., $n$ measurement points),

    (20)

    $\displaystyle\sup_{\mathbf{x}_{i},\cdots,\mathbf{x}_{n}\in\mathcal{X}}\text{\rm VRR}(\{\mathbf{x}_{i}\}_{i=1}^{n};y(\cdot))\leq\frac{\sum_{i=1}^{n}\lambda_{i}}{\sum_{i=1}^{\infty}\lambda_{i}},$

    where $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ are the eigenvalues of the correlation function $R_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$.

    Proposition 3.1 indicates that in the optimal case, the variance reduction due to the $n$ measurement points $\{\mathbf{x}_{i}\}_{i=1}^{n}$ would be the summation of the top $n$ eigenvalues $\sum_{i=1}^{n}\lambda_{i}$. For the PDE points, we have the following theorem describing a similar bound for the derivative process $y^{\mathcal{F}}(\cdot)$ with PDE points.

    Suppose that the Gaussian process model (1)-(2) holds and the differential operator $\mathcal{F}_{\mathbf{x}}[\cdot]$ has a linear form (5). Moreover, the correlation function is sufficiently smooth, i.e., $(\mathcal{F}_{\mathbf{x}}\times\mathcal{F}_{\mathbf{x}^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime})$ exists. Then the correlation function of $y^{\mathcal{F}}(\cdot)$ has the expansion:

    (21)

    $\displaystyle R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\prime})=\sum_{i=1}^{\infty}\xi_{i}\varepsilon_{i}(\mathbf{x})\varepsilon_{i}(\mathbf{x}^{\prime}),$

    where $\xi_{1}\geq\xi_{2}\geq\cdots\geq 0$ are eigenvalues and $\varepsilon_{i}(\cdot)$’s are the corresponding eigenfunctions. Furthermore, given $m$ PDE points and no measurement points,

    (22)

    $\displaystyle\sup_{\mathbf{x}^{\mathcal{F}}_{1},\cdots\mathbf{x}^{\mathcal{F}}_{m}\in\mathcal{X}^{m}}\text{\rm VRR}(\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m};y^{\mathcal{F}}(\cdot))\leq\frac{\sum_{j=1}^{m}\xi_{j}}{\sum_{j=1}^{\infty}\xi_{j}}.$

    Here, instead of the correlation function for the physical quantity $R_{\boldsymbol{\theta}}(\cdot,\cdot)$, we decompose the correlation function for the PDE process $R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\cdot,\cdot)$ by Mercer’s theorem. Similarly, $\{\xi_{j}\}_{j=1}^{\infty}$ is the set of eigenvalues with $\xi_{1}\geq\xi_{2}\geq\cdots\geq 0$ and $\{\varepsilon_{j}(\mathbf{x})\}_{j=1}^{\infty}$ is the set of corresponding orthonormal eigenfunctions in $L^{2}(\mathcal{X})$.
    Proposition 3.1 essentially says that the variance reduction by PDE points can also be bounded by the eigenvalues of the correlation function for the PDE process;
    A formal proof of this theorem can be found in Appendix E.

    According to Proposition 3.1 and Proposition 3.1, the variance reduction due to either measurement points or PDE points can be quantified by the corresponding eigenvalues.
    To ensure that the set of PDE points $\{\mathbf{x}^{\mathcal{F}}_{j}\}_{j=1}^{m}$ does contribute to the modeling while not overwhelm the information of the measurement data $\{\mathbf{x}_{i}\}_{i=1}^{n}$, we propose to select $m$ such that the VRRs of the two parts are similar:

    (23)

    $\displaystyle m^{*}=\operatorname*{\arg\!\min}_{m\in\mathbb{N}}\left[\frac{\sum_{i=1}^{n}\lambda_{i}}{\sum_{i=1}^{\infty}\lambda_{i}}-\frac{\sum_{j=1}^{m}\xi_{j}}{\sum_{j=1}^{\infty}\xi_{j}}\right]^{2}.$

    The proposed selection method is similar to selecting leading components in principal components analysis: We first compute the eigenvalues of the two correlation functions $R_{\boldsymbol{\theta}}(\cdot,\cdot)$ and $R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\cdot,\cdot)$, and then compare the two eigenvalue ratios and find the optimal $m$.

    ##
    4 APIK model estimation

    We present
    an expectation-maximization (EM) method for estimating the model parameters $\boldsymbol{\phi}=\{\boldsymbol{\theta},\boldsymbol{\beta},\sigma^{2}\}$ and the posterior moments of the latent variables $\mathbf{z}$, and actively selecting the PDE points $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$. Specifically, we fit a kriging model with only measurement data, and initialize the posterior distribution of the latent $\mathbf{z}$ by its posterior distribution in (12) and (13). We initialize the PDE points $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$ by an $m$-run Sobol’ sequence [38]. We then iterate the following three steps: (i) Expectation step – given the PDE points $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$, we take the expectation of the full likelihood with respect to the posterior distribution of $\mathbf{z}$; (ii) Maximization step – we optimize the expected likelihood and find the model parameters $\boldsymbol{\phi}$; (iii) Design step – given $\boldsymbol{\phi}$ and $\mathbf{z}$, we update $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$ via minimizing the IMSE criterion. Algorithm 1 summarizes the iterative method, with each step to be discussed in detail in this section.
    Note that the following discussion is for the APIK model with nonlinear PDEs (14); If the considered PDEs are linear (5), we discard the expectation step and iterate over the maximization step and the design step.

    Algorithm 1  LearnAPIK($\mathbf{y}_{1:n},\mathcal{F}_{\mathbf{x}}[\cdot],\text{MaxIter}$): Estimate APIK with nonlinear PDEs

    1:$\bullet$  obtain the size of PDE points $m$ via (23)

    2:$\bullet$  initialize $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}\leftarrow\texttt{Sobol sequence}(m,p)$

    3:$\bullet$  set $b_{j}\leftarrow b(\mathbf{x}_{j}^{\mathcal{F}})$ for $j=1,2,\cdots,m$ and set $\tilde{\mathbf{y}}_{n+m}\leftarrow[\mathbf{y}_{1:n}^{\top},b_{1},\cdots,b_{m}]^{\top}$

    4:$\bullet$  generate Monte Carlo input samples $\{\mathbf{x}^{MC}_{k}\}_{k=1}^{100n}$ in the input domain $\mathcal{X}$

    5:for $k=1,\cdots,\text{MaxIter}$ do

    6:    for $j=1,\cdots,m$ do

    7:         $\bullet$  set $\bar{\mathbf{z}}\leftarrow\mathbb{E}(\mathbf{z}|\tilde{\mathbf{y}}_{n+m})$ by (16)

    8:         $\bullet$  set $\boldsymbol{\Sigma}_{\mathbf{z}}\leftarrow\text{diag}(\mathbf{0}_{n+m},\text{Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m};\boldsymbol{\phi}])$ by (17)

    9:         $\bullet$  set function $Q(\boldsymbol{\theta},\boldsymbol{\beta},\sigma^{2})\leftarrow\mathbb{E}\left[l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})|\tilde{\mathbf{y}}_{n+m},\boldsymbol{\phi}\right]$ by (25) $\triangleright$ Expectation step

    10:         function ENLL($\boldsymbol{\theta}$)

    11:             $\bullet$  obtain $\boldsymbol{\beta}$ and $\sigma^{2}$ by (26) using $\boldsymbol{\theta}$

    12:             $\bullet$  return $Q(\boldsymbol{\theta},\boldsymbol{\beta},\sigma^{2})$

    13:         $\bullet$  optimize $\boldsymbol{\theta}^{*}\leftarrow\min\text{ENLL}(\boldsymbol{\theta})$ by L-BFGS $\triangleright$ Maximization step

    14:         $\bullet$  update $\boldsymbol{\beta}$ and $\sigma^{2}$ by (26) using $\boldsymbol{\theta}^{*}$

    15:         $\bullet$  set $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ by (16), and set $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}\leftarrow\boldsymbol{\Gamma}_{\boldsymbol{\theta}}[-j,-j]$ by deleting the $j$-th row and column

    16:         $\bullet$  compute $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}$ by inverting the matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}$

    17:         function IMSE($\mathbf{x}$)

    18:             $\bullet$  obtain $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$ by (30)

    19:             $\bullet$  compute the predictive variance $s^{2}_{k}\leftarrow s^{2}(\mathbf{x}^{MC}_{k})$ via (17) for $k=1,\cdots,100n$

    20:             $\bullet$  return the average $\sum_{k}(s^{2}_{k})/(100n)$

    21:         $\bullet$  optimize $\mathbf{x}_{j}^{\mathcal{F}}\leftarrow\min\text{IMSE}(\mathbf{x})$ by L-BFGS
    $\triangleright$ Design step

    22:         $\bullet$  update $b_{j}\leftarrow b(\mathbf{x}_{j}^{\mathcal{F}})$

    23:$\bullet$  return $\boldsymbol{\phi}=\{\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2}\}$, $\{\bar{\mathbf{z}},\boldsymbol{\Sigma}_{\mathbf{z}}\}$, and $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$

    ###
    4.1 Expectation step

    Given the PDE points $\{\mathbf{x}_{j}^{\mathcal{F}}\}_{j=1}^{m}$, we want to first find the likelihood for measurement data, PDE process outputs at PDE points, and latent variables $\mathbf{z}$. While the joint distribution of $[\tilde{\mathbf{y}}^{\top}_{n+m},\mathbf{z}^{\top}]^{\top}$ is not Gaussian, we have the full negative log-likelihood $l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})$ via the linearized PDE operator (see detailed derivation in Appendix D):

    (24)

    $\displaystyle\log(\det(\sigma^{2}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}))+\frac{1}{\sigma^{2}}\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\mathbf{z}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\right)^{\top}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\mathbf{z}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\right).$

    Here, $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is the $(n+m+mK)$ by $(n+m+mK)$ correlation matrix, and $\tilde{\mathbf{P}}$ is the $(n+m+mK)$ by $q$ correlation matrix, following from the PIK model with nonlinear PDEs (see Definition 2.5).
    The goal is to compute the expectation of the negative log-likelihood $\mathbb{E}\left[l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})|\tilde{\mathbf{y}}_{n+m},\boldsymbol{\phi}\right]$ with respect to the posterior distribution for $\mathbf{z}|\tilde{\mathbf{y}}_{n+m},\boldsymbol{\phi}$.

    Typically, this expectation should be evaluated via Monte Carlo samples due to the complex likelihood structure – note that the correlation matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is also a function of the latent variables $\mathbf{z}$ through the linearized differential operator $\bar{\mathcal{F}}[\cdot]$. However, the high dimensionality (i.e., $mK$, the number of PDE points by the number of derivatives) of $\mathbf{z}$ would require thousands of Monte Carlo samples for an effective approximation. This can be quite time-consuming in practice.
    To address that, we adopt a partially plugged-in method for evaluating the expectation, which would noticeably speed up the computation.

    Specifically, we propose to approximate the expectation by plugging the posterior mean $\bar{\mathbf{z}}=\mathbb{E}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m};\boldsymbol{\phi}]$ in only the correlation matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ due to the lack of closed-form expression, while computing the full expectation for the remaining part, i.e., the quadratic term. Therefore, the expectation $\mathbb{E}\left[l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})|\tilde{\mathbf{y}}_{n+m},\boldsymbol{\phi}\right]$ can be approximated by

    (25)

    $\displaystyle\log(\det(\sigma^{2}\bar{\boldsymbol{\Gamma}}_{\theta}))+\frac{1}{\sigma^{2}}\Big{[}\big{(}[\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\big{)}^{\top}\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\big{(}[\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\big{)}+\text{tr}\big{(}\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\boldsymbol{\Sigma}_{\mathbf{z}}\big{)}\Big{]}.$

    Here,
    $\bar{\boldsymbol{\Gamma}}_{\theta}$ is the correlation matrix with the expected value $\bar{\mathbf{z}}$ plugged in, and $\boldsymbol{\Sigma}_{\mathbf{z}}=\text{diag}(\mathbf{0}_{n+m},\break\text{Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m};\boldsymbol{\phi}])$ is the posterior covariance matrix already defined in Definition 2.5. Compared to the negative log-likelihood in (24), we notice an additional trace term $\text{tr}(\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\boldsymbol{\Sigma}_{\mathbf{z}})$ in (25) to account for the variability of the hidden $\mathbf{z}$. This is summarized in lines 7 - 9 in Algorithm 1.

    ###
    4.2 Maximization step

    The goal of the maximization step is to obtain the model parameters $\boldsymbol{\phi}=\{\boldsymbol{\theta},\boldsymbol{\beta},\sigma^{2}\}$ by minimizing the expectation $\mathbb{E}\left[l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})|\tilde{\mathbf{y}}_{n+m},\boldsymbol{\phi}\right]$ in (25). Note that this step can be further reduced to the optimization of only the correlation parameters $\boldsymbol{\theta}$, by plugging the following two first-order conditions in the objective function (25)

    (26)

    $\displaystyle\begin{split}\hat{\boldsymbol{\beta}}&=\left(\tilde{\mathbf{P}}^{\top}\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\tilde{\mathbf{P}}\right)^{-1}\tilde{\mathbf{P}}^{\top}\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\left[\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}\right]^{\top},\quad\\
    \widehat{\sigma^{2}}&=\frac{1}{N}\left[\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\hat{\boldsymbol{\beta}}\right)^{\top}\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\bar{\mathbf{z}}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\hat{\boldsymbol{\beta}}\right)+\text{tr}\left(\bar{\boldsymbol{\Gamma}}_{\theta}^{-1}\boldsymbol{\Sigma}_{\mathbf{z}}\right)\right],\end{split}$

    where $N=n+m+mK$. Note that the additional trace term in (25) also appears in the closed-form expression for $\widehat{\sigma^{2}}$.
    Then we treat these two variables operators of $\boldsymbol{\theta}$ and minimize $l(\boldsymbol{\theta},\hat{\boldsymbol{\beta}}(\boldsymbol{\theta}),\widehat{\sigma^{2}}(\boldsymbol{\theta},\hat{\boldsymbol{\beta}}(\boldsymbol{\theta}));\tilde{\mathbf{y}}_{n+m})$ over $\boldsymbol{\theta}$ by line search methods, e.g., the L-BFGS method [26]. This part is summarized in lines 10 - 14 in Algorithm 1.

    ###
    4.3 Design step

    Given the model parameters $\boldsymbol{\phi}=\{\boldsymbol{\theta},\boldsymbol{\beta},\sigma^{2}\}$, we want to update the PDE points via minimizing IMSE criterion (18).
    However, optimization (18) is extremely high-dimensional, containing $m$ PDE points each in $\mathbb{R}^{d}$. We therefore propose to adopt the sequential design strategy [25, 14]. Specifically,
    we update one PDE point $\mathbf{x}^{\mathcal{F}}_{j}$ at each iteration by sequentially minimizing IMSE criterion over one point with other PDE points fixed

    (27)

    $\displaystyle\begin{split}\hat{\mathbf{x}}^{\mathcal{F}}_{j}=\operatorname*{\arg\!\min}_{\mathbf{x}^{\mathcal{F}}_{j}}\int_{\mathcal{X}}\text{Var}\left[y(\mathbf{x}_{\rm new})|\mathbf{y}_{1:n},\mathbf{y}^{\mathcal{F}}_{1:j-1},b(\mathbf{x}^{\mathcal{F}}_{j}),\mathbf{y}^{\mathcal{F}}_{j+1:m}\right]\;d\mathbf{x}_{\rm new}.\end{split}$

    Here, the posterior variance is obtained in (17), and $b(\cdot)$ is the non-homogeneous term of the PDEs, i.e., the output of the PDE process. Note that the computation of objective function (27) involves an evaluation of the inverse correlation matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$ in the posterior variance $\text{Var}[y(\mathbf{x}_{\rm new})|\cdot]$.
    This can be time-consuming
    since we plan to use line search methods, which require many matrix inverse steps for updating one PDE point $\mathbf{x}^{\mathcal{F}}_{j}$. To address that, we have the following theorem for the efficient computation of $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$ in (27).

    The inverse $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$ can be efficiently computed as follows:

    (30)

    $\displaystyle\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}=\left[\begin{array}[]{cc}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}+\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})^{\top}s^{2}_{j}(\mathbf{x}^{\mathcal{F}}_{j})&\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\vspace{0.1cm}\\
    \mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})^{\top}&1/s^{2}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\end{array}\right],$

    where $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}$ is $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ without $j$-th row and column, $s^{2}_{j}(\mathbf{x})=\text{\rm Var}[y(\mathbf{x})|\mathbf{y}_{1:n},\mathbf{y}^{\mathcal{F}}_{1:j-1},\mathbf{y}^{\mathcal{F}}_{j+1:m}]$, and $\mathbf{g}_{j}(\mathbf{x})=\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}\boldsymbol{\gamma}_{-j}(\mathbf{x})/s^{2}_{j}(\mathbf{x})$.

    Here, $s^{2}_{j}(\mathbf{x})$ is the predictive variance given $(n+m-1)$ points by (11), and $\boldsymbol{\gamma}_{-j}(\mathbf{x})$ is the corresponding correlation vector with $n$ observed points and $(m-1)$ PDE points except for the $j$-th PDE point from (16). The proof of this theorem can be found in Appendix F.

    Section 4.3 provides an efficient way to solve design problem (27). Specifically, $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}$ does not depend on the $j$-th point $\mathbf{x}^{\mathcal{F}}_{j}$, and therefore, we can compute it ahead of the optimization iteration and use it in the whole optimization iteration. Combining both sequential strategy (27) and fast computation for the inverse (30), we can significantly speed up the design step.
    This part is summarized in lines 15 - 21 in Algorithm 1.

    ##
    5 Applications

    We present in this section four applications. We first look into two one-dimensional (1D) synthetic examples with a linear PDE and a nonlinear PDE, respectively. We then apply the proposed APIK method to two real-world case studies: Shock wave development in flow dynamics and laser heating process in wafer manufacturing.

    In all four applications, we use the posterior mean $\hat{y}(\mathbf{x}_{\rm new})$ for prediction, i.e., (10) with linear PDEs and (16) with nonlinear PDEs, and construct point-wise $2\sigma$ predictive intervals (PIs) via the posterior mean $\hat{y}(\mathbf{x}_{\rm new})$ and the posterior variance $s^{2}(\mathbf{x}_{\rm new})$, i.e., (11) with linear PDEs and (17) with nonlinear PDEs.
    To evaluate predictive accuracy of $\hat{y}(\mathbf{x}_{\rm new})$, we consider the root mean-squared error (RMSE) metric:

    (31)

    $\displaystyle\text{RMSE}=\sqrt{\sum_{i=1}^{n_{\rm test}}\frac{|y(\mathbf{x}_{i})-\hat{y}(\mathbf{x}_{i})|^{2}}{n_{\rm test}}}.$

    Here, $\{\mathbf{x}_{i}\}_{i=1}^{n_{\rm test}}$ are $n_{\rm test}$ test points.

    One way to quantify the performance of the constructed PIs is to use interval score (IS, [13]). For a $2\sigma$ PI $[\hat{y}_{l},\hat{y}_{u}]$, IS is defined as the width of the PI added by penalties of violation

    (32)

    $\displaystyle\text{IS}=(\hat{y}_{u}-\hat{y}_{l})+\frac{2}{\alpha}(\hat{y}_{l}-y^{*})_{+}+\frac{2}{\alpha}(y^{*}-\hat{y}_{u})_{+},$

    where $y^{*}$ is the ground truth, $(a)_{+}=\textrm{max}(a,0)$, and $\alpha=5\%$ is the corresponding type I error of $2\sigma$ PIs, assuming a Gaussian predictive distribution. We then compute the mean interval score (MIS) over the test set $\{\mathbf{x}_{i}\}_{i=1}^{n_{\rm test}}$ as a performance metric over the whole input space.

    ###
    5.1 1D linear PDE example

    We first illustrate the proposed APIK method with a 1D linear PDE. Suppose the physical quantity of interest has the following mean function

    (33)

    $\displaystyle y(x)=x\sin(11x+2),\quad x\in[0,1],$

    with measurement noise variance $\sigma_{\epsilon}^{2}=0.05^{2}$.
    The measurement points are selected via minimax designs [23].
    The test set $\{\mathbf{x}_{i}\}_{i=1}^{n_{\rm test}}$ contains $n_{\rm test}=500$ equally-spaced points over the input space $\mathcal{X}=[0,1]$. We would compute both RMSE and MIS on $\{\mathbf{x}_{i}\}_{i=1}^{n_{\rm test}}$.

    Further suppose the known PDE $\mathcal{F}_{{x}}[y](x)=b(x)$ is delivered in the form

    (34)

    $\displaystyle\mathcal{F}_{{x}}[y](x)=100y(x)+\frac{\partial^{2}y}{\partial x^{2}}(x),\quad b(x)=22\cos(11x+2),\quad\text{and}\quad x\in[0,1].$

    One can easily check that the ground truth function $y(x)$ is one solution of the known PDE $\mathcal{F}_{{x}}[y](x)=b(x)$, by plugging $y(x)$ into (34). The PDE data size $m$ is selected based on the size of measurement data $n$, using the method discussed in Section 3.2.

    Figure 2:  A 1D example with the linear PDE (34): A comparison of the standard kriging method (b), the PIK method (c), and the APIK method (d). The solid lines mark the predictive means, the vertical lines mark the locations of PDE points, and the shades outline the quantification of predictive uncertainty.

    We consider three kriging-based methods, with Gaussian correlation functions $R_{\theta}({x},{x}^{\prime})=\exp(-({x}-{x}^{\prime})^{2}/\theta)$ and constant mean $\mu({x})=\beta_{0}$. The first is the standard kriging method by R package DiceKriging [35]. The second is the proposed PIK model in Section 2 with equally-spaced PDE points. Note that, with linear PDEs, the PIK method can be viewed as an extension to the model in [43]. Finally, we consider the proposed APIK model in Section 3, where the PDE points are actively designed.

    Fig. 2 illustrates the predictions (solid lines) and uncertainty quantification (shades) with $n=5$ measurement points and $m=7$ PDE points, under the three considered methods. As shown in Fig. 2 (b), the standard kriging method estimates an overly-small scale parameter $\theta$ in the Gaussian correlation function $R_{\theta}(x,x^{\prime})$, leading to a clear tread of “regressing” to the constant mean and large PIs. For PIK with equally-spaced PDE points in Fig. 2 (c), the prediction is smoother since we also incorporate the PDE information (34). However, the predictive accuracy is still not good in some regions (e.g., $[0.4,0.6]$) due to the limited data and measurement noise. In addition, while the PIs by PIK are narrower compared to that by standard kriging, they do not cover the true function in input regions of $[0.75,0.9]$.

    The proposed APIK method, as shown in Fig. 2 (d), performs very well in both prediction and PIs. The posterior mean almost overlaps with the ground truth function, and the PIs cover the true function and are much narrower. This is because APIK leverages the PDE information by selecting PDE points that minimize predictive uncertainty. Particularly, the designed PDE points (vertical lines in Fig. 3 (d)) not only explore the whole input domain but also exploit the locations near the measurement points; Those locations appear to provide the most information in reducing the uncertainty. As a result, our APIK model achieves (RMSE, MIS) of (0.08, 0.26), which is much smaller than (0.29, 1.18) for standard kriging and (0.19, 1.05) for PIK.

    |
    | RMSE
    | MIS

    ---|---|---|---

    $n$
    | $m$
    | Kriging
    | PIK
    | APIK
    | Kriging
    | PIK
    | APIK

    4
    | 6
    | 0.5378
    | 0.4509
    | 0.3244
    | 12.232
    | 6.092
    | 3.6122

    5
    | 7
    | 0.2904
    | 0.1908
    | 0.0842
    | 1.179
    | 1.048
    | 0.2621

    7
    | 10
    | 0.1708
    | 0.1650
    | 0.0448
    | 0.6714
    | 0.2415
    | 0.2396

    10
    | 14
    | 0.0504
    | 0.0575
    | 0.0386
    | 0.2821
    | 0.3012
    | 0.2204

    15
    | 20
    | 0.0492
    | 0.0409
    | 0.0376
    | 0.2213
    | 0.2184
    | 0.2010

    Table 1: A comparison of the predictive performance (i.e., RMSE and MIS) under the three considered methods in the 1D example with linear PDE (34).

    |
    | RMSE
    | MIS

    ---|---|---|---

    $n$
    | $m$
    | Kriging
    | PIK
    | APIK
    | Kriging
    | PIK
    | APIK

    3
    | 5
    | 0.4369
    | 0.1910
    | 0.0671
    | 0.9058
    | 2.7520
    | 0.3521

    4
    | 6
    | 0.0904
    | 0.1019
    | 0.0409
    | 0.5111
    | 0.5995
    | 0.1985

    5
    | 7
    | 0.0565
    | 0.0787
    | 0.0292
    | 0.2264
    | 0.4806
    | 0.1977

    7
    | 10
    | 0.0477
    | 0.0463
    | 0.0288
    | 0.3974
    | 0.2809
    | 0.1825

    10
    | 14
    | 0.0583
    | 0.0541
    | 0.0253
    | 0.2155
    | 0.7376
    | 0.1960

    Table 2: A comparison of the predictive performance (i.e., RMSE and MIS) under the three considered methods in the 1D example with nonlinear PDE (35).

    We then consider different measurement data sizes $n\in\{4,5,7,10,15\}$. The corresponding sizes for PDE data are $m\in\{6,7,10,14,20\}$ by (23).
    Table 1 (left) shows the RMSEs using the three methods. First, we notice clear improvements of PIK in accuracy compared to standard kriging. This is again because of the incorporation of PDE information. Second, adopting APIK further improves the performance, since APIK not only incorporates the PDE information but also actively designs the PDE points to leverage the PDE information according to measurement points.
    As a result, we notice at least a $20\%$ reduction in RMSE compared to standard kriging.
    Finally, the improvement of APIK is particularly noticeable when the measurement data size is small, which aligns with target small data challenge in engineering applications.
    We see over $70\%$ improvements when data size $n\leq 7$.

    Table 1 (right) shows the MISs using the three considered methods with different data sizes.
    We notice that while PIK struggles in reducing MISs under several measurement sizes (e.g., $n=5$ and $10$), APIK consistently improves the uncertainty quantification performance with smaller MISs.
    This is because APIK actively introduces PDE information by selecting the PDE points with minimal predictive uncertainty. With improvements in both predictive accuracy and uncertainty quantification, the APIK model better incorporates the known PDE information into the kriging framework, and therefore, improves the learning performance with limited measurement data.

    ###
    5.2 1D nonlinear PDE example

    Now, consider a synthetic 1D example with a nonlinear PDE. Suppose the mean physical quantity of interest and the associated PDE are

    (35)

    $\displaystyle y(x)=\sin(5x+0.35),\quad\text{and}\quad y\frac{\partial y}{\partial x}=2.5\sin(10x+0.7),\quad x\in[0,1].$

    Similar to the example in Section 5.1, the measurement noise variance is $\sigma_{\epsilon}^{2}=0.05^{2}$, the measurement points are selected by maximin designs, and the test set $\{\mathbf{x}_{i}\}_{i=1}^{n_{\rm test}}$ contains $n_{\rm test}=500$ equally-spaced points for computing RMSE and MIS.

    Figure 3:  A 1D example with the nonlinear PDE (35): A comparison of the standard kriging method (b), PIK method (c), and the APIK method (d). The solid lines mark the predictive means, the vertical lines mark the locations of PDE points, and the shades outline the quantification of predictive uncertainty.

    We consider the same three kriging-based methods – standard kriging, PIK, and APIK.
    As shown in Fig. 3 (b), standard kriging performs reasonable; however, the prediction does not factorize in the PDE information at hand. The PIK model in Fig. 3 (c), though incorporating the PDE information via equally-spaced PDE points, still struggles in providing improvements. This is because the PDE points in PIK are not specifically selected according to the measurement points.

    As shown in Fig. 3 (d), APIK achieves noticeable improvements in accuracy and PI. The APIK method incorporates PDE information via a design procedure that considers the correlation between PDE points and measurement points. Therefore, the selected PDE points (vertical lines in Fig. 3 (d)) both explore the whole space and exploit the regions where the PDE information is useful based on the measurement points. For example, the region $[0.4,0.45]$ appears to be important to provide gradient information even with a measurement point nearby, whereas the region $[0.5,0.75]$ may not require additional gradient information since it is mostly linear.
    By smartly leveraging PDE information, our APIK model achieves (RMSE, MIS) of (0.041, 0.199), which is much smaller than (0.090, 0.511) for standard kriging and (0.102, 0.600) for PIK.

    Table 2 compares the predictive performance of the three considered methods, with measurement data sizes ranging from $n\in\{3,4,5,7,10\}$ and the corresponding PDE data sizes $m\in\{5,6,7,10,14\}$. For predictive accuracy, APIK achieves the smallest RMSEs with different data sizes. The improvements are particularly noticeable when the measurement sizes are small. For uncertainty quantification, APIK provides consistent improvements against both standard kriging and PIK. This is again because the proposed AIPK method incorporates and properly leverages the PDE information via carefully designed PDE points, and thereby demonstrating better learning performances.

    Figure 4:  Shack wave case study: (a) shows the ground truth velocity field;
    (b), (c), and (d) show the predicted velocity fields by the three considered methods, with maximum projection design, maximin Latin hypercube design, and minimax design, respectively. The black dots mark the measurement prints, and the blue crosses mark the PDE prints.

    ###
    5.3 Shock wave case study

    We consider the nonlinear inviscid Burger’s equation [2]:

    (36)

    $\displaystyle\frac{\partial y(\mathbf{x})}{\partial t}+y(\mathbf{x})\frac{\partial y(\mathbf{x})}{\partial z}=0,\quad\mathbf{x}=[t,z]\in[0,1]^{2}.$

    Here, $y(\mathbf{x})$ is the velocity field under consideration, with the input vector $\mathbf{x}$ containing the temporal and spatial coordinates. Burger’s equation is a fundamental PDE for conservation equations that can develop discontinuities (i.e., shock waves), occurring in various areas of applied mathematics, such as fluid mechanics, nonlinear acoustics, and traffic flow.

    | RMSE
    | MIS

    ---|---|---

    Design
    | Kriging
    | PIK
    | APIK
    | Kriging
    | PIK
    | APIK

    MaxPro
    | 0.0303
    | 0.0233
    | 0.0195
    | 0.3068
    | 0.2310
    | 0.1131

    MmLHS
    | 0.0251
    | 0.0229
    | 0.0193
    | 0.3312
    | 0.2046
    | 0.1266

    mM
    | 0.0389
    | 0.0365
    | 0.0266
    | 0.4185
    | 0.3514
    | 0.2051

    Table 3: A comparison of the predictive performance (i.e., RMSE and MIS) under the three considered methods in the shock wave case study.

    As exampled in Section 2.3, the inviscid Burger’s equation (36) is a nonlinear PDE due to an interaction term $y(\mathbf{x})\frac{\partial y(\mathbf{x})}{\partial z}$.
    The goal is to understand the velocity field $y(\mathbf{x})$ over the whole 2D spatial and temporal domain $[t,z]\in[0,1]^{2}$, with $n=8$ measurement points and $m=15$ PDE points.
    The actual sensor measurement locations are selected via three popular space-filling designs – the maximal projection design [24], the minimax design using clustering [27], and the maximin Latin hypercube design [30],
    thanks to the implementations in the corresponding R packages.

    The detailed setup is as follows. We first assume (unknown in reality, see Section 2.4) an initial condition $y([t=0,z])$, and obtain the whole field $y(\mathbf{x})$ by solving inviscid Burger’s equation (36) via numerical finite difference method [41]. The obtained velocity field $y(\mathbf{x})$ will be the ground truth function. The test set contains $201\times 201$ grid points over the input space $[0,1]^{2}$. We then simulate $n=8$ measurement data with $y_{i}=y(\mathbf{x}_{i})+\epsilon_{i},i=1,\cdots,8$, with additional measurement errors i.i.d. normally distributed $\epsilon_{i}\sim\mathcal{N}(0,1)$. The simulated measurement data $\{y_{i}\}_{i=1}^{n}$ are then used to estimate the standard kriging model. For PIK, we set the PDE points from an $m=15$-run Sobol’ sequence since it provides a simple way of constructing a space-filling design.
    The PIK model is estimated via an EM algorithm, similar to LearnAPIK (i.e., Algorithm 1) yet without the design step.
    Finally, we estimate the APIK model via LearnAPIK.

    Fig. 4 shows the predictive results with the three measurement data designs, using the three considered methods. Table 3 lists the corresponding RMSEs and MISs.
    Compared to standard kriging, PIK slightly improves the predictive accuracy thanks to the additional PDE information. In addition, APIK further improves the learning performance with at least $20\%$ smaller RMSEs and $60\%$ smaller MISs, compare to both baseline methods.
    This is again because our APIK method designs a PDE data set that minimizes the predictive uncertainty. For all three designs shown in Fig. 4 (b-d), the PDE points selected by APIK not only explore the whole input space but also exploit the locations where the PDE information is the most important. For example, the velocity change in the vertical coordinate is larger than that in the horizontal coordinate in this application. The APIK method accordingly samples clusters of PDE points covering a wide range of vertical coordinates yet with similar horizontal coordinates. By properly leveraging the PDE information via the PDE points, APIK achieves noticeable improvements compared to standard kriging and PIK.

    Figure 5:  An illustration of the laser heating process, and the assassinate temperature measurement in the wafer manufacturing case study.

    ###
    5.4 Laser heating case study

    Finally, we look into the laser heating process, which is an important step in the semiconductor wafer manufacturing process [3, 6]. Wafer manufacturing involves processing silicon wafers in a series of refinement stages to be used as circuit chips.
    Among these stages, thermal processing is one of the most important stages
    since it facilitates the necessary chemical reactions and allows for surface oxidation.
    Fig. 5 illustrates the typical thermal processing procedure: A laser beam (in orange) is moved back and forth over a rotating wafer. Here, industrial engineers wish to understand the temperature field of the wafer over the whole heating process, which would provide a better understanding of possible thermal stresses and thereby the quality of final circuit chips.

    During the heating process, the engineers can conduct temperature measurements using the heat measurement gun (also see Fig. 5). However, since the heating process is typically short, only $n=12$ measurement data are collected. Those measurement locations are selected by a minimax design for its flexibility in the design space shape and good predictive performance. Furthermore, the heat transfer process is governed by the Fourier equation:

    (37)

    $\displaystyle\frac{\partial y}{\partial t}-\beta\frac{\partial^{2}y}{\partial z_{1}^{2}}-\beta\frac{\partial^{2}y}{\partial z_{2}^{2}}=0.$

    According to the discussion in Section 3.2, we select $m=25$ PDE points. The detailed procedure is similar to that in Section 5.3. The only difference is that we use COMSOL multiphysics [9], a finite element analysis software, to obtain the temperature profile on the wafer, with no heat flux boundary condition assumed.

    Figure 6:  A comparison of the temperature contours at different time frames for the wafer manufacturing case study: (a) ground truth, (b) prediction by standard kriging, (c) prediction by PIK, and (d) prediction by APIK.

    Fig. 6 shows the four time frames of the ground truth wafer temperature profile during the heating process, and the predicted profiles via the three considered methods. While PIK already shows improvements against standard kriging by incorporating the PDE information, the proposed APIK method further improves the learning performance by exploring the correlation between measurement points and PDE points.
    Adopting a set of PDE points that smartly leverage PDE information, our APIK method achieves (RMSE, MIS) of (0.0308, 0.1173), which is much smaller than (0.0397, 0.1513) for standard kriging and (0.0333, 0.1376) for PIK.
    Furthermore, as shown in Fig. 6, APIK appears to better capture the divergence field of the temperature, which is known to be critical in introducing the thermal stress and affecting wafer quality [12].

    ##
    6 Conclusion

    In this work, we propose a new physics-informed learning method that combines partial differential equations (PDEs) and measurement data. We adopt the popular kriging method, i.e., Gaussian process regression, for its flexibility and closed-form expressions for both predictive mean and uncertainty quantification. We first propose a PDE Informed Kriging model (PIK) to incorporate both linear and nonlinear PDEs via a pre-selected set of PDE points. Specifically, for nonlinear PDEs, we propose introducing latent variables to linearize the PDEs, and therefore, PIK enjoys the elegant posterior prediction framework similar to standard kriging. We then develop an APIK framework to also actively design the PDE points. The selected PDE points not only explore the whole input domain but also exploit locations where the PDE information is complementary to the measurement data. To estimate the APIK model, we present an expectation-maximization algorithm with fast updates, and also provide a heuristic method for the size of PDE points.
    Finally, we demonstrate the effectiveness of the APIK method in two synthetic examples and two real-world case studies.

    ##
    Appendix A Proof for Proposition 2.3

    We derive here the prediction formula for PIK with linear PDEs. Following Proposition 2.3, we have the joint distribution for measurement data $\mathbf{y}_{1:n}$, outputs of the PDE process $\mathbf{y}^{\mathcal{F}}_{1:m}$, and the prediction $y(\mathbf{x}_{\rm new})$ at a new location $\mathbf{x}_{\rm new}$

    (38)

    $\displaystyle\begin{bmatrix}\mathbf{y}_{1:n}\vspace{0.1cm}\\
    \mathbf{y}^{\mathcal{F}}_{1:m}\vspace{0.1cm}\\
    y(\mathbf{x}_{\rm new})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{P}\vspace{0.1cm}\\
    \mathbf{P}_{\mathcal{F}}\vspace{0.1cm}\\
    \mathbf{p}(\mathbf{x}_{\rm new})\end{bmatrix}\boldsymbol{\beta},\;\sigma^{2}\begin{bmatrix}\mathbf{R}_{\mathcal{I}\mathcal{I}}&\mathbf{R}_{\mathcal{I}\mathcal{F}}&\mathbf{r}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{R}_{\mathcal{I}\mathcal{F}}^{\top}&\mathbf{R}_{\mathcal{F}\mathcal{F}}&\mathbf{r}_{\mathcal{F}}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{r}^{\top}(\mathbf{x}_{\rm new})&\mathbf{r}_{\mathcal{F}}^{\top}(\mathbf{x}_{\rm new})&1\end{bmatrix}\right),$

    Here, $\mathbf{r}_{\mathcal{F}}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$ is the correlation vector, matrices $\mathbf{R}_{\mathcal{I}\mathcal{I}}$, $\mathbf{R}_{\mathcal{I}\mathcal{F}}$, and $\mathbf{R}_{\mathcal{F}\mathcal{F}}$ are defined in Proposition 2.3. Using the conditional Gaussian formulation, we have

    (39)

    $\displaystyle\hat{y}(\mathbf{x}_{\rm new})$
    $\displaystyle:=\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]=\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}(\tilde{\mathbf{y}}_{n+m}-\tilde{\mathbf{P}}\boldsymbol{\beta}),$

    (40)

    $\displaystyle s^{2}(\mathbf{x}_{\rm new})$
    $\displaystyle:=\text{\rm Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]=\sigma^{2}(1-\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})),$

    where $\tilde{\mathbf{P}}=[\mathbf{P}^{\top},\mathbf{P}_{\mathcal{F}}^{\top}]^{\top}$ is the model matrix, and $\boldsymbol{\gamma}(\mathbf{x})=[\mathbf{r}^{\top}(\mathbf{x}),\mathbf{r}_{\mathcal{F}}^{\top}(\mathbf{x})]^{\top}$ is the correlation vector, and $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is the corresponding correlation matrix.

    ##
    Appendix B Proof for Corollary 2.4

    The goal is to find the predictive distribution of $y^{\mathcal{G}}(\cdot):=\mathcal{G}_{\mathbf{x}}[y](\cdot)$.
    Similarly to Appendix A, we have the joint distribution of the same three parts

    (41)

    $\displaystyle\begin{bmatrix}\mathbf{y}_{1:n}\vspace{0.1cm}\\
    \mathbf{y}^{\mathcal{F}}_{1:m}\vspace{0.1cm}\\
    y^{\mathcal{G}}(\mathbf{x}_{\rm new})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{P}\vspace{0.1cm}\\
    \mathbf{P}_{\mathcal{F}}\vspace{0.1cm}\\
    \mathbf{p}_{\mathcal{G}}(\mathbf{x}_{\rm new})\end{bmatrix}\boldsymbol{\beta},\;\sigma^{2}\begin{bmatrix}\mathbf{R}_{\mathcal{I}\mathcal{I}}&\mathbf{R}_{\mathcal{I}\mathcal{F}}&\mathbf{r}_{\mathcal{G}}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{R}_{\mathcal{I}\mathcal{F}}^{\top}&\mathbf{R}_{\mathcal{F}\mathcal{F}}&\mathbf{r}_{\mathcal{F}\mathcal{G}}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{r}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})&\mathbf{r}_{\mathcal{F}\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})&r_{\mathcal{G}}(\mathbf{x}_{\rm new})\end{bmatrix}\right),$

    Here, $\mathbf{r}_{\mathcal{G}}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{G}}(\mathbf{x},\mathbf{x}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{G}}(\mathbf{x},\mathbf{x}_{n})]^{\top}$ is the correlation vector for the measurement data, and $\mathbf{r}_{\mathcal{F}\mathcal{G}}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{G}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{G}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$ is the correlation vector for the PDE data, with the notation for derivative correlation functions

    (42)

    $\displaystyle R_{\boldsymbol{\theta}}^{\mathcal{I}\mathcal{G}}(\mathbf{x},\mathbf{x}^{\prime}):=(\mathcal{I}_{\mathbf{x}}\times\mathcal{G}_{z^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime}),\quad R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{G}}(\mathbf{x},\mathbf{x}^{\prime}):=(\mathcal{F}_{\mathbf{x}}\times\mathcal{G}_{z^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime})$

    following from (8). Furthermore, $\sigma^{2}r_{\mathcal{G}}(\mathbf{x})=\sigma^{2}R_{\boldsymbol{\theta}}^{\mathcal{G}\mathcal{G}}(\mathbf{x},\mathbf{x})$ is the prior process variance at $\mathbf{x}$.
    The posterior distribution of the prediction $y^{\mathcal{G}}(\mathbf{x}_{\rm new})$ at a new location $\mathbf{x}_{\rm new}$ is Gaussian with following mean and variance

    (43)

    $\displaystyle\mathbb{E}\left[y^{\mathcal{G}}(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$
    $\displaystyle=\mathbf{p}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}(\tilde{\mathbf{y}}_{n+m}-\tilde{\mathbf{P}}\boldsymbol{\beta}),$

    (44)

    $\displaystyle\text{\rm Var}\left[y^{\mathcal{G}}(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$
    $\displaystyle=\sigma^{2}(r_{\mathcal{G}}(\mathbf{x})-\boldsymbol{\gamma}_{\mathcal{G}}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}_{\mathcal{G}}(\mathbf{x}_{\rm new})).$

    Here, coefficients $\mathbf{p}_{\mathcal{G}}(\mathbf{x})=[\mathcal{G}_{\mathbf{x}}[p_{1}](\mathbf{x}),\cdots,\mathcal{G}_{\mathbf{x}}[p_{q}](\mathbf{x})]^{\top}$ consists of the derivatives of the $q$ basis functions, and $\boldsymbol{\gamma}_{\mathcal{G}}(\mathbf{x})=[\mathbf{r}_{\mathcal{G}}^{\top},\mathbf{r}_{\mathcal{F}\mathcal{G}}^{\top}]^{\top}$ is the correlation vector.

    ##
    Appendix C Proof for Definition 2.5

    Now, we derive the prediction formula for PIK with nonlinear PDEs. As discussed in Section 2.3, we have linearize the nonlinear differential operator $\mathcal{F}[\cdot]$ to a linear differential operator $\bar{\mathcal{F}}[\cdot]$ via a set of derivatives $\partial$ and the corresponding latent variable $\mathbf{z}$. We have the following joint distribution for measurement data $\mathbf{y}_{1:n}$, outputs of the linearized PDE process $\mathbf{y}^{\bar{\mathcal{F}}}_{1:m}$, the latent variable $\mathbf{z}$, and the prediction $y(\mathbf{x}_{\rm new})$ at a new location $\mathbf{x}_{\rm new}$

    (45)

    $\displaystyle\begin{bmatrix}\mathbf{y}_{1:n}\vspace{0.1cm}\\
    \mathbf{y}^{\bar{\mathcal{F}}}_{1:m}\vspace{0.1cm}\\
    \mathbf{z}\vspace{0.1cm}\\
    y(\mathbf{x}_{\rm new})\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}\mathbf{P}\vspace{0.1cm}\\
    \mathbf{P}_{\bar{\mathcal{F}}}\vspace{0.1cm}\\
    \mathbf{P}_{\partial}\vspace{0.1cm}\\
    \mathbf{p}(\mathbf{x}_{\rm new})\end{bmatrix}\boldsymbol{\beta},\;\sigma^{2}\begin{bmatrix}\mathbf{R}_{\mathcal{I}\mathcal{I}}&\mathbf{R}_{\mathcal{I}\bar{\mathcal{F}}}&\mathbf{R}_{\mathcal{I}\partial}&\mathbf{r}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{R}_{\mathcal{I}\bar{\mathcal{F}}}^{\top}&\mathbf{R}_{\bar{\mathcal{F}}\bar{\mathcal{F}}}&\mathbf{R}_{\bar{\mathcal{F}}\partial}&\mathbf{r}_{\bar{\mathcal{F}}}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{R}_{\mathcal{I}\partial}^{\top}&\mathbf{R}_{\bar{\mathcal{F}}\partial}^{\top}&\mathbf{R}_{\partial\partial}&\mathbf{r}_{\partial}(\mathbf{x}_{\rm new})\vspace{0.1cm}\\
    \mathbf{r}^{\top}(\mathbf{x}_{\rm new})&\mathbf{r}_{\bar{\mathcal{F}}}^{\top}(\mathbf{x}_{\rm new})&\mathbf{r}_{\partial}^{\top}(\mathbf{x}_{\rm new})&1\end{bmatrix}\right).$

    Here, $\mathbf{r}_{\partial}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{I}\partial}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{I}\partial}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$ is the correlation vector for the latent variable $\mathbf{z}$, $\mathbf{r}_{\bar{\mathcal{F}}}(\mathbf{x})=[R_{\boldsymbol{\theta}}^{\mathcal{I}\bar{\mathcal{F}}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{1}),\cdots,R_{\boldsymbol{\theta}}^{\mathcal{I}\bar{\mathcal{F}}}(\mathbf{x},\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$ is the correlation vector for the outputs of the linearized PDE process $\bar{\mathcal{F}}[\cdot]$, and $\mathbf{P}_{\partial}=[\mathbf{p}_{\partial}(\mathbf{x}^{\mathcal{F}}_{1}),\cdots,\mathbf{p}_{\partial}(\mathbf{x}^{\mathcal{F}}_{m})]^{\top}$ is the model matrix, with $\mathbf{p}_{\partial}(\mathbf{x})=[\partial[p_{1}](\mathbf{x}),\cdots,\partial[p_{q}](\mathbf{x})]^{\top}$ containing the $q$ derivatives of the basis functions. Note that, for notation simplicity, we consider $K=|\partial|=1$; For the cases that $K>1$, $\mathbf{r}_{\partial}(\mathbf{x})$ and $\mathbf{P}_{\partial}$ would be the concatenation of all the derivatives in the set $\partial$. Furthermore, the correlation matrices take the following forms:

    (46)

    $\displaystyle\mathbf{R}_{\mathcal{I}\bar{\mathcal{F}}}={[R_{\boldsymbol{\theta}}^{\mathcal{I}\bar{\mathcal{F}}}(\mathbf{x}_{i},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}\quad$
    $\displaystyle\text{with}\quad R_{\boldsymbol{\theta}}^{\mathcal{I}\bar{\mathcal{F}}}(\mathbf{x},\mathbf{x}^{\prime}):=(\mathcal{I}_{\mathbf{x}}\times\bar{\mathcal{F}}_{\mathbf{x}^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime});$

    (47)

    $\displaystyle\mathbf{R}_{\bar{\mathcal{F}}\bar{\mathcal{F}}}={[R_{\boldsymbol{\theta}}^{\bar{\mathcal{F}}\bar{\mathcal{F}}}(\mathbf{x}_{i},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}\quad$
    $\displaystyle\text{with}\quad R_{\boldsymbol{\theta}}^{\bar{\mathcal{F}}\bar{\mathcal{F}}}(\mathbf{x},\mathbf{x}^{\prime}):=(\bar{\mathcal{F}}\times\bar{\mathcal{F}}_{\mathbf{x}^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime});$

    (48)

    $\displaystyle\mathbf{R}_{\mathcal{I}\partial}={[R_{\boldsymbol{\theta}}^{\mathcal{I}\partial}(\mathbf{x}_{i},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}\quad$
    $\displaystyle\text{with}\quad R_{\boldsymbol{\theta}}^{\mathcal{I}\partial}(\mathbf{x},\mathbf{x}^{\prime}):=(\mathcal{I}_{\mathbf{x}}\times\partial_{z^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime});$

    (49)

    $\displaystyle\mathbf{R}_{\bar{\mathcal{F}}\partial}={[R_{\boldsymbol{\theta}}^{\bar{\mathcal{F}}\partial}(\mathbf{x}_{i},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}\quad$
    $\displaystyle\text{with}\quad R_{\boldsymbol{\theta}}^{\bar{\mathcal{F}}\partial}(\mathbf{x},\mathbf{x}^{\prime}):=(\bar{\mathcal{F}}_{\mathbf{x}}\times\partial_{z^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime});$

    (50)

    $\displaystyle\mathbf{R}_{\partial\partial}={[R_{\boldsymbol{\theta}}^{\partial\partial}(\mathbf{x}_{i},\mathbf{x}_{j}^{\mathcal{F}})]_{i=1}^{n}}_{j=1}^{m}\quad$
    $\displaystyle\text{with}\quad R_{\boldsymbol{\theta}}^{\partial\partial}(\mathbf{x},\mathbf{x}^{\prime}):=(\partial_{z}\times\partial_{z^{\prime}})[R_{\boldsymbol{\theta}}](\mathbf{x},\mathbf{x}^{\prime}).$

    Then the posterior mean of the prediction $y(\mathbf{x}_{\rm new})$ at a new input location $\mathbf{x}_{\rm new}$, conditional on both the measurement data $\mathbf{y}_{1:n}$ and the outputs of the linearized PDE process $\mathbf{y}^{\bar{\mathcal{F}}}_{1:m}$ can be evaluated by the law of total expectation:

    (51)

    $\displaystyle\hat{y}(\mathbf{x}_{\rm new})$
    $\displaystyle=\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]=\mathbb{E}\left[\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m},\mathbf{Z}\right]|\tilde{\mathbf{y}}_{n+m}\right]$

    (52)

    $\displaystyle=\mathbb{E}\left[\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}([\tilde{\mathbf{y}}_{n+m}^{\top},\mathbf{z}^{\top}]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta})|\tilde{\mathbf{y}}_{n+m}\right]$

    (53)

    $\displaystyle=\mathbf{p}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\beta}+\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}([\tilde{\mathbf{y}}^{\top}_{n+m},\mathbb{E}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}]^{\top}]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}).$

    Similarly, the posterior variance of the prediction $y(\mathbf{x}_{\rm new})$ at a new input location $\mathbf{x}_{\rm new}$, conditional on $\mathbf{y}_{1:n}$ and $\mathbf{y}^{\bar{\mathcal{F}}}_{1:m}$ is

    (54)

    $\displaystyle s^{2}(\mathbf{x}_{\rm new})$
    $\displaystyle=\text{Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m}\right]$

    (55)

    $\displaystyle=\mathbb{E}\left[\text{Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m},\mathbf{z}\right]|\tilde{\mathbf{y}}_{n+m}\right]+\text{Var}\left[\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m},\mathbf{z}\right]|\tilde{\mathbf{y}}_{n+m}\right]$

    (56)

    $\displaystyle=\text{Var}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m},\mathbf{z}\right]+\text{Var}\left[\mathbb{E}\left[y(\mathbf{x}_{\rm new})|\tilde{\mathbf{y}}_{n+m},\mathbf{z}\right]|\tilde{\mathbf{y}}_{n+m}\right]$

    (57)

    $\displaystyle=\sigma^{2}-\boldsymbol{\gamma}^{\top}(\mathbf{x}_{\rm new})\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\left[\sigma^{2}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})-\text{diag}(\mathbf{0}_{n+m},\text{Var}[\mathbf{z}|\tilde{\mathbf{y}}_{n+m}])\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\boldsymbol{\gamma}(\mathbf{x}_{\rm new})\right].$

    Note that the law of total variance is used in (55).

    ##
    Appendix D Likelihood function

    Without loss of generality, we set the size of measurement data $n=1$, i.e., $\mathbf{y}_{1:n}=y_{1}$, and the size of the PDE data $m=1$, i.e., $\mathbf{y}^{\mathcal{F}}_{1:m}=w_{1}$. Denote the random variable for the measurement observation $Y$, that for the latent variable $Z$, that for observation with the nonlinear PDE $Y^{\mathcal{F}}$ and that for the corresponding linearized PDE $\bar{Y}$.
    According to the definition of likelihood, we have

    (58)

    $\displaystyle l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})=-\log\{\mathbb{P}[Y=y_{1},Y^{\mathcal{F}}=w_{1},Z=z|\boldsymbol{\phi}]\}$

    Now introduce the linearized PDE with differential operator $\bar{\mathcal{F}}$ and the associated output random variable $\bar{Y}$:

    (59)

    $\displaystyle\mathbb{P}[Y=y_{1},Y^{\mathcal{F}}=w_{1},Z=z|\boldsymbol{\phi}]=\mathbb{P}[Y=y_{1},\bar{Y}=w_{1},Z=z|\boldsymbol{\phi}]$

    Note that thanks to the linearization, the joint probability density function of $[Y,\bar{Y},Z]$ is in Gaussian form with mean $\tilde{\mathbf{P}}\boldsymbol{\beta}$ and the correlation $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$. Note that the correlation $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is a function of the latent value $z$. Finally, we have the negative log-likelihood $l(\boldsymbol{\phi};\tilde{\mathbf{y}}_{n+m},\mathbf{z})$

    (60)

    $\displaystyle\log(\det(\sigma^{2}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}))+\frac{1}{\sigma^{2}}\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\mathbf{z}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\right)^{\top}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}\left(\left[\tilde{\mathbf{y}}^{\top}_{n+m},\mathbf{z}^{\top}\right]^{\top}-\tilde{\mathbf{P}}\boldsymbol{\beta}\right).$

    ##
    Appendix E Proof for Proposition 3.1

    Since the correlation function $R_{\boldsymbol{\theta}}$ is sufficiently smooth and $R^{\mathcal{F}\mathcal{F}}_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$ exists, by Theorem 2.2.2 in [1] and induction, we obtain that $y^{\mathcal{F}}$ also follows a Gaussian process with correlation function $R^{\mathcal{F}\mathcal{F}}_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$.

    Then, by Mercer’s expansion, we have

    $R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\mathbf{x},\mathbf{x}^{\prime})=\sum_{i=1}^{\infty}\xi_{i}\varepsilon_{i}(\mathbf{x})\varepsilon_{i}(\mathbf{x}^{\prime}),$

    where $e_{i}(x)^{\prime}$s are the orthonormal basis, i.e., $\int_{\mathcal{X}}\varepsilon_{i}(\mathbf{x})\varepsilon_{j}(\mathbf{x})d\mathbf{x}=1$ if $i=j$, otherwise $0$.

    Furthermore, by Proposition 3.1, we obtain the desired results.

    ##
    Appendix F Proof for Section 4.3

    The correlation matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ can be written in 2-by-2 block matrix form (with switches in rows and columns)

    (63)

    $\displaystyle\boldsymbol{\Gamma}_{\boldsymbol{\theta}}=\left[\begin{array}[]{cc}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}&\boldsymbol{\gamma}_{-j}\vspace{0.15cm}\\
    \boldsymbol{\gamma}_{-j}^{\top}&d_{j}\end{array}\right].$

    Here, $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}$ is the correlation matrix $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ without the $j$-th row and column, $\boldsymbol{\gamma}_{-j}(\mathbf{x})$ is the corresponding correlation vector with $n$ observed points and $(m-1)$ PDE points except for the $j$-th PDE point from (16), and $d_{j}=R_{\boldsymbol{\theta}}^{\mathcal{F}\mathcal{F}}(\mathbf{x}_{j}^{\mathcal{F}})$ is the value of the derivative of the correlation function at $\mathbf{x}_{j}^{\mathcal{F}}$. Adopting the formulation for the block matrix inversion, we have

    (68)

    $\displaystyle\left[\begin{array}[]{cc}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}&\boldsymbol{\gamma}_{-j}\vspace{0.15cm}\\
    \boldsymbol{\gamma}_{-j}^{\top}&d_{j}\end{array}\right]^{-1}=\left[\begin{array}[]{cc}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}+\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}\boldsymbol{\gamma}_{-j}d^{-}_{j}\boldsymbol{\gamma}_{-j}^{\top}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}&-\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}\boldsymbol{\gamma}_{-j}d^{-}_{j}\vspace{0.15cm}\\
    -\boldsymbol{\gamma}_{-j}^{\top}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}d^{-}_{j}&d^{-}_{j}\end{array}\right].$

    where scalar $d^{-}_{j}=1/(d_{j}-\boldsymbol{\gamma}_{-j}^{\top}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}\boldsymbol{\gamma}_{-j})$. Note that $d^{-}_{j}=1/s^{2}_{j}(\mathbf{x}_{j}^{\mathcal{F}})$ holds according to the definition of the posterior variance in (11). Denoting vector $\mathbf{g}_{j}(\mathbf{x})=\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}\boldsymbol{\gamma}_{-j}(\mathbf{x})/s^{2}_{j}(\mathbf{x})$, the inverse $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$ can be efficiently computed as follows:

    (71)

    $\displaystyle\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}=\left[\begin{array}[]{cc}\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}+\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})^{\top}s^{2}_{j}(\mathbf{x}^{\mathcal{F}}_{j})&\mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\vspace{0.15cm}\\
    \mathbf{g}_{j}(\mathbf{x}^{\mathcal{F}}_{j})^{\top}&1/s^{2}_{j}(\mathbf{x}^{\mathcal{F}}_{j})\end{array}\right].$

    Note that the computation of $\boldsymbol{\Gamma}_{\boldsymbol{\theta}-j}^{-1}$ can also be obtained in a similar manner from $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}^{-1}$. We leave this part out due to more cumbersome notation.

    ## References

      *
    [1]

    R. J. Adler, The Geometry of Random Fields, vol. 62, SIAM,
    Philadelphia, 2010.

      *
    [2]

    J. M. Burgers, A mathematical model illustrating the theory of
    turbulence, in Advances in Applied Mechanics, vol. 1, Elsevier, 1948,
    pp. 171–199.

      *
    [3]

    J. Chen, S. Mak, V. R. Joseph, and C. Zhang, Adaptive design for
    gaussian process regression under censoring, arXiv preprint
    arXiv:1910.05452, (2019).

      *
    [4]

    J. Chen, S. Mak, V. R. Joseph, and C. Zhang, Function-on-function
    kriging, with applications to three-dimensional printing of aortic tissues,
    Technometrics, (2020), pp. 1–12.

      *
    [5]

    J. Chen, Y. Xie, K. Wang, C. Zhang, M. A. Vannan, B. Wang, and Z. Qian,
    Active image synthesis for efficient labeling, IEEE Transactions on
    Pattern Analysis and Machine Intelligence, to appear (2020).

      *
    [6]

    Z. Chen, S. Mak, and C. F. J. Wu, A hierarchical expected
    improvement method for Bayesian optimization, arXiv preprint
    arXiv:1911.07285, (2019).

      *
    [7]

    C. Currin, T. Mitchell, M. Morris, and D. Ylvisaker, Bayesian
    prediction of deterministic functions, with applications to the design and
    analysis of computer experiments, Journal of the American Statistical
    Association, 86 (1991), pp. 953–963.

      *
    [8]

    S. Da Veiga and A. Marrel, Gaussian process modeling with inequality
    constraints, in Annales de la Faculté des Sciences de Toulouse:
    Mathématiques, vol. 21, 2012, pp. 529–555.

      *
    [9]

    E. J. Dickinson, H. Ekström, and E. Fontes, COMSOL
    Multiphysics®: Finite element software for electrochemical
    analysis. A mini-review, Electrochemistry Communications, 40 (2014),
    pp. 71–74.

      *
    [10]

    P. J. Diggle, J. A. Tawn, and R. A. Moyeed, Model-based
    geostatistics, Journal of the Royal Statistical Society: Series C (Applied
    Statistics), 47 (1998), pp. 299–350.

      *
    [11]

    L. Ding, S. Mak, and C. F. J. Wu, Bdrygp: a new Gaussian process
    model for incorporating boundary information, arXiv preprint
    arXiv:1908.08868, (2019).

      *
    [12]

    M. R. Eslami, R. B. Hetnarski, J. Ignaczak, N. Noda, N. Sumi, and
    Y. Tanigawa, Theory of Elasticity and Thermal Stresses, vol. 197,
    Springer, 2013.

      *
    [13]

    T. Gneiting and A. E. Raftery, Strictly proper scoring rules,
    prediction, and estimation, Journal of the American Statistical Association,
    102 (2007), pp. 359–378.

      *
    [14]

    R. B. Gramacy and D. W. Apley, Local Gaussian process
    approximation for large computer experiments, Journal of Computational and
    Graphical Statistics, 24 (2015), pp. 561–578.

      *
    [15]

    M. Gulian, M. Raissi, P. Perdikaris, and G. Karniadakis, Machine
    learning of space-fractional differential equations, SIAM Journal on
    Scientific Computing, 41 (2019), pp. 2485–2509.

      *
    [16]

    J. Han and L. Zhang, Integrating machine learning with physics-based
    modeling, arXiv preprint arXiv:2006.02619, (2020).

      *
    [17]

    Z.-H. Han, Y. Zhang, C.-X. Song, and K.-S. Zhang, Weighted
    gradient-enhanced kriging for high-dimensional surrogate modeling and design
    optimization, AIAA Journal, 55 (2017), pp. 4330–4346.

      *
    [18]

    O. Harari, D. Bingham, A. Dean, and D. Higdon, Computer experiments:
    Prediction accuracy, sample size and model complexity revisited, Statistica
    Sinica, (2018), pp. 899–919.

      *
    [19]

    X. He and P. Chien, On the instability issue of gradient-enhanced
    Gaussian process emulators for computer experiments, SIAM/ASA Journal on
    Uncertainty Quantification, 6 (2018), pp. 627–644.

      *
    [20]

    J. M. Hernández-Lobato, J. R. Lloyd, and D. Hernández-Lobato,
    Gaussian process conditional copulas with applications to financial time
    series, in Advances in Neural Information Processing Systems, 2013,
    pp. 1736–1744.

      *
    [21]

    C. Jidling, N. Wahlström, A. Wills, and T. B. Schön, Linearly constrained Gaussian processes, in Advances in Neural Information
    Processing Systems, 2017, pp. 1215–1224.

      *
    [22]

    C. Johnson, Numerical Solution of Partial Differential Equations by
    the Finite Element Method, Dover Publications, Incorporated, 2012.

      *
    [23]

    M. E. Johnson, L. M. Moore, and D. Ylvisaker, Minimax and maximin
    distance designs, Journal of Statistical Planning and Inference, 26 (1990),
    pp. 131–148.

      *
    [24]

    V. R. Joseph, E. Gul, and S. Ba, Maximum projection designs for
    computer experiments, Biometrika, 102 (2015), pp. 371–380.

      *
    [25]

    C. Q. Lam, Sequential Adaptive Designs in Computer Experiments for
    Response Surface Model Fit, PhD thesis, The Ohio State University, 2008.

      *
    [26]

    D. C. Liu and J. Nocedal, On the limited memory BFGS method for
    large scale optimization, Mathematical Programming, 45 (1989), pp. 503–528.

      *
    [27]

    S. Mak and V. R. Joseph, Minimax and minimax projection designs
    using clustering, Journal of Computational and Graphical Statistics, 27
    (2018), pp. 166–178.

      *
    [28]

    S. Mak, C.-L. Sung, X. Wang, S. T. Yeh, Y. H. Chang, V. R. Joseph,
    V. Yang, and C. J. Wu, An efficient surrogate model for emulation and
    physics extraction of large eddy simulations, Journal of the American
    Statistical Association, 113 (2018), pp. 1443–1456.

      *
    [29]

    G. Mathéron, Principles of geostatistics, Economic Geology, 58
    (1963), pp. 1246–1266.

      *
    [30]

    M. D. Morris and T. J. Mitchell, Exploratory designs for
    computational experiments, Journal of Statistical Planning and Inference, 43
    (1995), pp. 381–402.

      *
    [31]

    Z. Qian, K. Wang, S. Liu, X. Zhou, V. Rajagopal, C. Meduri, J. R. Kauten,
    Y.-H. Chang, C. Wu, and C. Zhang, Quantitative prediction of
    paravalvular leak in transcatheter aortic valve replacement based on
    tissue-mimicking 3D printing, JACC: Cardiovascular Imaging, 10 (2017),
    pp. 719–731.

      *
    [32]

    M. Raissi, P. Perdikaris, and G. E. Karniadakis, Machine learning of
    linear differential equations using Gaussian processes, Journal of
    Computational Physics, 348 (2017), pp. 683–693.

      *
    [33]

    M. Raissi, P. Perdikaris, and G. E. Karniadakis, Numerical
    Gaussian processes for time-dependent and nonlinear partial differential
    equations, SIAM Journal on Scientific Computing, 40 (2018), pp. A172–A198.

      *
    [34]

    J. Riihimäki and A. Vehtari, Gaussian processes with
    monotonicity information, in Proceedings of the Thirteenth International
    Conference on Artificial Intelligence and Statistics, 2010, pp. 645–652.

      *
    [35]

    O. Roustant, D. Ginsbourger, and Y. Deville, Dicekriging: Kriging
    methods for computer experiments, R Journal, 1 (2010).

      *
    [36]

    J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn, Design and
    analysis of computer experiments, Statistical Science, 4 (1989),
    pp. 409–423.

      *
    [37]

    T. J. Santner, B. J. Williams, and W. I. Notz, The Design and
    Analysis of Computer Experiments, Springer Science & Business Media, 2018.

      *
    [38]

    I. M. Sobol’, On the distribution of points in a cube and the
    approximate evaluation of integrals, Zhurnal Vychislitel’noi Matematiki i
    Matematicheskoi Fiziki, 7 (1967), pp. 784–802.

      *
    [39]

    A. Solin and M. Kok, Know your boundaries: Constraining gaussian
    processes by variational harmonic features, in The 22nd International
    Conference on Artificial Intelligence and Statistics, 2019, pp. 2193–2202.

      *
    [40]

    T. Sottinen and E. Valkeila, Fractional Brownian Motion as a Model
    in Finance, Department of Mathematics, University of Helsinki, 2001.

      *
    [41]

    J. C. Strikwerda, Finite Difference Schemes and Partial Differential
    Equations 2nd., SIAM, 2004.

      *
    [42]

    S. Ulaganathan, I. Couckuyt, T. Dhaene, J. Degroote, and E. Laermans,
    Performance study of gradient-enhanced kriging, Engineering with
    Computers, 32 (2016), pp. 15–34.

      *
    [43]

    X. Wang and J. O. Berger, Estimating shape constrained functions
    using Gaussian processes, SIAM/ASA Journal on Uncertainty Quantification,
    4 (2016), pp. 1–25.

      *
    [44]

    M. W. Wheeler, D. B. Dunson, S. P. Pandalai, B. A. Baker, and A. H.
    Herring, Mechanistic hierarchical gaussian processes, Journal of the
    American Statistical Association, 109 (2014), pp. 894–904.

      *
    [45]

    J. Willard, X. Jia, S. Xu, M. Steinbach, and V. Kumar, Integrating
    physics-based modeling with machine learning: A survey, arXiv preprint
    arXiv:2003.04919, (2020).

      *
    [46]

    C. K. Williams and C. E. Rasmussen, Gaussian Processes for Machine
    Learning, vol. 2, MIT press Cambridge, MA, 2006.