Gaussian Process Regression with Soft Equality Constraints

1. Introduction

In many machine learning applications, measuring complex systems or evaluating computational models can be challenging due to the high costs, time requirements, or intensive computational demands. GP regression offers a solution by creating a surrogate model that approximates the behavior of these systems. As a Bayesian approach to supervised learning, GP regression is commonly applied in tasks such as regression and classification. It defines a probability distribution over a set of functions, allowing for predictions based on observed data through Bayesian conditioning [1,2].

A key advantage of GP regression arises from its ability to incorporate prior information through kernels, enabling the model to make predictions while also providing uncertainty estimates [3]. However, in many cases, the prior information reflects underlying physical laws [1]. For example, physical properties such as temperature, density, and viscosity must adhere to non-negativity constraints [4]. The standard GP regression framework may not account for these constraints, potentially leading to predictions that violate fundamental physical principles, resulting in unbounded or infeasible values. Therefore, the integration of physical constraints acts as a form of regularization, improving the accuracy and reliability of a GP model in real-world applications [5].

For example, enforcing inequality constraints in GP regression is challenging, as the conditional process does not retain typical GP properties [6]. Approaches to address these type of problems include data augmentation to enforce constraints at multiple locations [7] or using block covariance kernels [8]. Implicitly constrained GP methods show that linear constraints can be satisfied if met by the training data [9], while other methods demonstrate that imposing linear inequality constraints results in a compound GP with a truncated Gaussian mean [10]. Most approaches assume constraints on a finite set of points, approximating the posterior distribution accordingly [6,11,12]. The study in [13] extends the work in [6] by approximating the posterior distribution by a truncated multinormal distribution, with various Markov chain Monte Carlo (MCMC) methods, including Gibbs, Metropolis-Hastings (MH), and Hamiltonian Monte Carlo (HMC) used for sampling. While truncated Gaussian methods provide high accuracy and the flexibility to satisfy multiple inequality conditions, they can be time-consuming, particularly when dealing with large datasets or high-dimensional problems, despite addressing some limitations identified in [6]. In [14], the QHMC algorithm is employed to train the GP model and enforce inequality and monotonicity constraints in a probabilistic manner. This method addresses the computational challenges associated with high-dimensional problems and large datasets. Unlike the truncated Gaussian techniques used in [13] for handling inequality constraints or the additive GP approach in [15] for imposing monotonicity constraints, the proposed method is capable of maintaining efficiency even in higher-dimensional settings.

Furthermore, constraints may also be more tightly integrated with the underlying physics: the GP can be constrained to satisfy constraints which represent physical laws expressed as partial different equations (PDE) [5]. An earlier method for constraining GP with differential equations involves constructing a specialized covariance kernel so that the GP satisfies the constraint across the entire domain. This technique traces back to the divergence-free kernel [16] for vector-valued GPs. It not only enforces the constraint strictly but also avoids the computational cost associated with the four-block covariance matrix. However, this approach has limited applicability and requires specialized implementation since it involves analytically deriving a kernel with the necessary constraining properties. Later, a method for designing the covariance function of a multivariate GP subject to known linear operator constraints on the target function is presented in [17]. The method will by construction guarantee that any sample drawn from the resulting process will obey the constraints in all points. While there are other methods for constructing the kernel function, like the one proposed in [18], these approaches often involve specific conditions. This can make it more challenging to generalize the methods across different scenarios.

The traditional approaches for solving differential equations has long relied on analytical and numerical methods under well-posed conditions, constrained by the quality of boundary and forcing conditions. However, recent advancements are fundamentally shifting this paradigm by integrating probabilistic machine learning techniques with differential equations. For example, data-driven algorithms that use GP priors tailored for integro-differential operators enable the solving of general linear equations. These algorithms can work with sparse, noisy, multi-fidelity data that do not need to be located at the domain boundaries [19,20]. This approach not only provides predictive posterior distributions to quantify uncertainty but also enables adaptive solution refinement through active learning.

