Blending physics with data using an efficient Gaussian process regression with soft inequality and monotonicity constraints

Abstract

In this work, we propose a new Gaussian process (GP) regression framework that enforces the physical constraints in a probabilistic manner. Specifically, we focus on inequality and monotonicity constraints. This GP model is trained by the quantum-inspired Hamiltonian Monte Carlo (QHMC) algorithm, which is an efficient way to sample from a broad class of distributions by allowing a particle to have a random mass matrix with a probability distribution. Integrating the QHMC into the inequality and monotonicity constrained GP regression in the probabilistic sense, our approach enhances the accuracy and reduces the variance in the resulting GP model. Additionally, the probabilistic aspect of the method leads to reduced computational expenses and execution time. Further, we present an adaptive learning algorithm that guides the selection of constraint locations. The accuracy and efficiency of the method are demonstrated in estimating the hyperparameter of high-dimensional GP models under noisy conditions, reconstructing the sparsely observed state of a steady state heat transport problem, and learning a conservative tracer distribution from sparse tracer concentration measurements.

1 Introduction

In many real-world applications, measuring complex systems or evaluating computational models can be time-consuming, costly or computationally intensive. Gaussian process (GP) regression is one of the Bayesian techniques that addresses this problem by building a surrogate model. It is a supervised machine learning framework that has been widely used in regression and classification tasks. A GP can be interpreted as a suitable probability distribution on a set of functions, which can be conditioned on observations using Bayes’ rule (Lange-Hegermann, 2021). GP regression has found applications in various challenging practical problems including multi-target regression problems Nabati et al. (2022), biomedical applications Dürichen et al. (2014), Pimentel et al. (2013), robotics Williams and Rasmussen (2006) and mechanical engineering applications Song et al. (2021), Li et al. (2023), etc. The recent research demonstrate that a GP regression model can make predictions incorporating prior information (kernels) and generate uncertainty measures over predictions. However, prior knowledge often includes physical laws, and using the standard GP regression framework may lead to an unbounded model in which some points can take infeasible values that violate physical laws (Lange-Hegermann, 2021). For example, non-negativity is a requirement for various physical properties such as temperature, density and viscosity (Pensoneault et al., 2020). Incorporating physical information in GP framework can regularize the behaviour of the model and provide more realistic uncertainties, since the approach concurrently evaluates problem data and physical models (Swiler et al., 2020; Ezati et al., 2024).

A significant amount of research has been conducted to incorporate physical information in GP framework, resulting in various techniques and methodologies (Swiler et al., 2020). For example, a probit model for the likelihood of derivative information can be employed to enforce monotonicity constraints (Riihimäki and Vehtari, 2010). Although this approach can also be used to enforce convexity in one dimension, an additional requirement on Hessian is incorporated for higher dimensions (Da Veiga and Marrel, 2012). In (López-Lopera et al., 2022) an additive GP approach is introduced to account for monotonicity constraints. Although posterior sampling step can be challenging, the additive GP framework enables to satisfy the constraints everywhere in the input space, and it is scalable to higher dimensions. The work presented in Gulian et al. (2022) presents a framework in which spectral decomposition covariance kernels and differential equation constraints are used in a co-kriging setup to perform GP regression constrained by boundary value problems. With their inherent advantages, physics-informed GP models that incorporate physical constraints has applications in diverse areas, such as manufacturing Qiang et al. (2023), forecasting in power grids Mao et al. (2020) or urban flooding models Kohanpur et al. (2023), mimicing drivers’ behavior Wang et al. (2021), monitoring intelligent tire systems Barbosa et al. (2021), predicting fuel flow rate Chati and Balakrishnan (2017), designing wind turbines Wilkie and Galasso (2021), etc. Due to their flexibility, physics-informed GP models can be combined with several approaches to enhance the accuracy of model predictions. These works show that integrating physical knowledge into the prediction process provides accurate results.

Enforcing inequality constraints into a GP is typically challenging as the conditional process, subject to these constraints, does not retain the properties of a GP (Maatouk and Bay, 2017). One of the approaches to handle thi s problem is a data augmentation approach in which the inequality constraints are enforced at various locations and approximate samples are drawn from the predictive distribution (Abrahamsen and Benth, 2001), or using a block covariance kernel (Raissi et al., 2017). Implicitly constrained GP regression method proposed in (Salzmann and Urtasun, 2010) shows that the mean prediction of a GP implicitly satisfies linear constraints, if the constraints are satisfied by the training data. A similar approach shows that when we impose linear inequality constraints on a finite set of points in the domain, the resulting process is a compound Gaussian Process with a truncated Gaussian mean (Agrell, 2019). Most of the approaches assume that the inequalities are satisfied on a finite set of input locations. Based on that assumption, the methods approximate the posterior distribution given those constraint input points. The approach introduced in (Da Veiga and Marrel, 2012) is an example of these methods, where maximum likelihood estimation of GP hyperparameters are investigated under the constraint assumptions. In practice, this should also limit the number of constraint points needed for an effective discrete-location approximation. In addition, the method is not efficient on high-dimensional datasets as it takes a large amount of time to train the GP model.

To the best of our knowlege, the first Gaussian method that satisfies certain inequalities at all the input space is proposed by Maatouk and Bay (2017). The GP approximation of the samples are performed in the finite-dimensional space functions, and a rejection sampling method is used for approximating the posterior. The convergence properties of the method is investigated in (Maatouk et al., 2015). Although using the rejection sampling to obtain posterior helps convergence, it might be computationally expensive. Similar to the previous approaches in which a set of inputs satisfy the constraints, this method also suffers from the curse of dimensionality. Later, the truncated Gaussian approach (López-Lopera et al., 2018) extends the framework in (Maatouk and Bay, 2017) to general sets of linear inequalities. Building upon the approaches in (Maatouk and Bay, 2017; Maatouk et al., 2015), the work presented in (López-Lopera et al., 2018) introduces a finite-dimensional approach that incorporates inequalities for both data interpolation and covariance parameter estimation. In this work, the posterior distribution is expressed as a truncated multinormal distribution. The method uses different Markov Chain Monte Carlo (MCMC) methods and exact sampling methods to obtain the posterior distribution. Among the various MCMC sampling techniques including Gibbs, Metropolis-Hastings (MH) and Hamiltonian Monte Carlo (HMC), the results indicate that HMC sampling is the most efficient one. The truncated Gaussian approaches offer several advantages, including the ability to achieve high accuracy and the flexibility in satisfying multiple inequality conditions. However, although those types of methods address the limitations in (Maatouk and Bay, 2017), they might be time consuming particularly in applications with large datasets or high-dimensional spaces.

In this work, we use QHMC algorithm to train the GP model, and enforce the inequality and monotonicity constraints in a probabilistic manner. Our work addresses the computational limitations caused by high dimensions or large datasets. Unlike truncated Gaussian methods in (López-Lopera et al., 2018) for inequality constraints, or additive GP (López-Lopera et al., 2022) with monotonicity constraints, the proposed method can maintain its efficiency on higher dimensions. Further, we adopt an adaptive learning algorithm that selects the constraint locations. The efficiency and accuracy of the QHMC algorithms are demonstrated on inequality and monotonicity constrained problems. Inequality constrained examples include lower and higher dimensional synthetic problems, a conservative tracer distribution from sparse tracer concentration measurements and a three-dimensional heat transfer problem, while monotonicity constrained examples provide lower and higher dimensional synthetic problems. Our contributions can be summarized in three key points: i) QHMC reduces difference between posterior mean and the ground truth, ii) utilizing QHMC in a probabilistic sense decreases variance and uncertainty, and iii) the proposed algorithm is a robust, efficient and flexible method applicable to a wide range of problems. We implemented QHMC sampling in the truncated Gaussian approach to enhance accuracy and efficiency while working with the QHMC algorithm.

2 Gaussian process under inequality constraints

2.1 Standard GP regression framework

Suppose we have a target function represented by values , where are observations at locations . Here, represents -dimensional vectors in the domain . Using the framework provided in Kuss and Rasmussen (2003), we approximate the target function by a GP, denoted as . We can express aswhere is the mean function and is the covariance and mean functions in Equation 1 are defined as

Typically, the standard squared exponential covariance kernel can be used as a kernel function in Equation 2:where is the signal variance, is the Kronecker delta function and is the length-scale. We then assume that the observation includes an additive independent identically distributed (i.i.d.) Gaussian noise term and having zero mean and variance . We denote the hyperparameters in Equation 3 by , and estimate them using the training data. The parameters can be estimated by minimizing the negative marginal log-likelihood Kuss and Rasmussen (2003), Stein (1988), Zhang (2004):

The following section shows how the parameter updates are performed using the QHMC method.

2.2 Quantum-inspired Hamiltonian Monte Carlo

QHMC is an enhanced version of the HMC algorithm that incorporates a random mass matrix for the particles, following a probability distribution. In conventional HMC, the position is represented by the original variables , while Gaussian momentum is represented by auxiliary variables . Utilizing the energy-time uncertainty relation of quantum mechanics, QHMC allows a particle to have a random mass matrix with a probability distribution. Consequently, in addition to the position and momentum variables, a mass variable is introduced within the QHMC framework. Having a third variable offers the advantage of exploring various landscapes in the state-space. As a result, unlike standard HMC or conventional sampling methods such as MH and Gibbs, QHMC can perform well on discontinuous, non-smooth and spiky distributions Barbu and Zhu (2020), Liu and Zhang (2019). In particular, while the performance of HMC and MH sampling degrade when the distribution is ill-conditioned or multi-modal, the performance of QHMC does not have these limitations. Moreover, QHMC maintains its performance with almost zero additional cost of resampling the mass variable. Due to its efficiency and adaptibility, QHMC can easily integrate with other techniques, or be modified to enhance its performance based on specific objectives and applications. For example, stochastic versions of QHMC can yield accurate solutions with increased efficiency, and the approach is applicable to various scenarios involving missing data Liu and Zhang (2019), Kochan et al. (2022).

The quantum nature of QHMC can be understood by considering a one-dimensional harmonic oscillator example provided in Liu and Zhang (2019). Let us consider a ball with a fixed mass attached to a spring at the origin. Assuming is the displacement, the magnitude of the restoring force that pulls back the ball to the origin is , and the ball oscillates around the origin with period . In contrast to standard HMC where the mass is fixed at 1, QHMC incorporates a time-varying mass, allowing the ball to experience acceleration and explore various distribution landscapes. That is, QHMC has the capability to employ a short time period , corresponding to a small mass , to efficiently explore broad but flat regions. Conversely, in spiky regions, it can switch to a larger time period , i.e., larger , to ensure thorough exploration of all corners of the landscape Liu and Zhang (2019).

The implementation of QHMC is straightforward: i) construct a stochastic process for the mass, and at each time , ii) sample from a distribution . Resampling the positive-definite mass matrix is the only additional step to the standard HMC procedure. In practice, assuming that is independent of and , a mass density function with mean and variance can be where is the identity matrix. QHMC framework simulates the following dynamical system:

In this setting, the potential energy function of the QHMC system is , i.e., the negative of marginal log-likelihood. Algorithm 1 summarizes the steps of QHMC sampling, and, here, we consider the location variables in GP model as the position variables in Algorithm 1. The method evolves the QHMC dynamics in Equation 5 to update the locations . In this work, we implement the QHMC method for inequality constrained GP regression in a probabilistic manner.

2.3 Proposed method

Instead of enforcing all constraints strictly, the approach introduced in Pensoneault et al. (2020) minimizes the negative marginal log-likelihood function in Equation 4 while allowing constraint violations with a small probability. For example, for non-negativity constraints, the following requirement is imposed to the problem:where .

In contrast to enforcing the constraint via truncated Gaussian assumption Maatouk and Bay (2017) or performing inference based on the Laplace approximation and expectation propagation Jensen et al. (2013), the proposed method preserves the Gaussian posterior of the standard GP regression. The method uses a slight modification of the existing cost function. Given a model that follows a Gaussian distribution, the constraint in Equation 6 can be re-expressed by the posterior mean and posterior standard deviation:where stands for the posterior mean, is the standard deviation and is the cumulative distribution function of a Gaussian random variable. Following the work in Pensoneault et al. (2020), in this study was set to for demonstration purposes. As a result, , indicating that two standard deviations below the mean is still nonnegative. Then, using Equation 7 the formulation of the optimization problem is given asIn this particular form of the optimization problem, a functional constraint described by Equation 8 is existent. It can be prohibitive or impossible to satisfy this constraint at all points across the entire domain. Therefore, we adopt a strategy where Equation 8 is enforced only on a selected set of constraint points denoted as . The optimization problem can be reformulated aswhere hyperparameters are estimated to enforce bounds. Solving this optimization problem can be very challenging, and hence, in Pensoneault et al. (2020) additional regularization terms are added. Rather than directly solving the optimization problem, this work adopts the soft-QHMC method, which introduces inequality constraints with a high probability (e.g., 95%) by selecting a specific set of constraint points in the domain. Then non-negativity on the posterior GP is enforced at these selected points. The log-likelihood in Equation 4 is minimized using the QHMC algorithm. Leveraging the Bayesian estimation Gelman et al. (2014), we can approximate the posterior distribution by log-likelihood function and prior probability distribution as shown in the following:The QHMC training flow starts with the Bayesian learning shown in Equation 10 this Bayesian learning and proceeds with an MCMC procedure for drawing samples generated by the Bayesian framework. A general sampling procedure at step is given as

The workflow of soft inequality-constrained GP regression with the sampling procedure in Equation 11 is summarized in Algorithm 2, where QHMC sampling (provided in Algorithm 1) is used as a GP training method. In this version of non-negativity enforced GP regression, the constraint points are located where the posterior variance is highest.

2.3.1 Enforcing monotonicity constraints

Monotonicity constraints on a GP can be enforced using the likelihood of derivative observations. After the selection of active constraints, non-negativity constraints are incorporated in the partial derivative, i.e.,where is a vector of latent values. In the soft-constrained GP method, we introduce the non-negativity information in Equation 12 on a set of selected points, and apply the same procedure as in Equation 9.

Since the derivative is also a GP with mean and covariance matrix Riihimäki and Vehtari (2010):the new posterior distribution using the parameters in Equation 13 is given aswhere and denote the predictions at location . The general sampling procedure for the posterior distribution in Equation 14 is the same as in Equation 11.

3 Theoretical analysis of the method

In this section, employing Bayes’ Theorem, we demonstrate how QHMC is capable of producing a steady-state distribution that approximates the actual posterior distribution. Then, we examine the convergence characteristics of the probabilistic approach on the optimization problem outlined in Equation 9.

3.1 Convergence of QHMC training

The study presented in Liu and Zhang (2019) demonstrates that the QHMC framework can effectively capture a correct steady-state distribution that describes the desired posterior distribution via Bayes’ rule. The joint probability density of can be calculated by Bayesian theorem:The conditional distribution in Equation 15 is approximated as follows:Then, using Equation 16, can be written asEquation 17 shows that the marginal steady distribution approaches the true posterior distribution Liu and Zhang (2019).

3.2 Convergence properties of probabilistic approach

In this section, we show that satisfying the constraints on a set of locations in the domain preserves convergence. Recall the following inequality-constrained optimization problem:

Now, it is necessary to demonstrate that the result obtained by using the selected set of input locations as in Equation 18 converge to the value of the regression model’s output. This convergence ensures that probabilistic approach will eventually result in a model that satisfy the desired conditions.

Note that throughout the proof, it is assumed that is finite. The proof can be constructed for the cases whether the domain is countable or uncountable.

First, let us provide distance metrics used throughout the proof. Following the definitions in Guo et al. (2015), letdenote the distance from a point to a set . Then, the distance of two compact sets and can be defined asThen, using the definitions in Equations 21, 22, the Hausddorff distance between and is defined as . Finally, we define a pseudo-metric to describe the distance between two probability distributions and aswhere is the domain specified in Section 3.2.

4 Numerical examples

In this section, we evaluate the performance of the proposed algorithms on various examples including synthetic and real data. The evaluations consider the size and dimension of the datasets. Several versions of QHMC algorithms are introduced and compared depending on the selection of constraint point locations and probabilistic approach.