Building upon the works in [14,19], in this work, we integrate the QHMC training into equality-constrained GP regression in a probabilistic manner. Additionally, an adaptive learning algorithm is employed to optimize the selection of constraint locations. This work addresses the computational challenges associated with high-dimensional spaces and large datasets. It offers a flexible alternative to conventional numerical methods, while also being scalable to high-dimensional problems. The effectiveness and precision of the QHMC algorithms are demonstrated across various scenarios, including both low-dimensional and high-dimensional synthetic problems, and one real example for 3-dimensional heat transfer problem under PDE constraints. The main contributions of this work can be outlined in three points: (i) QHMC helps minimize the gap between the posterior mean and the true values, (ii) applying QHMC within a probabilistic framework reduces variance and uncertainty, and (iii) the proposed algorithm serves as a robust, efficient, and flexible approach that can be applied across a wide range of problems.

2. Gaussian Process Under Equality Constraints

2.1. Standard GP Regression Framework

Consider a target function represented by the values $\mathbf{y}={(y^{\left(1\right)},y^{\left(2\right)},...,y^{\left(T\right)})}^N$, where each $y^{\left(i\right)}\in \mathbb{R}$ corresponds to an observation at the locations $\mathbf{X}={\left\{x^{\left(i\right)}\right\}}_{i=1}^N$. Here, each $x^{\left(i\right)}$ is a d-dimensional vector within the domain $\mathcal{D}\in {\mathbb{R}}^d$. In line with the framework from [2], we approximate the target function using a GP, denoted as $Y(.,.):D\times \Omega \to \mathrm{R}$. The GP Y can be written as

$$\begin{array}{c}\hfill Y\left(x\right):=GP[\mu \left(x\right),K(x,x^{\prime })],\end{array}$$

where $\mu (.)$ is the mean function and $K(x,x^{\prime })$ is the covariance function, each specified by

$$\mu \left(x\right)=\mathbb{E}\left[Y\left(x\right)\right],\mathrm{and}K(x,x^{\prime })=\mathbb{E}[Y\left(x\right)-\mu \left(x\right)][Y\left(x^{\prime }\right)-\mu \left(x^{\prime }\right)].$$

A typical choice for the covariance kernel function is the standard squared exponential covariance kernel, which is described as

$$K(x,x^{\prime })={\sigma }^{2}exp\left(-{\displaystyle \frac{\left|\right|x-x^{\prime }{\left|\right|}_2^2}{2l^2}}\right)+{\sigma }_n^2{\delta }_{x,x^{\prime }},$$

where ${\sigma }^2$ represents the signal variance, ${\delta }_{x,x^{\prime }}$ is the Kronecker delta function and l is the length-scale. It is assumed that the observations include an additive independent identically distributed (i.i.d.) Gaussian noise term $ϵ$ and having zero mean and variance ${\sigma }_n^2$. The set of hyperparameters, $\theta =(\sigma ,l,{\sigma }_n)$, is estimated using the training data, typically by minimizing the negative marginal log-likelihood [2,21,22]:

$$\begin{array}{c}\hfill -log\left[p\left(\mathbf{Y}\right|\mathbf{X},\theta )\right]={\displaystyle \frac{1}{2}}{[(\mathbf{y}-\mu )}^{\mathrm{T}}{K}^{-1}(\mathbf{y}-\mu )+log\left|K\right|+Nlog\left(2\pi \right)].\end{array}$$

(1)

In this work, we estimate the set of parameters by performing the QHMC updates.

2.2. Quantum-Inspired Hamiltonian Monte Carlo

is an advanced variation of the traditional HMC algorithm that introduces an additional random mass matrix for the particles, governed by a probability distribution. In standard HMC, the system’s state is described by two main components: the position, represented by the original variables (x), and the Gaussian momentum, modeled through auxiliary variables (q). QHMC expands on this by incorporating principles from quantum mechanics, particularly the energy-time uncertainty relation, to allow the mass of the particle to vary randomly according to a probability distribution. This modification introduces a new variable, the mass (m), alongside the position and momentum variables. The introduction of this third variable offers the key benefit of enabling the QHMC algorithm to explore a broader range of potential landscapes in the state-space, making it more versatile. Consequently, QHMC can more effectively navigate complex distributions, including those that are discontinuous, non-smooth, or have spiky features, areas where traditional HMC may struggle. This enhanced exploration capability gives QHMC an edge in sampling from challenging distributions, as it can adapt to varying local structures in the target distribution [23,24].

The quantum nature of QHMC can be illustrated through the example of a one-dimensional harmonic oscillator, as described in [24]. Consider a ball with a fixed mass m attached to a spring at the origin. Let x represent the displacement of the ball from the origin. The restoring force acting on the ball is given by $F=-kx$, which pulls the ball back toward the origin, causing it to oscillate with a period $T=2\pi \sqrt{{\displaystyle \frac{m}{k}}}$. In standard HMC, the mass m is fixed at a value of 1. QHMC, however, introduces a key enhancement by allowing the mass to vary with time. This enables the ball to have different accelerations and, as a result, explore a wider range of distribution landscapes. Specifically, QHMC can use a shorter time period T, which corresponds to a smaller mass m, to efficiently traverse broad and relatively flat regions of the distribution. On the other hand, when dealing with regions that are spiky or have sharp features, QHMC can adapt by employing a longer time period T, associated with a larger mass m, ensuring that all areas of the landscape are thoroughly explored. This flexibility in adjusting the mass and time scale allows QHMC to better navigate complex and varied distributions compared to standard HMC [24].

The implementation of QHMC is relatively simple: a stochastic process $M\left(t\right)$ is constructed for the mass, and at each time t, the mass $M\left(t\right)$ is sampled from a distribution $P_M\left(M\right)$. The only addition to the standard HMC process is resampling the positive-definite mass matrix at each step. In practice, we assume that $P_M\left(M\right)$ is independent of x and q, and select a mass density function $P_M\left(M\right)$ with mean ${\mu }_m$ and variance ${\sigma }_m^2$ as $logm\sim \mathcal{N}({\mu }_m,{\sigma }_m^2),M=mI,$ where I is the identity matrix. Then, the QHMC framework simulates the following dynamical system:

$$\begin{array}{c}\hfill d\left(\begin{array}{c}x\\ q\end{array}\right)=dt\left(\begin{array}{c}M{\left(t\right)}^{-1}q\\ -\nabla U\left(x\right)\end{array}\right).\end{array}$$

3. Proposed Method

In the QHMC context provided in the previous section, the potential energy function for the QHMC system is given by $U\left(x\right)=-log\left[p\right(\mathbf{Y}|\mathbf{X},\theta )]$, i.e., the negative of marginal log-likelihood. Algorithm 1 summarizes the steps of QHMC sampling, and, here, we consider the location variables ${\left\{x^{\left(i\right)}\right\}}_{i=1}^N$ in GP model as the position variables x in Algorithm 1. The method evolves the QHMC dynamics to update the locations x. In this work, we implement the QHMC method for equality constrained GP regression in a probabilistic manner.

Rather than strictly enforcing all constraints, the methods proposed in [4,25] focuses on minimizing the negative marginal log-likelihood function from Equation 1, while allowing occasional violations of the constraints with a small probability. For example, an equality constrained optimization problem can be modified by introducing the following condition to the problem:

$$P\left[\right(\mathbf{Y}\left(x\right)|x,\theta )\ne c]\le \eta ,\mathrm{for}\mathrm{all}x\in \mathcal{D},$$

(2)

where $0<\eta <<1$. Unlike methods that enforce constraints using a truncated Gaussian assumption [6] or by applying inference techniques such as the Laplace approximation and expectation propagation [26], the proposed approach maintains the Gaussian posterior typical of standard GP regression. This method introduces only a slight modification to the existing cost function. For a model that follows a Gaussian distribution, the constraint can be reformulated in terms of the posterior mean and posterior standard deviation:

$$\begin{array}{c}\hfill y^{*}\left(x\right)+{\varphi }^{-1}\left(\eta \right)s\left(x\right)=0,\mathrm{for}\mathrm{all}x\in \mathcal{D},\end{array}$$

(3)

where $y^{*}\left(x\right)$ stands for the posterior mean, s is the standard deviation and $\varphi$ is the cumulative distribution function of a Gaussian random variable. In this work, we set $\eta$$2.2\%$ for demonstration purposes, as in [4]. This gives ${\varphi }^{-1}\left(\eta \right)=-2$, indicating that the violation is allowed for two standard deviations. Then, the reformulation of the optimization problem is given as

$$\begin{array}{cc}\hfill {\mathrm{argmin}}_{\theta }& -log\left[p\right(\mathbf{Y}|\mathbf{X},\theta )]\mathrm{such}\mathrm{that}\hfill \\ \hfill & c-2s\left(x\right)\le y^{*}\left(x\right)\le c+2s\left(x\right).\hfill \end{array}$$

(4)

In this specific formulation of the optimization problem, there is a functional constraint outlined in Equation 4. It may be impractical or even impossible to satisfy this constraint at every point throughout the entire domain. Consequently, we implement a strategy that enforces Equation 4 only at a chosen set of m constraint points denoted as ${\mathbf{X}}_c={x_c^{\left(i\right)}}_{i=1}^m$. We reformulate the optimization problem:

$$\begin{array}{cc}\hfill {\mathrm{argmin}}_{\theta }& -log\left[p\right(\mathbf{Y}|\mathbf{X},\theta )]\mathrm{such}\mathrm{that}\hfill \\ \hfill & c-2s\left(x_c^{\left(i\right)}\right)\le y^{*}\left(x_c^{\left(i\right)}\right)\le c+2s\left(x_c^{\left(i\right)}\right)\mathrm{for}\mathrm{all}i=1,2,...,m,\hfill \end{array}$$

(5)

where hyperparameters are estimated to enforce bounds. Solving this optimization problem can be quite challenging, leading to adding regularization terms in the approach presented in [4]. Instead of directly solving the optimization problem, this work employs the soft-QHMC method first introduced in [25], which introduces inequality constraints with a high probability (for example, $95\%$) by selecting a specific set of m constraint points within the domain. Non-negativity is then enforced on the posterior GP at these selected points.

The QHMC algorithm is used to minimize the log-likelihood from Equation 1. By utilizing Bayesian estimation techniques [27], we can approximate the posterior distribution based on the log-likelihood function and the prior probability distribution, as outlined below:

$$\begin{array}{c}\hfill p(\mathbf{X},\theta |\mathbf{Y})\propto p(\mathbf{X},\theta ,\mathbf{Y})=p(\theta \left)p\right(\mathbf{X}\left|\theta \right)p\left(\mathbf{Y}\right|\mathbf{X},\theta ).\end{array}$$

(6)

The QHMC training process begins with this Bayesian learning and switch to an MCMC procedure to draw samples generated within the Bayesian framework. A general sampling procedure at step t is given as

$$\begin{array}{cc}\hfill X^{(t+1)}& \sim \pi \left(X\right|\theta )=p\left(X\right|{\theta }^{\left(t\right)},Y),\hfill \\ \hfill {\theta }^{(t+1)}& \sim \pi \left(\theta \right|X)=p\left(\theta \right|X^{(t+1)},Y).\hfill \end{array}$$

(7)

The steps of soft equality-constrained GP regression is outlined in Algorithm 2, where QHMC sampling (provided in Algorithm 1) is used as a GP training method. In this version of equality constrained GP regression, the constraint points are located where the posterior variance is highest.

Algorithm 1 QHMC Training for GP with Equality Constraints