We represent the constraint-adaptive, hard-constrained approach as QHMCad and its soft-constrained counterpart as QHMCsoftad. Similarly, QHMCvar refers to the method focusing on variance, while QHMCsoftvar corresponds to its soft-constrained version. The combination of the two approaches with hard and soft constraints are denoted by QHMCboth and QHMCsoftboth, respectively. For the sake of comparison, truncated Gaussian algorithms using an HMC sampler (tnHMC) and a QHMC sampler (tnQHMC) for inequality-constrained examples are implemented, while additive GP (additiveGP) algorithm is adapted for monotonicity-constrained examples.

For the synthetic examples, the time and accuracy performances of the algorithms are evaluated while simultaneously changing the dataset size and noise level in the data. Following Pensoneault et al. (2020), as our metric, we calculate the relative error between the posterior mean and the true value of the target function on a set of test points :

We solve the log-likelihood minimization problem in MATLAB, using the GPML package Rasmussen and Nickisch (2010). For the constrained optimization, we use the fmincon from the MATLAB Optimization Toolbox based on the built-in interior point algorithm. Additionally, in order to highlight the advantage of QHMC over HMC, the proposed approach is implemented with using the standard HMC procedure. The relative error, (calculated as in Equation 35), posterior variance and execution time of each version of QHMC and HMC algorithms are presented.

4.1 Inequaltiy constraints

This section provides two synthetic examples and two real-life application examples to demonstrate the effectiveness of QHMC algorithms on inequality constraints. Synthetic examples compare the performance QHMC approach with truncated Gaussian methods for a 2-dimensional and a 10-dimensional problems. For the 2-dimensional example, the primary focus is on enforcing the non-negativity constraints within the GP model. In the case of the 10-dimensional example, we generalize our analysis to satisfy a different inequality constraint, and evaluate the performances of truncated Gaussian, QHMC and soft-QHMC methods. Third example considers conservative transport in a steady-state velocity field in heterogeneous porous media. Despite being a two-dimensional problem, the non-homogeneous structure of the solute concentration introduces complexity and increases the level of difficulty. The last example is a 3-dimensional heat transfer problem in a hallow sphere.

4.1.1 Example 1

Consider the following 2D function under non-negativity constraints:where . In this example, the GP model is trained via QHMC over 20 randomly selected locations.

In Table 1, the comparison between QHMC and HMC algorithms with a dataset size of 200 is presented. The relative error values indicate that QHMC yields approximately more accurate results than HMC, and it achieves this with a shorter processing time. Consequently, QHMC demonstrates both higher accuracy and efficiency compared to HMC.

Figure 1 presents the relative error values of the algorithms with respect to two parameters: the size of the dataset and signal-to-noise ratio (SNR). It can be seen that the most accurate results without adding any noise are provided by QHMCboth and tnQHMC algorithms with around relative error. However, upon introducing the noise to the data and increasing its magnitude, a distinct pattern is observed. The QHMC methods exhibit relative error values of approximately within the SNR range of to . In contrast, the relative error of the truncated Gaussian methods increases to within the same noise range. This pattern demonstrates that QHMC methods can tolerate noise and maintain higher accuracy under these conditions.

FIGURE 1

Further, we compare the time performances of the algorithms in Figure 2 which demonstrates that QHMC methods, especially the probabilistic QHMC approaches can perform much faster than the truncated Gaussian methods. In this simple 2D problem in Equation 36, the presence of noise does not significantly impact the running times of the QHMC algorithms. In contrast, truncated Gaussian algorithms are slower under noisy environment even when the dataset size is small. Additionally, it can be observed in Figure 3 that the QHMC algorithms, especially QHMCvar and QHMCboth are the most robust ones, as their small relative error comes with a small posterior variance. In contrast, the posterior variance values of the truncated Gaussian methods are higher than QHMC posterior variances even when there is no noise, and gets higher along with the relative error (see Figure 1) when the SNR levels increase. Combining all of these experiments, the inference is that QHMC methods achieve higher accuracy within a shorter time frame. Consequently, these methods prove to be more efficient and robust as they can effectively tolerate changes in parameters. Additionally, it is worth noting that a slight improvement is achieved in the performance of truncated Gaussian algorithms by implementing tnQHMC. Based on the numerical results obtained by tnQHMC, it can be concluded that employing tnQHMC not only yields higher accuracy but also saves some computational time compared to tnHMC.

FIGURE 2

FIGURE 3

4.1.2 Example 2

Consider the 10D Ackley function Eriksson and Poloczek (2021) defined as follows:where , , and . We study the performance of the algorithms on the domain while enforcing the function in Equation 37 to be greater than 5. Table 2 shows the comparison between QHMC and HMC algorithms with 200 data points. Similar to the previous example, the results indicate that QHMC yields approximately 20% more accurate results than HMC in a shorter amount of time.

Figure 4 illustrates that QHMCboth, QHMCsoftboth and truncated Gaussian algorithms yield the lowest error when there is no noise in the data. However, as the noise level increases, truncated Gaussian methods fall behind all QHMC approaches. Specifically, both the QHMCboth and QHMCsofthboth algorithms demonstrate the ability to tolerate noise levels up to with an associated relative error of approximately . However, other variants of QHMC methods display greater noise tolerance when dealing with larger datasets. With fewer than 100 data points, the error rate reaches around , but it decreases to when the number of data points exceeds 100.

FIGURE 4

Figure 5 illustrates the time comparison of the algorithms, where QHMC methods provide around time efficiency for the datasets larger than a size of 150. Combining this time advantage with the higher accuracy of QHMC indicates that both soft and hard constrained QHMC algorithms outperform truncated Gaussian methods across various criteria. QHMC methods offer the flexibility to employ one of the algorithms depending on the priority of the experiments. For example, if speed is the primary consideration, QHMCsoftvar is the fastest method while maintaining a good level of accuracy. If accuracy is the most important metric, employing QHMCboth would be a wiser choice, as it still offers significant time savings compared to other methods.

FIGURE 5

Figure 6 presents that the posterior variance values of truncated Gaussian methods are significantly higher than that of the QHMC algorithms, especially when the noise levels are higher than . As expected, QHMCvar and QHMCsoftvar algorithms offer the lowest variance, while QHMCboth and QHMCsoftboth follow them. A clear pattern is shown in the figure, in which QHMC approaches can tolerate higher noise levels especially when the dataset is large. It is notable that our method demonstrates a significant increase in efficiency as the dimension increases. When comparing this 10D example to the 2D case, the execution times of the truncated Gaussian methods are notably impacted by the dimension, even in the absence of noise in the datasets. Although their relative error levels can remain low without noise, it takes 1.5 times longer than the QHMC methods to offer those accuracy. Additionally, this observation holds only for cases where the data is noise-free. As soon as noise is present, the accuracy of truncated Gaussian methods deteriorates, whereas QHMC methods can withstand the noise and yield good results in a shorter time span. In all tables, bold values indicate the best performance in each metric.

FIGURE 6

4.1.3 Example 3: solute transport in heterogeneous porous media

Following the example in Yang et al. (2019), we examine conservative transport within a constant velocity field in heterogeneous porous media. Let us denote the solute concentration by , and suppose that the measurements of are available at various locations at different times. Conservation laws can be used to describe the processes of flow and transport. Specifically, Darcy flow equation describes the flow by Yang et al. (2019).where is the hydraulic head, is the simulation domain with and , and are known boundary head values and is the unknown hydraulic conductivity field. The field is represented as a stochastic process, with the distribution of values described by a log-normal distribution. Specifically, it is expressed as , where is a second-order stationary GP with a known exponential covariance function, where variance and correlation length . The solute transport by the advection-dispersion equation Emmanuel and Berkowitz (2005), Lin and Tartakovsky (2009), Yang et al. (2019) can be described by

In this context, represents the solute concentration defined over , denotes the fluid velocity given by with being porosity; is the diffusion coefficient, stands for the tortuosity, and is the dispersivity tensor, with diagonal components and . In this study, the transport parameters for Equations 38, 39 are defined as: and . Lastly, the solute is instantaneously injected at at with the intensity Yang et al. (2019). In Figure 7, the ground truth with observation locations and constraint locations are presented to provide an insight into the structure of solute concentration.

FIGURE 7

Table 3 presents a comparison of all versions of QHMC and HMC methods, along with the truncated Gaussian algorithms. Similar to the results observed with synthetic examples, the QHMCboth, QHMCsoftboth, and tnQHMC algorithms demonstrate the most accurate predictions with a relative error of . Notably, QHMCsoftboth emerges as the fastest among the methods while achieving higher accuracy. For instance, the error value obtained by QHMCsoftboth is 0.14, whereas tnQHMC’s error is 0.15. However, QHMCsoftboth delivers a time efficiency gain with slightly superior accuracy. In Figure 8, a comprehensive comparison of the algorithms is presented. The decrease in relative error values is noticeable as constraints are gradually added, following the adopted adaptive approach. Initially, the error is 0.5 and gradually decreases to approximately 0.13. Furthermore, it is evident that the QHMCboth and QHMCsoftboth methods consistently deliver the most accurate results at each step, whereas the performance of the QHMCsoftvar method is outperformed by other approaches.

FIGURE 8

4.1.4 Example 4: heat transfer in a hallow sphere

This 3-dimensional example considers a heat transfer problem in a hallow sphere. Let represent a ball centered at 0 with radius . Defining the hallow sphere as , the equations are given as Yang et al. (2021).

In this context, denotes the normal vector pointing outward, while and represent the azimuthal and elevation angles, respectively, of points within the sphere. We determine the precise heat conductivity in Equation 40 using . The quadratic elements with 12,854 degrees of freedom are employed, and we set to solve the partial differential equations (PDE). Starting with 6 initial locations at 0 on the surface, 6 new constraint locations are introduced based on the active learning approach of the QHMC version. In Figure 9, the decrease is evident in relative error while the constraints are added step by step. In addition, Figure 10 shows the ground truth and the GP result obtained by QHMCsoftboth algorithm, where QHMCsoftboth matches the reference model. The constraint locations of the result are shown in Figure 11. Moreover, its posterior variance is small based on the results shown in Table 2. The table also provides the error, posterior variance and time performances of QHMC and HMC algorithms, where the advantages of QHMC over HMC in all categories, even with the truncated Gaussian algorithm are observed. Although all of the algorithms complete the GP regression in less than 1 min, comparing the truncated Gaussian method with QHMC-based algorithms, time efficiency along with compatible accuracy of QHMC algorithms is achieved. In addition to the time and accuracy performances, it is shown that the posterior variance values are smallest with QHMCvar and QHMCboth approaches, followed by tnQHMC and QHMCad approaches. Using HMC sampling in all methods generates larger posterior variances.

FIGURE 9

FIGURE 10

FIGURE 11

4.2 Monotonicity constraints

This section provides two numerical examples to investigate the effectiveness of our algorithms on monotonicity constraints. As stated in Section 2.3.1, the monotonicity constraints are enforced in the direction of active variables. Similar to the comparisons in previous section, we illustrate the advantages of QHMC over HMC, and then compare the performance of QHMC algorithms with additive GP approach introduced in López-Lopera et al. (2022) with respect to the same criteria as in the previous section.