Input: Initial point $x_0$, step size $ϵ$, number of simulation steps L, mass distribution parameters ${\mu }_m$ and ${\sigma }_m$. 1: for$t=1,2,...$do 2:    Resample $M_t\sim P_M\left(M\right)$    Resample $q_t\sim N(0,M_t)$    $(x_0,q_0)=(x^{\left(t\right)},q^{\left(t\right)})$    $q_0\leftarrow q_0-{\displaystyle \frac{ϵ}{2}}\nabla U\left(x_0\right)$ 3:    for $i=1,2,...,L-1$ do 4:       $x_i\leftarrow x_{i-1}+ϵM_t^{-1}{q}_{i-1}$       $q_i\leftarrow q_{i-1}-{\displaystyle \frac{ϵ}{2}}\nabla U\left(x_i\right)$ 5:    end for    $x_L\leftarrow x_{L-1}+ϵM_t^{-1}{q}_{L-1}$    $q_L\leftarrow q_{L-1}-{\displaystyle \frac{ϵ}{2}}\nabla U\left(x_L\right)$    $(\widehat{x},\widehat{q})=(x_L,q_L)$    MH step:$u\sim$ Uniform$[0,1]$;    $\rho =e^{-H(\widehat{x},\widehat{q})+H(x^{\left(t\right)},q^{\left(t\right)})}$; 6:    if $u<min(1,\rho )$ then 7:      $(x^{(t+1)},q^{(t+1)})=(\widehat{x},\widehat{q})$ 8:    else 9:      $(x^{(t+1)},q^{(t+1)}=(x^{\left(t\right)},q^{\left(t\right)})$ 10:    end if 11: end for Output:$\{x^{\left(1\right)},x^{\left(2\right)},...\}$

3.1. Adaptive Learning Mechanism

Instead of randomly selecting m constraint points, the proposed algorithm initiates its process with an empty set of constraints and adaptively determines the locations of the constraint points one by one. This adaptive approach allows for a more sophisticated selection process that can lead to improved performance. Throughout this iterative process, several strategies are employed to effectively add constraints based on the behavior of the model.

The first strategy is known as the constraint-adaptive approach. Specifically, the algorithm continuously examines whether the constraints are being satisfied at various pre-decided locations within the domain. For each candidate location, the function value is calculated, and if a violation of the constraint is detected, a constraint point is then added at that specific location. This real-time evaluation ensures that the algorithm can quickly respond to areas where the constraints are not met, allowing for a more dynamic and responsive optimization process.

The second strategy employed is the variance-adaptive approach. This technique involves calculating the posterior variance within the test set to identify which locations exhibit the highest uncertainty. By focusing on these high-variance areas, the algorithm aims to add constraint points that will help reduce overall prediction variance and thereby increase the stability of the model’s outputs. This approach ensures that the algorithm prioritizes regions that may contribute most significantly to uncertainty, effectively allowing it to enhance the reliability of its predictions.

Additionally, the algorithm employs a hybrid strategy that combines both constraint and variance adaptation. In this combined approach, a threshold value is established for the variance. The algorithm identifies constraint points at locations where the highest prediction variance is observed, aiming to address and stabilize these areas first. Once the variance for these points is successfully reduced to the predetermined threshold value, the algorithm switches back to the constraint-adaptive approach. In this phase, it focuses on locating constraint points in areas where violations are occurring, ensuring that the model remains robust and compliant with the necessary constraints while continually refining its predictions. By integrating these strategies, the algorithm effectively navigates the challenges posed by the optimization problem, leading to more accurate and reliable results.

An example of the variance-adaptive approach can be seen in the workflow of soft equality-constrained GP regression, as outlined in Algorithm 2. In this implementation, QHMC sampling, described in Algorithm 1, serves as the training method for the Gaussian Process (GP). In this specific version of non-negativity enforced GP regression, constraint points are strategically positioned at locations where the posterior variance is at its highest.

Algorithm 2 GP Regression with Soft Constraints

1: Specify m constraint points denoted by ${\mathbf{X}}_c={x_c^{\left(i\right)}}_{i=1}^m$, where corresponding observation $y^{*}{\left(x_c\right)}^{\left(i\right)}.$ 2: for  $i=1,2,...,m$  do 3:    Compute the MSE of $s^2\left(x_c^{\left(i\right)}\right)$ of MLE prediction $y^{*}\left(x_c\right)$ for $x_c\in \mathcal{D}$.    Obtain observation $y^{*}{\left(x_c\right)}^{\left(i\right)}$ at $x_c^{\left(i\right)}$    Locate $x_c^{(i+1)}$ for the maximum of $s^2\left(x_c^{\left(i\right)}\right)$ for $x_c\in \mathcal{D}$. 4: end for Construct the MLE prediction of $y^{*}\left(x\right)$ using QHMC training.

4. Numerical Examples

In this section, we evaluate the performance of the proposed algorithms on various examples including synthetic and real data. Synthetic examples are generated to demonstrate the efficacy of the proposed method on larger datasets and higher dimensions. In these numerical examples, the time and accuracy performances of the algorithms are assessed while varying the dataset size and noise level in the data. We consider the relative $l_2$ error between the posterior mean $y^{*}$ and the true value of the target function $f\left(x\right)$ on a set of test points ${\mathbf{X}}_t={\left\{x_T^{\left(i\right)}\right\}}_{i=1}^{N_t}$ [4]:

$$E=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^{N_t}{[y^{*}\left(x_T^{\left(i\right)}\right)-f\left(x_T^{\left(i\right)}\right)]}^2}{{\sum }_{i=1}^{N_t}f{\left(x_T^{\left(i\right)}\right)}^2}}}.$$