4.2.1 Example 1

Consider the following 5D function with monotonicity constraints López-Lopera et al. (2022):

In this example, we enforce the function in Equation 41 to be non-decreasing with x ∈ [0, 1]^5. Table 5 shows the performances of HMC and QHMC algorithms, where we observe that QHMC achieves higher accuracy with lower variance in a shorter amount of time. The comparison proves that each version of QHMC is more efficient than HMC In addition, Figure 12 shows the relative error values of QHMC and additive GP algorithms with respect to the change in SNR and dataset size. Based on the results, it is clear that QHMCboth and QHMCsoftboth provide the most accurate results under every different condition, while the difference is more remarkable for the cases in which noise is higher. Although QHMCboth and QHMCsoftboth provides the most accurate results, other QHMC versions also generate more accurate results then additive GP method. Moreover, Figure 13 shows that the soft-constrained QHMC approaches are faster than the hard-constrained QHMC, while hard-constrained QHMC versions are still faster than additive GP algorithm. Additionally, we can see in Figure 14 that QHMC, both hard and soft constrained versions, can reduce the posterior variance. Unlike the additiveGP method, QHMC approaches generate small variances even under high noise.

FIGURE 12

FIGURE 13

FIGURE 14

4.2.2 Example 2

We provide a 20-dimensional example to indicate the applicability and effectiveness of QHMC algorithms on higher dimensions with monotonicity constraint. We consider the target function used in López-Lopera et al. (2022), Bachoc et al. (2022).with . We enforce monotonicity constraint on the 20D function in Equation 42.

Table 6 illustrates accuracy and time advantages of QHMC over HMC. For each version of QHMC and HMC, using QHMC sampling in a specific version accelerates the process while increasing the accuracy. Overall comparison shows that among all versions with QHMC and HMC sampling, QHMCboth is the most accurate approach, while QHMCsoftboth is the fastest and ranked second in accuracy. Figures 15, 16 show the relative error and time performances of QHMC-based algorithms, HMCsoftboth and additive GP algorithm, respectively. In this final example with the highest dimension, the same phenomenon is observed as in previous results: soft-constrained versions demonstrate greater efficiency, while hard-constrained QHMC approaches remain faster than additive GP across different conditions, including high noise levels. Based on Figure 15, QHMCboth can tolerate noise levels up to with the smallest error, and it can still provide good accuracy (error is around 0.15) even when the SNR is higher than . It is also worth to mention that although the error values generated by HMCsoftboth and additiveGP are pretty close, HMCsoftboth performs faster than additiveGP, especially when the dataset is larger and noise level is higher. Furthermore, QHMC reduces the posterior variance as in shown in Figure 17. The behavior of the algorithms follows the same trend as 5-dimensional example, where the methods can tolerate the noise in the data, especially with larger datasets.

FIGURE 15

FIGURE 16

FIGURE 17

4.3 Discussion

5 Conclusion

Leveraging the accuracy of QHMC training and the efficiency of probabilistic approach, this work introduced a soft-constrained QHMC algorithm to enforce inequality and monotonicity constraints on the GP. The proposed algorithm reduces the difference between ground truth and the posterior mean in the resulting GP model, while increasing the efficiency by attaining the accurate results in a short amount of time. To further enhance the performance of the QHMC algorithms across various scenarios, modified versions of QHMC are implemented adopting adaptive learning. These versions provide flexibility in selecting the most suitable algorithm based on the specific priorities of a given problem.

We provided the convergence of QHMC by showing that its steady-state distribution approach the true posterior density, and theoretically justified that the probabilistic approach preserves convergence. Finally, we have implemented our methods to solve several types of optimization problems. Each experiment initially outlined the benefits of QHMC sampling in comparison to HMC sampling. These advantages remained consistent across all cases, resulting in approximately a time-saving and higher accuracy. Having demonstrated the advantages of QHMC sampling, further evaluation on the performances of the algorithms across various scenarios was performed. The examples cover higher-dimensional problems featuring both inequality and monotonicity constraints. Furthermore, the evaluations include real-world applications where injecting physical properties is essential, particularly in cases involving inequality constraints.

In the context of inequality-constrained Gaussian processes (GPs), we explored 2-dimensional and 10-dimensional synthetic problems, along with two real applications involving 2-dimensional and 3-dimensional data. For synthetic examples, the relative error, posterior variance and execution time of the algorithms were compared while gradually increasing the noise level and dataset size. Overall, QHMC-based algorithms outperformed the truncated Gaussian methods. Although the truncated Gaussian methods provide high accuracy in the absence of noise and are compatible with QHMC approaches, their relative error and posterior variances increase as the noise appeared and increased. Moreover, the advantages of soft-constrained QHMC became more evident, particularly in higher-dimensional cases, when compared to truncated Gaussian and even hard-constrained QHMC. The time comparison of the algorithms underscores that the truncated Gaussian methods are significantly impacted by the curse of dimensionality and large datasets, exhibiting slower performance under these conditions. In real-world application scenarios featuring 2-dimensional and 3-dimensional data, the findings were consistent with those observed in the synthetic examples. Although the accuracy level may not reach the highest levels observed in the synthetic examples and 3-dimensional heat equation problem, the observed trend remains consistent. The lower accuracy observed in the latter problem can be attributed to the non-homogeneous structure of solute concentration.