(8)

We solve the log-likelihood minimization problem in Equation 5 in MATLAB. Furthermore, to demonstrate the advantage of QHMC over standard HMC, the proposed method is also implemented using the traditional HMC procedure. The comparison includes the relative error, posterior variance, and execution time for both the QHMC and HMC algorithms, with results presented for each version.

The results for both soft-constrained and hard-constrained versions of QHMC are presented for comparison. The constraint-adaptive, hard-constrained version is labeled as QHMCad, while its soft-constrained counterpart is referred to as QHMCsoftad. Similarly, QHMCvar denotes the variance-focused approach, with its soft-constrained version called QHMCsoftvar. The combination of the two methods, incorporating both hard and soft constraints, is represented by QHMCboth and QHMCsoftboth, respectively.

4.1. Poisson Equation in 2D

Consider the following differential equation [28]

$${\displaystyle \frac{{\partial }^2}{\partial x_1^2}}u\left(x\right)+{\displaystyle \frac{{\partial }^2}{\partial x_2^2}}u\left(x\right)=f\left(x\right).$$

(9)

We generate a synthetic dataset that includes a single-fidelity collection of noise-free observations for the forcing term $f\left(x\right)=-2{\pi }^{2}sin\left(\pi x_1\right)sin\left(\pi x_2\right)$. Additionally, the initial noise-free locations are generated from the exact solution $u\left(x\right)=sin\left(\pi x_1\right)sin\left(\pi x_2\right)$. Beginning with this initial training set, we initiate an active learning iteration loop. In each iteration, a new observation is added to the training set, based on the selected sampling policy.

Table 1 shows the error and time comparison of all algorithms alongside their HMC counterparts, based on an example of a non-noisy data with 50 samples. All the QHMC algorithms outperform their HMC equivalents, with QHMCboth delivering the highest accuracy. Meanwhile, QHMCsoftboth provides a balance, offering the second most accurate results while also being the second fastest. Figure 1 presents more comprehensive accuracy comparison of algorithms with respect to the SNR and datasize. The figure indicates that the proposed algorithms can tolerate the noise in the data, especially when the dataset is larger. Although all the algorithms provide a really small error in this simple 2D example, QHMCboth and QHMCsoftboth are the best ones, while QHMCvar and QHMCsoftvar are outperformed by QHMCad and QHMCsoftad. However, Figure 2 shows that the variance based QHMC approach is faster than both QHMCad and QHMCboth algorithms.

4.2. Poisson Equation in 10D

We consider the 10-dimensional Poisson example in [28] given by

$$f_2\left(x\right)=-8{\pi }^{2}sin\left(2\pi x_1\right)sin\left(2\pi x_3\right),$$

(10)

with the noise term $\{{\mathbf{x}}_{\mathbf{2}},{\mathbf{y}}_{\mathbf{2}}\}$, ${\mathbf{y}}_2=f\left({\mathbf{x}}_2\right)+{ϵ}_2$, where ${ϵ}_2\sim N(0,0.05I)$. The training set contains $n=80$ observations.

Table 2 presents the error and time comparison of QHMC and HMC versions of all approaches with $n=80$ observations, where QHMC approaches provides a significant improvement in accuracy. Additionally, Figure 3 shows the change in relative error of the algorithms for increasing dataset size and noise. Although the error values are fluctuating around $2\%$ for all the methods, QHMCboth results are the most accurate and robust. Based on the comparison in Figure 3, QHMCboth can tolerate noise levels up to $20\%$ with the smallest error, and it can still provide good accuracy (error is around $0.02$). Since this is a 10-dimensional example, the time performances are more important than the previous example, which are shown in Figure 4. Here, we can clearly observe the time advantage of soft-constrained approaches, especially for QHMCsoftvar and QHMCsoftboth. Considering the accuracy and time efficiency together, QHMCsoftboth provides the most efficient and robust method for especially larger datasets.

4.3. Heat Transfer in a Hallow Sphere

This 3-dimensional example considers a heat transfer problem in a hallow sphere. Let $B_r\left(0\right)$ represent a ball centered at 0 with radius r. Defining the hallow sphere as $D=B_4\left(0\right)-B_2\left(0\right)$, the equations are given as [29]

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\displaystyle \frac{\partial u(\mathbf{x},t)}{\partial t}}-\nabla ·(\kappa \nabla u(\mathbf{x},t))=0,\hfill & \mathbf{x}\in D,\hfill \\ \kappa {\displaystyle \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{n}}}={\theta }^2{(\pi -\theta )}^2{\varphi }^2{(\pi -\varphi )}^2,\hfill & {\mathrm{if}\parallel \mathbf{x}|}_2=4\mathrm{and}\varphi \ge 0,\hfill \\ \kappa {\displaystyle \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{n}}}=0,\hfill & {\mathrm{if}\parallel \mathbf{x}\parallel }_2=4\mathrm{and}\varphi <0,\hfill \\ u(\mathbf{x},t)=0,\hfill & {\mathrm{if}\parallel \mathbf{x}\parallel }_2=2.\hfill \end{array}\right.\end{array}$$

(11)

In this context, $\mathbf{n}$ denotes the normal vector pointing outward, while $\theta$ and $\varphi$ represent the azimuthal and elevation angles, respectively, of points within the sphere. We determine the precise heat conductivity using $\kappa =1.0+exp\left(0.05u\right)$. The quadratic elements with 12,854 degrees of freedom are employed, and we set $y\left(\mathbf{x}\right)=u(\mathbf{x},10)$ to solve the partial differential equations (PDE).

We assume that the accurate observation data are available at six locations. We start with those initial locations, and then add six new constraint locations one by one, using the adaptive learning approach. Table 3 shows the final error and time results for all types of algorithms, including HMC versions. For each variant of QHMC and HMC, employing QHMC sampling in a particular version speeds up the process and enhances accuracy. Overall, the comparison reveals that among all the QHMC and HMC variants, QHMCboth delivers the highest accuracy, while QHMCsoftboth is the fastest and ranks second in terms of accuracy. The comparison among QHMC versions are shown in detail in Figure 5, where step-by-step decrease shows the effect of adaptive approaches and the superiority of QHMCboth and QHMCsoftboth algorithms. In this real data example, we observe the same behaviour as in the Poisson examples: QHMCbith is the most accurate, while QHMCsoftboth is the fastest. On the other hand, QHMCsoftboth ranks second for accuracy, therefore it should be the most efficient, since the error values are close to each other.

5. Discussion