In the case of monotonicity-constrained GP, we addressed 5-dimensional and 20-dimensional examples, utilizing the same configuration as employed for inequality-constrained GP. A comprehensive comparison was conducted between all versions of QHMC algorithms and the additive GP method. The results indicate that QHMC-based approaches hold a notable advantage, particularly in scenarios involving noise and large datasets. While additive GP proves to be a strong method suitable for high-dimensional cases, QHMC algorithms performed faster and yielded lower variances.

In conclusion, the work has demonstrated that soft-constrained QHMC is a robust, efficient and flexible method that can be applicable to higher dimensional cases and large datasets. Numerical results have shown that soft-constrained QHMC is promising to be generalized to various applications with different physical properties.

References

• 1

  AbrahamsenP.BenthF. E. (2001). Kriging with inequality constraints. Math. Geol.33, 719–744. 10.1023/a:1011078716252

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AgrellC. (2019). Gaussian processes with linear operator inequality constraints. arXiv Prepr. arXiv:1901, 03134.

  • Google Scholar
  • View reference in article

• 3

  BachocF.López-LoperaA. F.RoustantO. (2022). Sequential construction and dimension reduction of Gaussian processes under inequality constraints. SIAM J. Math. Data Sci.4, 772–800. 10.1137/21m1407513

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BarbosaB. H. G.XuN.AskariH.KhajepourA. (2021). Lateral force prediction using Gaussian process regression for intelligent tire systems. IEEE Trans. Syst. Man, Cybern. Syst.52, 5332–5343. 10.1109/tsmc.2021.3123310

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BarbuA.ZhuS.-C. (2020) Monte Carlo methods, 35. Springer.

  • Google Scholar
  • View reference in article

• 6

  ChatiY. S.BalakrishnanH. (2017). “A Gaussian process regression approach to model aircraft engine fuel flow rate,” in Proceedings of the 8th international conference on cyber-physical systems, 131–140.

  • Google Scholar
  • View reference in article

• 7

  Da VeigaS.MarrelA. (2012). Gaussian process modeling with inequality constraints. Ann. la Fac. Sci. Toulouse Mathématiques21, 529–555. 10.5802/afst.1344

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  DürichenR.PimentelM. A.CliftonL.SchweikardA.CliftonD. A. (2014). “Multi-task Gaussian process models for biomedical applications,” in IEEE-EMBS international Conference on Biomedical and health informatics (BHI) (IEEE), 492–495.

  • Google Scholar
  • View reference in article

• 9

  EmmanuelS.BerkowitzB. (2005). Mixing-induced precipitation and porosity evolution in porous media. Adv. water Resour.28, 337–344. 10.1016/j.advwatres.2004.11.010

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  ErikssonD.PoloczekM. (2021). Scalable constrained Bayesian optimization, In Proceedings of the 24th International Conference on Artificial Intelligence and Statistics (PMLR 130:730–738, 2021), 730–738606.

  • Google Scholar
  • View reference in article

• 11

  EzatiM.EsmaeilbeigiM.KamandiA. (2024). Novel approaches for hyper-parameter tuning of physics-informed Gaussian processes: application to parametric pdes. Eng. Comput.40, 3175–3194. 10.1007/s00366-024-01970-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  GelmanA.CarlinJ. B.SternH. S.DunsonD. B.VehtariA.RubinD. B. (2014). Bayesian data analysis. Philadelphia, Pennsylvania: Tyler & Francis Group, Inc.

  • Google Scholar
  • View reference in article

• 13

  GulianM.FrankelA.SwilerL. (2022). Gaussian process regression constrained by boundary value problems. Comput. Methods Appl. Mech. Eng.388, 114117. 10.1016/j.cma.2021.114117

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  GuoS.XuH.ZhangL. (2015). Convergence analysis for mathematical programs with distributionally robust chance constraint. SIAM J. Optim.27, 784–816. (to appear). 10.1137/15m1036592

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  HessC. (1999). Conditional expectation and martingales of random sets. Pattern Recognit.32, 1543–1567. 10.1016/s0031-3203(99)00020-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  JensenB. S.NielsenJ. B.LarsenJ. (2013). “Bounded Gaussian process regression,” in 2013 IEEE international workshop on machine learning for signal processing (MLSP) (IEEE), 1–6.

  • Google Scholar
  • View reference in article

• 17

  KochanD.ZhangZ.YangX. (2022). “A quantum-inspired Hamiltonian Monte Carlo method for missing data imputation,” in Mathematical and scientific machine learning (PMLR), 17–32.

  • Google Scholar
  • View reference in article

• 18

  KohanpurA. H.SaksenaS.DeyS.JohnsonJ. M.RiasiM. S.YeghiazarianL.et al (2023). Urban flood modeling: uncertainty quantification and physics-informed Gaussian processes regression forecasting. Water Resour. Res.59, e2022WR033939. 10.1029/2022wr033939

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  KussM.RasmussenC. (2003). Gaussian processes in reinforcement learning. Adv. neural Inf. Process. Syst.16.

  • Google Scholar
  • View reference in article

• 20

  Lange-HegermannM. (2021). Linearly constrained Gaussian processes with boundary conditions. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics (PMLR628 130: 1090–1098).

  • Google Scholar
  • View reference in article

• 21

  LiX.-Q.SongL.-K.ChoyY.-S.BaiG.-C. (2023). Fatigue reliability analysis of aeroengine blade-disc systems using physics-informed ensemble learning. Philosophical Trans. R. Soc. A381, 20220384. 10.1098/rsta.2022.0384

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  LinG.TartakovskyA. M. (2009). An efficient, high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media. Adv. Water Resour.32, 712–722. 10.1016/j.advwatres.2008.09.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LiuZ.ZhangZ. (2019). Quantum-inspired Hamiltonian Monte Carlo for bayesian sampling. arXiv Prepr. arXiv:1912.01937.

  • Google Scholar
  • View reference in article

• 24

  López-LoperaA.BachocF.RoustantO. (2022). High-dimensional additive Gaussian processes under monotonicity constraints. Adv. Neural Inf. Process. Syst.35, 8041–8053.

  • Google Scholar
  • View reference in article

• 25

  López-LoperaA. F.BachocF.DurrandeN.RoustantO. (2018). Finite-dimensional Gaussian approximation with linear inequality constraints. SIAM/ASA J. Uncertain. Quantification6, 1224–1255. 10.1137/17m1153157

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  MaatoukH.BayX. (2017). Gaussian process emulators for computer experiments with inequality constraints. Math. Geosci.49, 557–582. 10.1007/s11004-017-9673-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  MaatoukH.RoustantO.RichetY. (2015). Cross-validation estimations of hyper-parameters of Gaussian processes with inequality constraints. Procedia Environ. Sci.27, 38–44. 10.1016/j.proenv.2015.07.105

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  MaoZ.JagtapA. D.KarniadakisG. E. (2020). Physics-informed neural networks for high-speed flows. Comput. Methods Appl. Mech. Eng.360, 112789. 10.1016/j.cma.2019.112789

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  NabatiM.GhorashiS. A.ShahbazianR. (2022). Jgpr: a computationally efficient multi-target Gaussian process regression algorithm. Mach. Learn.111, 1987–2010. 10.1007/s10994-022-06170-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  PensoneaultA.YangX.ZhuX. (2020). Nonnegativity-enforced Gaussian process regression. Theor. Appl. Mech. Lett.10, 182–187. 10.1016/j.taml.2020.01.036

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  PimentelM. A.CliftonD. A.CliftonL.TarassenkoL. (2013). “Probabilistic estimation of respiratory rate using Gaussian processes,” in 2013 35th annual international conference of the IEEE engineering in medicine and biology society (EMBC) (IEEE), 2902–2905.

  • Google Scholar
  • View reference in article

• 32

  QiangB.ShiK.LiuN.RenJ.ShiY. (2023). Integrating physics-informed recurrent Gaussian process regression into instance transfer for predicting tool wear in milling process. J. Manuf. Syst.68, 42–55. 10.1016/j.jmsy.2023.02.019

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  RaissiM.PerdikarisP.KarniadakisG. E. (2017). Machine learning of linear differential equations using Gaussian processes. J. Comput. Phys.348, 683–693. 10.1016/j.jcp.2017.07.050

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  RasmussenC. E.NickischH. (2010). Gaussian processes for machine learning (gpml) toolbox. J. Mach. Learn. Res.11, 3011–3015.

  • Google Scholar
  • View reference in article

• 35

  RiihimäkiJ.VehtariA. (2010). “Gaussian processes with monotonicity information,” in Proceedings of the thirteenth international conference on artificial intelligence andf statistics (JMLR Workshop and Conference Proceedings), 645–652.

  • Google Scholar
  • View reference in article

• 36

  SalzmannM.UrtasunR. (2010). Implicitly constrained Gaussian process regression for monocular non-rigid pose estimation. Adv. neural Inf. Process. Syst.23.

  • Google Scholar
  • View reference in article

• 37

  SongF.ShiK.LiK. (2021). “A nonlinear finite element-based supervised machine learning approach for efficiently predicting collapse resistance of wireline tool housings subjected to combined loads,” in ASME International Mechanical Engineering Congress and Exposition (American Society of Mechanical Engineers). 10.1115/imece2021-7222285680V012T12A057

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  SteinM. L. (1988). Asymptotically efficient prediction of a random field with a misspecified covariance function. Ann. Statistics16, 55–63. 10.1214/aos/1176350690

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  SwilerL. P.GulianM.FrankelA. L.SaftaC.JakemanJ. D. (2020). A survey of constrained Gaussian process regression: approaches and implementation challenges. J. Mach. Learn. Model. Comput.1, 119–156. 10.1615/jmachlearnmodelcomput.2020035155

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  WangY.WangZ.HanK.TiwariP.WorkD. B. (2021). “Personalized adaptive cruise control via Gaussian process regression,” in 2021 IEEE international intelligent transportation systems conference (ITSC) (IEEE), 1496–1502.

  • Google Scholar
  • View reference in article

• 41

  WilkieD.GalassoC. (2021). Gaussian process regression for fatigue reliability analysis of offshore wind turbines. Struct. Saf.88, 102020. 10.1016/j.strusafe.2020.102020

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  WilliamsC. K.RasmussenC. E. (2006). Gaussian Processes for Machine Learning2. MA: MIT press Cambridge.

  • Google Scholar
  • View reference in article

• 43

  YangX.Barajas-SolanoD.TartakovskyG.TartakovskyA. M. (2019). Physics-informed cokriging: a Gaussian-process-regression-based multifidelity method for data-model convergence. J. Comput. Phys.395, 410–431. 10.1016/j.jcp.2019.06.041

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  YangX.TartakovskyG.TartakovskyA. M. (2021). Physics information aided kriging using stochastic simulation models. SIAM J. Sci. Comput.43, A3862–A3891. 10.1137/20m1331585

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  ZhangH. (2004). Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics. J. Am. Stat. Assoc.99, 250–261. 10.1198/016214504000000241

  • CrossRef
  • Google Scholar
  • View reference in article