Building up on the work in [14], this work evaluate both the advantages and limitations of applying a soft-constrained approach to physics-informed GP regression, focusing on equality constraints. To validate the superiority of the approach, comparisons are also made between modified versions of the proposed algorithm. The key findings and the potential underlying reasons are summarized as follows:

• To demonstrate robustness and effectiveness of the proposed methods, synthetic examples are designed with a focus on two main factors: the size of the dataset and the signal-to-noise ratio (SNR). The QHMC-based algorithms were tested across varying levels of SNR, up to $20\%$. The findings, illustrated in Figure 1 and Figure 3 demonstrate a consistent ability of the method —both in its soft-constrained and hard-constrained forms— to handle noise effectively, especially when the noise remains below $15\%$. Furthermore, the algorithms demonstrated improved performance with larger datasets, indicating that increased data size helps mitigate the impact of noise. These observations from the synthetic examples collectively confirm the robustness of the proposed method under varying conditions.
• The numerical results for the synthetic examples also provide execution times as the SNR and dataset size increase in each case, aiming to highlight the efficiency of the proposed algorithm. Figure 2, Figure 4 demonstrate the time-saving benefits of the algorithms, with the soft-constrained versions showing particularly notable advantages.
• To test the ability of the algorithms to handle higher-dimensional problems while maintaining robustness and efficiency, the evaluations focused on different dimensional settings: 2-dimensions and 10-dimensions. The results confirmed that the proposed methods could consistently deliver accurate results, even as dimensionality increased, and did so in relatively short execution times. This demonstrates the scalability of the algorithms across a range of problem complexities.
• To demonstrate the potential of the proposed method to generalize across different types of problems, we selected a real example: a 3D heat transfer problem that requires solving partial differential equations. Unlike the synthetic examples, in this part we use a fixed dataset size and contain no injected Gaussian noise. A thorough comparison of all methods, including the step-by-step change in the relative error is presented, validating the success of all versions.
• It is also worth noting that while the numerical results demonstrate the robustness and efficiency of the current QHMC algorithm, further investigation is needed regarding the effect of the probability of constraint violation. The experiments were conducted with a relatively low constraint release probability (approximately $5\%$), and accuracy was maintained under these conditions. However, increasing the allowance for violations could introduce limitations.

6. Conclusions

This work combines the accurate training capabilities of QHMC with the efficiency of a probabilistic approach to create a new soft-constrained QHMC algorithm. The algorithm is specifically designed to enforce equality constraints on GP regression. The proposed approach ensures more accurate predictions while also boosting efficiency by delivering successful results in a relatively short time. To further optimize the performance of QHMC algorithms across different problem settings, this work introduces several modified versions of QHMC that incorporate adaptive learning techniques. These adaptations allow the algorithm to dynamically adjust its behavior, improving convergence speed and robustness under varying conditions. By implementing these adaptive versions, the approach provides greater flexibility, enabling users to choose the most appropriate algorithm based on the specific characteristics and priorities of the problem at hand, such as balancing accuracy, computational time, or the handling of noisy data. This tailored approach broadens the potential applications of the algorithms and ensures more consistent results across a wide range of constrained optimization tasks.

The methods are implemented to solve a variety of problems, demonstrating their versatility and effectiveness. In each experiment, we started by demonstrating the benefits of QHMC sampling compared to traditional HMC sampling, with QHMC consistently outperforming HMC. Across all tested cases, QHMC sampling delivered approximately a $20\%$ reduction in computation time and a $15\%$ improvement in accuracy. Building on these initial findings, we conducted a more extensive evaluation of the performance of the algorithms across different scenarios. The examples included higher-dimensional problems, illustrating the capability of the algorithms in tackling complex tasks. Additionally, we tested the methods on a real-world application where we solve a PDE. These evaluations highlight the adaptability of the proposed approach and confirm its effectiveness in a broad range of practical and theoretical settings. Overall, this work has shown that the soft-constrained QHMC method is a robust, efficient, and adaptable approach, capable of handling high-dimensional problems and large datasets. The numerical results indicate that the soft-constrained QHMC has strong potential for generalization across a wide range of applications, including those with varying physical properties.