∎

11institutetext: André L. Marchildon 22institutetext: University of Toronto
Institute for Aerospace Studies, Toronto, ON, Canada
22email<EMAIL_ADDRESS>33institutetext: David W. Zingg
44institutetext: University of Toronto Institute for Aerospace Studies,
Toronto, ON, Canada

# A Solution to the Ill-Conditioning of Gradient-Enhanced Covariance Matrices
for Gaussian Processes

André L. Marchildon  David W. Zingg

(Received: date / Accepted: date)

###### Abstract

Gaussian processes provide probabilistic surrogates for various applications
including classification, uncertainty quantification, and optimization. Using
a gradient-enhanced covariance matrix can be beneficial since it provides a
more accurate surrogate relative to its gradient-free counterpart. An acute
problem for Gaussian processes, particularly those that use gradients, is the
ill-conditioning of their covariance matrices. Several methods have been
developed to address this problem for gradient-enhanced Gaussian processes but
they have various drawbacks such as limiting the data that can be used,
imposing a minimum distance between evaluation points in the parameter space,
or constraining the hyperparameters. In this paper a new method is presented
that applies a diagonal preconditioner to the covariance matrix along with a
modest nugget to ensure that the condition number of the covariance matrix is
bounded, while avoiding the drawbacks listed above. Optimization results for a
gradient-enhanced Bayesian optimizer with the Gaussian kernel are compared
with the use of the new method, a baseline method that constrains the
hyperparameters, and a rescaling method that increases the distance between
evaluation points. The Bayesian optimizer with the new method converges the
optimality, i.e. the $\ell_{2}$ norm of the gradient, an additional 5 to 9
orders of magnitude relative to when the baseline method is used and it does
so in fewer iterations than with the rescaling method. The new method is
available in the open source python library GpGradPy, which can be found at
https://github.com/marchildon/gpgradpy/tree/paper_precon. All of the figures
in this paper can be reproduced with this library.

###### Keywords:

Gaussian process Covariance matrix Condition number Bayesian optimization

###### MSC:

15A12, 60G15, 65K99

##  1 Introduction

In diverse fields and for various applications, such as uncertainty
quantification, classification, and optimization, an expensive function of
interest must be repeatedly evaluated zingg_comparative_2008 ;
shahriari_taking_2016 ; schulz_tutorial_2018 . To minimize the computational
cost it is desirable to minimize the number of expensive function evaluations.
One way to achieve this is by constructing a surrogate that approximates the
function of interest and is inexpensive to evaluate. Various methods to
construct surrogates can be utilized, including using a Gaussian process (GP).
This method provides a nonparametric probabilistic surrogate. The
nonparametric component of the GP indicates that it does not depend on a
parametric functional form, unlike a polynomial surrogate where the order of
the basis function must be selected a priori. As for the probabilistic
component of the GP, this enables the surrogate to provide an estimate for the
function of interest and to quantify the uncertainty in its estimate
rasmussen_gaussian_2006 . A Gaussian process requires a mean and a covariance
function ababou_condition_1994 . A constant is often used for the former and
its value is set by maximizing the marginal log-likelihood toal_kriging_2008 ;
toal_adjoint_2009 ; toal_development_2011 ; ollar_gradient_2017 . For the
covariance function, it is popular to use kernels, of which many are available
davis_six_1997 ; rasmussen_gaussian_2006 . The most popular kernel is the
Gaussian kernel, which is also known as the squared exponential kernel
rasmussen_gaussian_2006 ; shahriari_taking_2016 ; wu_exploiting_2018 . The
desirable properties of this kernel include its hyperparameters that can be
tuned, its simplicity, and its smoothness. This final property enables the
surrogate to be constructed using gradient evaluations, which makes the
surrogate more accurate dalbey_efficient_2013 ; eriksson_scaling_2018 ;
wu_exploiting_2018 .

Gradient-enhanced GPs use both the value and gradient of the function of
interest to construct the probabilistic surrogate. By using gradients with the
GP, a more accurate surrogate is constructed that matches both the value and
gradient of the function of interest where it has been evaluated in the
parameter space osborne_gaussian_2009 ; ulaganathan_performance_2016 ;
wu_bayesian_2017 . This is particularly useful in high-dimensional parameter
spaces since a single gradient evaluation provides much more information than
a single function evaluation. The gradient-enhanced covariance matrix can be
constructed either with the direct method or the indirect method
zimmermann_maximum_2013 . The former modifies the structure of the gradient-
free covariance matrix while the latter does not. The direct method is much
more common han_improving_2013 ; dalbey_efficient_2013 ; wu_exploiting_2018 ;
laurent_overview_2019  and is used in this paper. A drawback of using
gradient-enhanced GPs is that the covariance matrix is larger than its
gradient-free counterpart and is thus more expensive to invert. Various
strategies have been developed to mitigate this additional cost by using
random Fourier features hung_random_2021 , or by exploiting the structure of
the gradient-enhanced covariance matrix de_roos_high-dimensional_2021 .

A ubiquitous problem in the use of GPs is the ill-conditioning of their
covariance matrices ababou_condition_1994 ; kostinski_condition_2000 ;
zimmermann_condition_2015 . This problem is present with the use of many
kernels, including the Gaussian kernel. Various factors have been identified
that exacerbate the ill-conditioning, such as having the data points too close
together davis_six_1997 . The ill-conditioning of the covariance matrix is
problematic since it can cause the Cholesky factorization to fail
higham_cholesky_2009 , and it also increases the numerical error. Adding a
small positive nugget to the diagonal of the gradient-free covariance matrix
is sufficient to ensure that the condition number of the matrix is below a
user-set threshold mohammadi_analytic_2017 .

The ill-conditioning of the gradient-enhanced covariance matrix is even more
acute than the gradient-free case, and the addition of a nugget is
insufficient on its own to alleviate this problem he_instability_2018 ;
dalbey_efficient_2013 . Various approaches have been attempted to mitigate the
ill-conditioning problem, such as removing certain data points until the
condition number is sufficiently low march_gradient-based_2011 ;
dalbey_efficient_2013 , or imposing a minimum distance constraint between data
points in the parameter space osborne_gaussian_2009 . Both methods have
significant drawbacks since they restrict the data available to construct the
surrogate. Furthermore, neither method guarantees that the condition number of
the covariance matrix remains below a user-set threshold as the
hyperparameters are optimized. There is one recent method that does ensure
that the condition number of the gradient-enhanced covariance matrix remains
below a user-set threshold when the Gaussian kernel is used marchildon_non-
intrusive_2023 . This method uses non-isotropic rescaling of the data in order
to have a set minimum distance between the data points. While data points
cannot be collocated, they can get arbitrarily close in the parameter space,
and the method allows all of the data points to be kept in the construction of
the gradient-enhanced covariance matrix. However, the drawback of this method
is that, in some cases, the rescaling needs to be done iteratively, which
requires the hyperparameters to be optimized again. This adds additional
complexity and computational cost.

The new method presented in this paper shares the same benefits as the
rescaling method from marchildon_non-intrusive_2023 , i.e. all of the data
points can be used, there is no minimum distance constraint between the data
points in the parameter space, and the condition number of the gradient-
enhanced covariance matrix is bounded. The new method also has two additional
benefits: it only requires a single optimization of the hyperparameters, i.e.
it is not iterative, and there is no need for a constraint on the condition
number for the optimization of the hyperparameters. This simplifies the
implementation of the new method and reduces its computational cost.

The new and rescaling methods are available in the open source python library
GpGradPy, which can be accessed at
https://github.com/marchildon/gpgradpy/tree/paper_precon. This library
contains the Gaussian, Matérn $\frac{5}{2}$, and rational quadratic kernels.

The notation used in this paper is presented in Section 2. In Section 3 the GP
is presented along with the Gaussian kernel and the covariance matrix. A
modified covariance matrix is derived in Section 4. In Section 5 it is
demonstrated how the condition number of the modified covariance matrix can be
bounded with the use of a nugget. Details on the implementation of the new
covariance matrix can be found in Section 6 and optimization results are
provided in Section 7. Finally, the conclusions of the paper are presented in
Section 8.

##  2 Notation

Sans-serif capital letters are used for matrices. For example, $\mathsf{I}$ is
the identity matrix, and $\mathsf{X}$ is an $n_{x}\times d$ matrix that holds
the location of $n_{x}$ evaluation points in a $d$ dimensional parameter
space. Vectors are denoted in lowercase bold font. For instance,
$\boldsymbol{x}$ and $\boldsymbol{y}$ are vectors of length $d$ denoting
arbitrary points in the parameter space. The $i$-th row of $\mathsf{X}$ is
denoted as $\boldsymbol{x}_{i:}$ and its $j$-th column is indicated as
$\boldsymbol{x}_{:j}$. Finally, scalars are denoted in lowercase letters such
as $x_{ij}$, which is the entry at the $i$-th row and $j$-th column of
$\mathsf{X}$. The symbols $\textbf{0}_{d}$ and $\textbf{1}_{d}$ are vectors of
length $d$ with all of their entries equal to zero and one, respectively.
Variations of these symbols such as $\bar{\textbf{1}}$ or $\hat{\textbf{1}}$
are used to indicate a matrix of ones or a vector where some of its entries
are zero, respectively. These will be clarified when they appear in the paper.

##  3 Gaussian process

###  3.1 Gradient-free covariance matrix

To fully define a GP we require a mean function and a covariance function. The
mean function is selected here to be the constant $\beta$, which is a
hyperparameter that is selected by maximizing the marginal log-likelihood
function that will be presented in Section 3.3. For the covariance function we
use the popular Gaussian kernel

$\displaystyle k(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\gamma})$
$\displaystyle=e^{-\frac{1}{2}\sum_{i=1}^{d}\gamma_{i}^{2}(x_{i}-y_{i})^{2}},$
(1)

where $\gamma_{i}>0\,\forall\,i\in\\{1,\ldots,d\\}$ are hyperparameters. The
Gaussian kernel is typically presented with
$\boldsymbol{\theta}=\boldsymbol{\gamma}^{2}/2$ as its hyperparameters but it
is simpler in the later derivations to use $\boldsymbol{\gamma}$ instead. The
Gaussian kernel is a stationary kernel since it depends only on
$\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{y}$, i.e. the relative location of
$\boldsymbol{y}$ to $\boldsymbol{x}$. The gradient-free Gaussian kernel matrix
is

$\mathsf{K}=\mathsf{K}(\mathsf{X};\boldsymbol{\gamma})=\begin{bmatrix}1&k(\boldsymbol{x}_{1:},\boldsymbol{x}_{2:},;\boldsymbol{\gamma})&\ldots&k(\boldsymbol{x}_{1:},\boldsymbol{x}_{n_{x}:};\boldsymbol{\gamma})\\\
k(\boldsymbol{x}_{2:},\boldsymbol{x}_{1:};\boldsymbol{\gamma})&1&\ldots&k(\boldsymbol{x}_{2:},\boldsymbol{x}_{n_{x}:};\boldsymbol{\gamma})\\\
\vdots&\vdots&\ddots&\vdots\\\
k(\boldsymbol{x}_{n_{x}:},\boldsymbol{x}_{1:};\boldsymbol{\gamma})&k(\boldsymbol{x}_{n_{x}:},\boldsymbol{x}_{2:};\boldsymbol{\gamma})&\ldots&1,\end{bmatrix},$
(2)

where its diagonal entries are all unity. In general, the $i$-th diagonal
entry of $\mathsf{K}$ is
$k(\boldsymbol{x}_{i:},\boldsymbol{x}_{i:};\boldsymbol{\gamma})$. The
gradient-free Gaussian kernel matrix $\mathsf{K}$ is also a correlation matrix
since it satisfies all of the properties of the following definition.

###### Definition 1

A correlation matrix must satisfy all of the following conditions:

  1. 1.

All of the entries in the square matrix are real and between $-1$ and $1$

  2. 2.

The diagonal entries of the matrix are all unity

  3. 3.

The matrix is positive semidefinite

The noise-free regularized gradient-free covariance matrix is given by

$\Sigma(\mathsf{X};\boldsymbol{\gamma},\eta_{\mathsf{K}})=\sigma_{\mathsf{K}}^{2}\left(\mathsf{K}(\mathsf{X};\boldsymbol{\gamma})+\eta_{\mathsf{K}}\mathsf{I}\right),$
(3)

where the hyperparameter $\sigma_{\mathsf{K}}^{2}$ is the variance of the
stationary residual error, and the nugget $\eta_{\mathsf{K}}$ is used to have
$\kappa(\Sigma)\leq\kappa_{\max}$, where $\kappa_{\max}>1$ is the maximum
allowed condition number and is set by the user. The nugget is discussed
further in Section 5.1.

###  3.2 Gradient-enhanced covariance matrix

Constructing the gradient-enhanced kernel matrix requires the derivatives of
the kernel with respect to its inputs:

$\displaystyle\frac{\partial k(\boldsymbol{x},\boldsymbol{y})}{\partial
x_{i}}$
$\displaystyle=-\gamma_{i}^{2}\left(x_{i}-y_{i}\right)k(\boldsymbol{x},\boldsymbol{y})$
(4) $\displaystyle\frac{\partial k(\boldsymbol{x},\boldsymbol{y})}{\partial
y_{j}}$
$\displaystyle=\gamma_{j}^{2}\left(x_{j}-y_{j}\right)k(\boldsymbol{x},\boldsymbol{y})$
(5) $\displaystyle\frac{\partial^{2}k(\boldsymbol{x},\boldsymbol{y})}{\partial
x_{i}\partial y_{j}}$
$\displaystyle=\left(\delta_{ij}\gamma_{i}^{2}-\gamma_{i}^{2}\gamma_{j}^{2}\left(x_{i}-y_{i}\right)\left(x_{j}-y_{j}\right)\right)k(\boldsymbol{x},\boldsymbol{y}),$
(6)

where $\delta_{ij}$ is the Kronecker delta. The gradient-enhanced kernel
matrix is given by

$\displaystyle\mathsf{K}_{\nabla}(\mathsf{X};\boldsymbol{\gamma})$
$\displaystyle=\begin{bmatrix}\mathsf{K}&\frac{\partial\mathsf{K}}{\partial
y_{1}}&\ldots&\frac{\partial\mathsf{K}}{\partial y_{d}}\\\
\frac{\partial\mathsf{K}}{\partial
x_{1}}&\frac{\partial^{2}\mathsf{K}}{\partial x_{1}\partial
y_{1}}&\ldots&\frac{\partial^{2}\mathsf{K}}{\partial x_{1}\partial y_{d}}\\\
\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial\mathsf{K}}{\partial
x_{d}}&\frac{\partial^{2}\mathsf{K}}{\partial x_{d}\partial
y_{1}}&\ldots&\frac{\partial^{2}\mathsf{K}}{\partial x_{d}\partial
y_{d}}\end{bmatrix}$ (7)
$\displaystyle=\begin{bmatrix}\mathsf{K}&\gamma_{1}^{2}\mathsf{R}_{1}\odot\mathsf{K}&\ldots&\gamma_{d}^{2}\mathsf{R}_{d}\odot\mathsf{K}\\\
-\gamma_{1}^{2}\mathsf{R}_{1}\odot\mathsf{K}&\left(\gamma_{1}^{2}\bar{\textbf{1}}-\gamma_{1}^{4}\mathsf{R}_{1}^{(2)}\right)\odot\mathsf{K}&\ldots&-\gamma_{1}^{2}\gamma_{d}^{2}\mathsf{R}_{1}\odot\mathsf{R}_{d}\odot\mathsf{K}\\\
\vdots&\vdots&\ddots&\vdots\\\
-\gamma_{d}^{2}\mathsf{R}_{d}\odot\mathsf{K}&-\gamma_{1}^{2}\gamma_{d}^{2}\mathsf{R}_{1}\odot\mathsf{R}_{d}\odot\mathsf{K}&\ldots&\left(\gamma_{d}^{2}\bar{\textbf{1}}-\gamma_{d}^{4}\mathsf{R}_{d}^{(2)}\right)\odot\mathsf{K}\\\
\end{bmatrix},$ (8)

where $\bar{\textbf{1}}$ is a matrix of ones of size $n_{x}\times n_{x}$, the
operator $\odot$ is the Hadamard product for elementwise multiplication, and
$\mathsf{R}_{i}$ is a skew-symmetric matrix given by

$\displaystyle\mathsf{R}_{i}(\mathsf{X})$
$\displaystyle=\boldsymbol{x}_{:i}\textbf{1}_{n_{x}}^{\top}-\textbf{1}_{n_{x}}\boldsymbol{x}_{:i}^{\top}$
$\displaystyle=\begin{bmatrix}0&x_{1i}-x_{2i}&\ldots&x_{1i}-x_{n_{x}i}\\\
x_{2i}-x_{1i}&0&\ldots&x_{2i}-x_{n_{x}i}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{n_{x}i}-x_{1i}&x_{n_{x}i}-x_{2i}&\ldots&0\end{bmatrix}.$ (9)

Unlike $\mathsf{K}$, $\mathsf{K}_{\nabla}$ is not a correlation matrix since
it does not satisfy the first and second conditions in Definition 1. This is
clear from checking the diagonal of $\mathsf{K}_{\nabla}$:

$\operatorname{diag}(\mathsf{K}_{\nabla})=[\underbrace{1,\ldots,1}_{n_{x}},\underbrace{\gamma_{1}^{2},\ldots\gamma_{1}^{2}}_{n_{x}},\ldots,\underbrace{\gamma_{d}^{2},\ldots\gamma_{d}^{2}}_{n_{x}}].$
(10)

Definition 1 would only be satisfied for $\mathsf{K}_{\nabla}$ if
$\gamma_{1}=\ldots=\gamma_{d}=1$. However, the hyperparameters
$\boldsymbol{\gamma}$ are set by maximizing the marginal log-likelihood, which
is introduced in the following subsection. The noise-free regularized
gradient-enhanced covariance matrix is

$\Sigma_{\nabla}(\boldsymbol{\gamma};\eta_{\mathsf{K}_{\nabla}})=\sigma_{\mathsf{K}}^{2}\left(\mathsf{K}_{\nabla}(\boldsymbol{\gamma})+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right),$
(11)

where the nugget $\eta_{\mathsf{K}_{\nabla}}$ is a regularization term that is
used to ensure that $\kappa(\Sigma_{\nabla})\leq\kappa_{\max}$, as detailed in
Section 5.

###  3.3 Evaluating the Gaussian process

The mean and variance of the gradient-enhanced Gaussian process are evaluated
with wu_bayesian_2017

$\displaystyle\mu_{f}(\boldsymbol{x})$
$\displaystyle=\beta+\boldsymbol{k}_{\nabla}^{\top}(\boldsymbol{x})\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)$
(12) $\displaystyle\sigma_{f}^{2}(\boldsymbol{x})$
$\displaystyle=\sigma_{\mathsf{K}}^{2}\left(k(\boldsymbol{x},\boldsymbol{x})-\boldsymbol{k}_{\nabla}^{\top}(\boldsymbol{x})\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\boldsymbol{k}_{\nabla}(\boldsymbol{x})\right),$
(13)

where
$\hat{\textbf{1}}=[\textbf{1}_{n_{x}}^{\top},\textbf{0}_{n_{x}d}^{\top}]^{\top}$,
and

$\displaystyle\boldsymbol{k}_{\nabla}(\boldsymbol{x};\mathsf{X})=\begin{bmatrix}\boldsymbol{k}(\mathsf{X},\boldsymbol{x})\\\
\frac{\partial\boldsymbol{k}(\mathsf{X},\boldsymbol{x})}{\partial x_{1}}\\\
\vdots\\\ \frac{\partial\boldsymbol{k}(\mathsf{X},\boldsymbol{x})}{\partial
x_{d}}\end{bmatrix},\quad\boldsymbol{f}_{\nabla}(\mathsf{X})=\begin{bmatrix}\boldsymbol{f}(\mathsf{X})\\\
\frac{\partial\boldsymbol{f}(\mathsf{X})}{\partial x_{1}}\\\ \vdots\\\
\frac{\partial\boldsymbol{f}(\mathsf{X})}{\partial x_{d}}\end{bmatrix},$ (14)

where $\boldsymbol{f}(\mathsf{X})$ is the function of interest evaluated at
all of the rows in $\mathsf{X}$. In this paper the gradient of the function of
interest is calculated analytically, but it could also be calculated with
algorithmic differentiation or approximated with finite differences.

Prior to evaluating the GP, its hyperparameters must first be selected, which
is commonly done by maximizing the marginal likelihood toal_kriging_2008 ;
toal_adjoint_2009 ; toal_development_2011 ; ollar_gradient_2017

$\displaystyle
L(\boldsymbol{\gamma},\beta,\sigma_{\mathsf{K}}^{2};\mathsf{X},\boldsymbol{f}_{\nabla},\eta_{\mathsf{K}_{\nabla}})$
$\displaystyle=\frac{e^{-\frac{\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)^{\top}\Sigma_{\nabla}^{-1}\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)}{2}}}{\left(2\pi\right)^{\frac{n_{x}(d+1)}{2}}\sqrt{\det\left(\Sigma_{\nabla}\right)}}$
$\displaystyle=\frac{e^{-\frac{\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)^{\top}\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)}{2\sigma_{\mathsf{K}}^{2}}}}{\left(2\pi\sigma_{\mathsf{K}}^{2}\right)^{\frac{n_{x}(d+1)}{2}}\sqrt{\det\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)}}.$

For the noise-free case being considered it is straightforward to get closed-
form solutions for $\beta$ and $\sigma_{\mathsf{K}}^{2}$ that maximize
$\ln(L)$ toal_kriging_2008 :

$\displaystyle\beta(\boldsymbol{\gamma};\mathsf{X},\boldsymbol{f}_{\nabla},\eta_{\mathsf{K}_{\nabla}})$
$\displaystyle=\frac{\hat{\textbf{1}}^{\top}\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\boldsymbol{f}_{\nabla}}{\hat{\textbf{1}}^{\top}\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\hat{\textbf{1}}}$
(15)
$\displaystyle\sigma_{\mathsf{K}}^{2}(\boldsymbol{\gamma};\mathsf{X},\boldsymbol{f}_{\nabla},\eta_{\mathsf{K}_{\nabla}},\beta)$
$\displaystyle=\frac{\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)^{\top}\left(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}\right)^{-1}\left(\boldsymbol{f}_{\nabla}-\beta\hat{\textbf{1}}\right)}{n_{x}(d+1)}.$
(16)

Substituting these solutions for $\beta$ and $\sigma_{\mathsf{K}}^{2}$ into
$\ln(L)$ and dropping the constant terms gives

$\ln(L(\boldsymbol{\gamma};\mathsf{X},\eta_{\mathsf{K}_{\nabla}},\sigma_{\mathsf{K}}))=-\frac{n_{x}(d+1)\ln(\sigma_{\mathsf{K}}^{2})+\ln(\det(\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}))}{2}.$
(17)

The hyperparameters $\boldsymbol{\gamma}$ are set by maximizing Eq. (17)
numerically with the bound $\boldsymbol{\gamma}>0$.

(a) Gradient-free GP

(b) Gradient-enhanced GP

Figure 1: GPs with and without gradients that are approximating the function
from Eq. (18) with $\beta=-0.62$, $\sigma_{\mathsf{K}}^{2}=1.07$, and
$\gamma=1.7$.

To highlight the benefit of using gradients we compare a gradient-free and a
gradient-enhanced GP that approximate the following function:

$f(x)=\sin(x)+\sin\left(\frac{10x}{3}\right),$ (18)

which was evaluated at $\boldsymbol{x}=[3.5,4.5,5.5,6.5]^{\top}$. Maximizing
Eq. (17) numerically for the gradient-enhanced GP provides the following
hyperparameters: $\beta=-0.62$, $\sigma_{\mathsf{K}}^{2}=1.07$, and
$\gamma=1.7$. The gradient-free and gradient-enhanced GPs can be seen in Figs.
1(a) and 1(b), respectively. For Fig. 1(b) the black line is Eq. (18), the
solid green line is the mean of the surrogate from Eq. (12), and the light
green area represents $\pm 2\sigma_{f}(x)$, where $\sigma_{f}$ comes from Eq.
(13). For Fig. 1(a) the mean and standard deviation of the gradient-free GP
are calculated with equations analogous to Eqs. (12) and (13) that omit the
gradient evaluations and use the gradient-free kernel matrix.

It is clear from Fig. 1 that the use of gradients to construct the GP
significantly improves its accuracy and also reduces its uncertainty, i.e.
$\sigma_{f}$. The benefit of using gradients to construct a GP is even greater
for higher-dimensional parameter spaces since the gradient provides more
information as the number of parameters increases. However, a significant
problem for gradient-enhanced GPs is that their covariance matrices becomes
extremely ill-conditioned, which is addressed in the following section.

##  4 Modified gradient-enhanced covariance matrix

A modified gradient-enhanced kernel matrix $\tilde{\mathsf{K}}_{\nabla}$ is
derived that can be regularized with a modest nugget such that its condition
number is bounded below the user-set threshold $\kappa_{\max}$. From Eq. (11)
the noise-free and unregularized gradient-enhanced covariance matrix for the
vector $\boldsymbol{f}_{\nabla}$ from Eq. (14) is
$\Sigma_{\nabla}=\sigma_{\mathsf{K}}^{2}\mathsf{K}_{\nabla}$. A correlation
matrix can be formed from a covariance matrix by normalizing by the standard
deviation of the random variables, i.e. $\boldsymbol{f}_{\nabla}$, which are
the square root of the values along the diagonal of $\Sigma_{\nabla}$. Our
gradient-enhanced correlation matrix $\tilde{\mathsf{K}}_{\nabla}$ is thus
given by

$\displaystyle\tilde{\mathsf{K}}_{\nabla}$
$\displaystyle=\left(\frac{1}{\sigma_{\mathsf{K}}}\mathsf{P}^{-1}\right)^{-1}\Sigma_{\nabla}\left(\frac{1}{\sigma_{\mathsf{K}}}\mathsf{P}^{-1}\right)^{-1}$
$\displaystyle=\mathsf{P}^{-1}\mathsf{K}_{\nabla}\mathsf{P}^{-1},$ (19)

where

$\displaystyle\mathsf{P}$
$\displaystyle=\operatorname{diag}\left(\sqrt{\operatorname{diag}\left(\mathsf{K}_{\nabla}\right)}\right)$
(20)
$\displaystyle=\operatorname{diag}(\textbf{1}_{n_{x}},\underbrace{\gamma_{1},\ldots,\gamma_{1}}_{n_{x}},\ldots,\underbrace{\gamma_{d},\ldots,\gamma_{d}}_{n_{x}}).$
(21)

From Eq. (8), $\mathsf{K}_{\nabla}$ is formed with the block matrices composed
of $\mathsf{K}$, as well as its first and second derivatives with respect to
the entries of the dimensional vectors $\boldsymbol{x}$ and $\boldsymbol{y}$
from Eq. (1). Similarly $\tilde{\mathsf{K}}_{\nabla}$ is formed with block
matrices composed of $\mathsf{K}$ as well as its first and second derivatives
with respect to the dimensionless variables $\tilde{x}_{i}=\gamma_{i}x_{i}$
and $\tilde{y}_{i}=\gamma_{i}y_{i}$ for $i\in\\{1,\ldots,d\\}$. Consider the
following chain rule

$\displaystyle\frac{\partial
k(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\gamma})}{\partial\tilde{x}_{i}}$
$\displaystyle=\frac{\partial
k(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\gamma})}{\partial
x_{i}}\frac{\partial x_{i}}{\partial\tilde{x}_{i}}$
$\displaystyle=\frac{\partial
k(\boldsymbol{x},\boldsymbol{y};\boldsymbol{\gamma})}{\partial
x_{i}}\frac{1}{\gamma_{i}},$ (22)

where the term $\frac{\partial x_{i}}{\partial\tilde{x}_{i}}=\gamma_{i}^{-1}$
is contained along the diagonal of $\mathsf{P}^{-1}$. The modified gradient-
enhanced kernel matrix can thus be calculated with

$\displaystyle\tilde{\mathsf{K}}_{\nabla}(\mathsf{X};\boldsymbol{\gamma})$
$\displaystyle=\mathsf{P}^{-1}\mathsf{K}_{\nabla}\mathsf{P}^{-1}$
$\displaystyle=\begin{bmatrix}\mathsf{K}&\frac{\partial\mathsf{K}}{\partial\tilde{y}_{1}}&\ldots&\frac{\partial\mathsf{K}}{\partial\tilde{y}_{d}}\\\
\frac{\partial\mathsf{K}}{\partial\tilde{x}_{1}}&\frac{\partial^{2}\mathsf{K}}{\partial\tilde{x}_{1}\partial\tilde{y}_{1}}&\ldots&\frac{\partial^{2}\mathsf{K}}{\partial\tilde{x}_{1}\partial\tilde{y}_{d}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{\partial\mathsf{K}}{\partial\tilde{x}_{d}}&\frac{\partial^{2}\mathsf{K}}{\partial\tilde{x}_{d}\partial\tilde{y}_{1}}&\ldots&\frac{\partial^{2}\mathsf{K}}{\partial\tilde{x}_{d}\partial\tilde{y}_{d}}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\mathsf{K}&\tilde{\mathsf{R}}_{1}\odot\mathsf{K}&\ldots&\tilde{\mathsf{R}}_{d}\odot\mathsf{K}\\\
-\tilde{\mathsf{R}}_{1}\odot\mathsf{K}&\left(\bar{\textbf{1}}-\tilde{\mathsf{R}}_{1}\odot\tilde{\mathsf{R}}_{1}\right)\odot\mathsf{K}&\ldots&-\tilde{\mathsf{R}}_{1}\odot\tilde{\mathsf{R}}_{d}\odot\mathsf{K}\\\
\vdots&\vdots&\ddots&\vdots\\\
-\tilde{\mathsf{R}}_{d}\odot\mathsf{K}&-\tilde{\mathsf{R}}_{1}\odot\tilde{\mathsf{R}}_{d}\odot\mathsf{K}&\ldots&\left(\bar{\textbf{1}}-\tilde{\mathsf{R}}_{d}\odot\tilde{\mathsf{R}}_{d}\right)\odot\mathsf{K}\\\
\end{bmatrix},$ (23)

where $\bar{\textbf{1}}=\textbf{1}_{n_{x}\times n_{x}}$,
$\tilde{\mathsf{R}}_{i}=\gamma_{i}\mathsf{R}_{i}(\mathsf{X})=\mathsf{R}_{i}(\tilde{X})$,
and $\tilde{X}=\mathsf{X}\mathsf{P}$. The modified gradient-enhanced
covariance matrix is

$\tilde{\Sigma}_{\nabla}(\mathsf{X};\boldsymbol{\gamma},\eta_{\tilde{\mathsf{K}}_{\nabla}})=\sigma_{\mathsf{K}}^{2}\left(\tilde{\mathsf{K}}_{\nabla}(\mathsf{X};\boldsymbol{\gamma})+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}\right).$
(24)

(a) Correlations in one dimension

(b) Correlation between $\frac{\partial
f(\boldsymbol{x})}{\partial\tilde{x}_{1}}$ and $\frac{\partial
f(\boldsymbol{y})}{\partial\tilde{x}_{2}}$

Figure 2: Correlations for the Gaussian kernel from Eq. (1) for the function
evaluations of a generic function $f(\boldsymbol{x})$ and its gradient. The
term $k(\tilde{r})$ is the correlation between $f(\boldsymbol{x})$ and
$f(\boldsymbol{y})$, $\frac{\partial k(\tilde{r})}{\partial\tilde{r}}$ is the
correlation between $f(\boldsymbol{x})$ and $\frac{\partial
f(\boldsymbol{y})}{\partial x}$, and the correlation between $\frac{\partial
f(\boldsymbol{x})}{\partial x_{i}}$ and $\frac{\partial
f(\boldsymbol{y})}{\partial x_{j}}$ is
$-\frac{\partial^{2}k(\tilde{r})}{\partial\tilde{r}_{i}\partial\tilde{r}_{j}}$.

The correlations for the entries in $\boldsymbol{f}_{\nabla}$ can be seen in
Fig. 2. Having $\tilde{\mathsf{K}}_{\nabla}$ as a correlation matrix makes the
GP easier to interpret. Values in $\tilde{\mathsf{K}}_{\nabla}$ close to $-1$
or $1$ indicate near perfect inverse or direct correlation, respectively. On
the other hand, values close to $-1$ and $1$ in $\mathsf{K}_{\nabla}$ indicate
negative and positive relations, respectively, but provide little insight on
the strength of the relations between the evaluation points.

##  5 Bounding the condition number of the covariance matrix

###  5.1 The use of a nugget

A common approach to alleviate the ill-conditioning of a matrix is to add a
nugget to its diagonal. For a GP, the addition of a nugget to the covariance
matrix is analogous to having noisy data rasmussen_gaussian_2006 ;
ameli_noise_2022 . When the nugget is zero, the surrogate from the GP will
match the function of interest exactly at all points where it has been
evaluated. The same applies to the evaluated gradients if a gradient-enhanced
covariance matrix is used. However, if a positive nugget is used, the
surrogate will generally not match the function of interest exactly at points
where it has been sampled. It is therefore desirable to use the smallest
nugget value required to ensure the condition number of the covariance matrix
is below a desired threshold.

The eigenvalues of $\mathsf{K}$, $\mathsf{K}_{\nabla}$, and
$\tilde{\mathsf{K}}_{\nabla}$ are real since they are symmetric matrices. We
derive the required nugget to have
$\kappa(\mathsf{K}+\eta_{\mathsf{K}}\mathsf{I})\leq\kappa_{\max}$ when the
condition number is based on the $\ell_{2}$ norm:

$\displaystyle\kappa(\mathsf{K}+\eta_{\mathsf{K}}\mathsf{I})$
$\displaystyle=\frac{\lambda_{\text{max}}+\eta_{\mathsf{K}}}{\lambda_{\text{min}}+\eta_{\mathsf{K}}}\leq\frac{\lambda_{\text{max}}}{\eta_{\mathsf{K}}}+1\leq\kappa_{\max}$
$\displaystyle\eta_{\mathsf{K}}$
$\displaystyle\geq\frac{\lambda_{\text{max}}}{\kappa_{\max}-1},$ (25)

where $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ are the smallest and
largest eigenvalues of $\mathsf{K}$, respectively. The results are analogous
for $\mathsf{K}_{\nabla}$ and $\tilde{\mathsf{K}}_{\nabla}$. For positive
semidefinite matrices, such as $\mathsf{K}$, $\mathsf{K}_{\nabla}$, and
$\tilde{\mathsf{K}}_{\nabla}$, we have
$\lambda_{\text{max}}(\mathsf{K})\leq\text{tr}\left(\mathsf{K}\right)$. From
Eq. (25) it thus follows that sufficient nugget values to bound the condition
numbers of the kernel matrices below $\kappa_{\max}$ are

$\displaystyle\eta_{\mathsf{K}}$ $\displaystyle=\frac{n_{x}}{\kappa_{\max}-1}$
(26) $\displaystyle\eta_{\mathsf{K}_{\nabla}}(\boldsymbol{\gamma})$
$\displaystyle=\frac{n_{x}(\textbf{1}^{\top}\boldsymbol{\gamma}^{2}+1)}{\kappa_{\max}-1}$
(27) $\displaystyle\eta_{\tilde{\mathsf{K}}_{\nabla}}$
$\displaystyle=\frac{n_{x}(d+1)}{\kappa_{\max}-1}.$ (28)

These are sufficient but not necessary conditions to ensure that the condition
numbers of the covariance matrices are smaller than $\kappa_{\max}$ since the
bound ${\lambda_{\text{max}}<\text{tr}\left(\mathsf{K}\right)}$ is not tight
and Eq. (25) was derived with the worst case ${\lambda_{\text{min}}=0}$.

Eq. (26) provides a sufficiently small $\eta_{\mathsf{K}}$ to ensure that
${\kappa(\Sigma)\leq\kappa_{\max}}$. However, $\eta_{\mathsf{K}_{\nabla}}$
from Eq. (27) is undesirable since it depends on $\boldsymbol{\gamma}$.
Consequently, as $\boldsymbol{\gamma}$ gets larger,
$\eta_{\mathsf{K}_{\nabla}}$ must also get bigger to ensure that
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\leq\kappa_{\max}$. The method
from marchildon_non-intrusive_2023  to address this is summarized in Section
5.3. Finally, $\eta_{\tilde{\mathsf{K}}_{\nabla}}$ from Eq. (28) ensures that
$\kappa(\tilde{\Sigma}_{\nabla})\leq\kappa_{\max}$ with no dependence on
$\boldsymbol{\gamma}$ and could be used on its own. However, a tighter bound
on $\eta_{\tilde{\mathsf{K}}_{\nabla}}$ is derived in Section 5.4. that scales
with $\mathcal{O}(n_{x}\sqrt{d})$ instead of $\mathcal{O}(n_{x}d)$ when Eq.
(28) is used. In the remaining subsections three methods are introduced that
bound the condition number of the gradient-enhanced covariance matrix.

###  5.2 Baseline method: constrained optimization of $\boldsymbol{\gamma}$

One previous approach that has been used to ensure that
$\kappa(\Sigma_{\nabla})\leq\kappa_{\max}$ is to add a constraint to the
maximization of the marginal log-likelihood from Eq. (17) won_maximum_2006 .
The hyperparameters $\boldsymbol{\gamma}$ are thus selected by solving the
following constrained optimization problem:

$\boldsymbol{\gamma}^{*}=\operatornamewithlimits{argmax}_{\boldsymbol{\gamma}>0}L(\boldsymbol{\gamma})\quad\text{s.t.}\quad\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\leq\kappa_{\max}.$
(29)

There will always be a feasible solution to Eq. (29) if
$\eta_{\mathsf{K}_{\nabla}}\geq\frac{n_{x}}{\kappa_{\max}-1}$. This can be
verified from Eq. (27) with $\boldsymbol{\gamma}\rightarrow\textbf{0}_{d}$.
Solving Eq. (29) to set the hyperparameters thus ensures that
$\kappa(\Sigma_{\nabla}(\mathsf{X},\boldsymbol{\gamma}))\leq\kappa_{\max}\forall\,\mathsf{X}\in\mathbb{R}^{n_{x}\times
d}$. However, the constraint in Eq. (29) may result in selecting
hyperparameters that provide a significantly lower marginal log-likelihood.
This impacts the accuracy of the surrogate and is shown in Section 7.3 to be
detrimental to the optimization.

###  5.3 Rescaling method

A short overview of the rescaling method from marchildon_non-intrusive_2023
is provided in this section; the proofs can be found in marchildon_non-
intrusive_2023 . To ensure that
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma};\eta_{\mathsf{K}_{\nabla}}))\leq\kappa_{\max}$
when $\gamma_{1}=\ldots=\gamma_{d}$ and $\Sigma_{\nabla}$ is not diagonally
dominant, the parameter space is rescaled such that the minimum Euclidean
distance between evaluation points is $v_{\text{min,set}}$, where

$\displaystyle v_{\text{min,set}}(d,n_{x})$
$\displaystyle=\min\left(2\sqrt{d},\frac{2+\sqrt{4+2e^{2}\ln\left(\frac{(n_{x}-1)(1+2\sqrt{d})}{2}\right)}}{e}\right).$
(30)

The condition that $\gamma_{1}=\ldots=\gamma_{d}$ can be achieved by rescaling
the data non-isotropically marchildon_non-intrusive_2023 . Meanwhile, the
condition that $\Sigma_{\nabla}$ is not diagonally dominant is required since
the condition number of $\Sigma_{\nabla}$ is unbounded as
$\boldsymbol{\gamma}$ tends to infinity, regardless of the selected nugget.
However, since the condition on the diagonal dominance of $\Sigma_{\nabla}$
only applies for large values of $\boldsymbol{\gamma}$, it is not in practice
limiting since there is little correlation between evaluation points and thus
the marginal log-likelihood is unlikely to be maximized at these values of
$\boldsymbol{\gamma}$ marchildon_non-intrusive_2023 .

For the initial isotropic rescaling we use

$\displaystyle\mathsf{X}$ $\displaystyle=\tau\mathsf{X}_{\text{initial}}$
$\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_{i}}$
$\displaystyle=\frac{1}{\tau}\left(\frac{\partial f(\boldsymbol{x})}{\partial
x_{i}}\right)_{\text{initial}},$

where $\tau=\frac{v_{\text{min,set}}}{v_{\text{min,initial}}}$. The required
nugget to bound the condition number of $\Sigma_{\nabla}$ with the rescaling
method is

$\eta_{\mathsf{K}_{\nabla}}(d,n_{x})=\frac{1+(n_{x}-1)\frac{2\sqrt{d}}{v_{\text{req}}}e^{\frac{v_{\text{req}}}{2\sqrt{d}}-1}}{\kappa_{\max}-1}.$
(31)

The hyperparameters are set by solving the optimization problem from Eq. (29).
Since the bound
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma};\eta_{\mathsf{K}_{\nabla}}))\leq\kappa_{\max}$
is only ensured when $\gamma_{1}=\ldots=\gamma_{d}$, the maximization of the
marginal log-likelihood needs to contain a constraint on the condition number
to ensure that it does not exceed $\kappa_{\max}$, which could cause the
Cholesky decomposition to fail. If the constraint is not active, even if all
the hyperparameters are not equal, then no further iterations are required and
the optimized hyperparameters can be used. However, if the constraint on the
condition number is active, then a non-isotropic rescaling of $\mathsf{X}$ and
the gradients can be performed to get an unconstrained solution to Eq. (29).
The details on the non-isotropic rescaling can be found in marchildon_non-
intrusive_2023 .

###  5.4 Preconditioning method with
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$

Eq. (28) can be used to select $\eta_{\tilde{\mathsf{K}}_{\nabla}}$ but it
results in $\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}d)$. In this
section a smaller sufficient nugget is derived such that
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$ is sufficient
to ensure that
$\kappa(\tilde{\Sigma}_{\nabla}(\boldsymbol{\gamma};\eta_{\tilde{\mathsf{K}}_{\nabla}}))\leq\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$.
The derivation uses the Gershgorin circle theorem, which bounds the largest
eigenvalue of a symmetric matrix $\mathsf{A}$ by

$\lambda_{\text{max}}(\mathsf{A})\leq\max_{i}(a_{ii}+\sum_{j\neq i}|a_{ij}|).$
(32)

The two following propositions provide an upper bound on the sum of the
absolute values of the off-diagonal entries of $\tilde{\mathsf{K}}_{\nabla}$
when it is constructed with the Gaussian kernel from Eq. (1).

###### Proposition 1

For $n_{x},d\in\mathbb{Z}^{+}$ the sum of the absolute values for the off-
diagonal entries for any of the first $n_{x}$ rows of
$\tilde{\mathsf{K}}_{\nabla}$ from Eq. (23) using the Gaussian kernel is
bounded by $u_{\text{G}}$, where

$u_{\text{G}}(n_{x},d)=(n_{x}-1)\frac{1+\sqrt{1+4d}}{2}e^{-\frac{1+2d-\sqrt{1+4d}}{4d}}.$
(33)

###### Proof

We derive an upper bound for the sum of the absolute values for the off-
diagonal entries for any of the first $n_{x}$ rows of
$\tilde{\mathsf{K}}_{\nabla}$ from Eq. (23). The derivation is the same for
any of the rows and we thus consider the $a$-th row, where $1\leq a\leq n_{x}$
and

$\displaystyle\sum_{\begin{subarray}{c}i=1\\\ i\neq
a\end{subarray}}^{n_{x}}\left(\tilde{\mathsf{K}}_{\nabla}\right)_{ai}$
$\displaystyle=\sum_{\begin{subarray}{c}i=1\\\ i\neq
a\end{subarray}}^{n_{x}}\left(1+\sum_{j=1}^{d}\left|\tilde{x}_{aj}-\tilde{x}_{ij}\right|\right)e^{-\frac{\|\tilde{x}_{a:}-\tilde{x}_{i:}\|_{2}^{2}}{2}}$
$\displaystyle\leq(n_{x}-1)\max_{i}\left(1+\sum_{j=1}^{d}\left|\tilde{x}_{aj}-\tilde{x}_{ij}\right|\right)e^{-\frac{\|\tilde{x}_{a:}-\tilde{x}_{i:}\|_{2}^{2}}{2}}$
$\displaystyle\leq(n_{x}-1)\max_{\boldsymbol{w}\geq
0}\left(\left(1+\boldsymbol{w}^{\top}\textbf{1}_{d}\right)e^{-\frac{\boldsymbol{w}^{\top}\boldsymbol{w}}{2}}\right),$

where $w_{j}=|\tilde{x}_{aj}-\tilde{x}_{ij}|$ and we denote the expression
inside the max function as $g(\boldsymbol{w})$. To identify the maximum of
$g(\boldsymbol{w})$ we calculate its derivative and set it to zero:

$\displaystyle\frac{\partial g(\boldsymbol{w})}{\partial w_{i}}$
$\displaystyle=\left(1-w_{i}(1+\boldsymbol{w}^{\top}\textbf{1}_{d})\right)e^{-\frac{\boldsymbol{w}^{\top}\boldsymbol{w}}{2}}=0$
$\displaystyle w_{i}$
$\displaystyle=\frac{1}{1+\boldsymbol{w}^{\top}\textbf{1}}\,\forall\,i\in\\{1,\ldots,d\\}.$

It is clear that the gradient of $g(\boldsymbol{w})$ is zero if and only if
all of the entries in $\boldsymbol{w}$ are equal. We thus use
$\boldsymbol{w}=\alpha\textbf{1}_{d}$ and solve for the value of $\alpha$ that
maximizes $g(\alpha\textbf{1}_{d})$:

$\displaystyle\frac{\partial
g(\alpha\textbf{1}_{d})}{\partial\alpha}=d\left(1-\alpha(1+d\alpha)\right)e^{-\frac{d\alpha^{2}}{2}}$
$\displaystyle=0$ $\displaystyle\alpha^{*}$
$\displaystyle=\frac{-1+\sqrt{1+4d}}{2d},$

where we only kept the positive root since $\boldsymbol{w}\geq 0$ and it is
straightforward to verify that this provides the maximum of
$g(\boldsymbol{w})$. Eq. (33) is recovered by evaluating
$g(\alpha^{*}\textbf{1}_{d})$, which completes the proof.

###### Proposition 2

The sum of the absolute values for the off-diagonal entries for any of the
last $n_{x}d$ rows of $\tilde{\mathsf{K}}_{\nabla}$ using the Gaussian kernel
is smaller than $u_{\text{G}}(n_{x},d)$ from Eq. (33) for
$n_{x},d\in\mathbb{Z}^{+}$.

###### Proof

The proof can be found in Section A.1.

As a result of Propositions 1 and 2 and the Gershgorin circle theorem we have
$\lambda_{\text{max}}(\tilde{\mathsf{K}}_{\nabla})\leq
1+u_{\text{G}}(n_{x},d)$, where $u_{\text{G}}(n_{x},d)$ is from Eq. (33).
Therefore, we can ensure that
$\kappa(\tilde{\mathsf{K}}_{\nabla}(\boldsymbol{\gamma})+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I})\leq\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$
with Eqs. (25) and (33):

$\eta_{\tilde{\mathsf{K}}_{\nabla}}(n_{x},d)=\frac{1+(n_{x}-1)\frac{1+\sqrt{1+4d}}{2}e^{-\frac{1+2d-\sqrt{1+4d}}{4d}}}{\kappa_{\max}-1}.$
(34)

The following lemma proves that Eq. (34) provides
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$.

###### Lemma 1

From Eq. (34) we have
$\eta_{\tilde{\mathsf{K}}_{\nabla}}(n_{x},d)=\mathcal{O}(n_{x}\sqrt{d})$ for
$n_{x},d\in\mathcal{Z}^{+}$ since

$\eta_{\tilde{\mathsf{K}}_{\nabla}}(n_{x},d)<\frac{1+(n_{x}-1)(1+\sqrt{d})e^{-\frac{3-\sqrt{5}}{4}}}{\kappa_{\max}-1}.$
(35)

###### Proof

The proof can be found in Section A.2.

From Lemma 1 it is clear that the use of Eq. (34) to calculate
$\eta_{\tilde{\mathsf{K}}_{\nabla}}$ is advantageous for high dimensional
problems since it provides
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$ instead of
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}d)$ from Eq. (28).

##  6 Implementation

Since the inverse of
$\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}$ is needed to
calculate $\mu_{f}$, $\sigma_{f}^{2}$, and $\ln(L)$ along with their
gradients, it is desirable to calculate its Cholesky decomposition. Once the
Cholesky decomposition has been calculated, it becomes inexpensive to evaluate
$\mu_{f}(\boldsymbol{x})$ and $\sigma_{f}^{2}(\boldsymbol{x})$ for various
$\boldsymbol{x}$. However, doing so directly may cause the decomposition to
fail since the condition number of
$\mathsf{K}_{\nabla}+\eta_{\mathsf{K}_{\nabla}}\mathsf{I}$ cannot be bounded
with a finite $\eta_{\mathsf{K}_{\nabla}}$, as explained in Section 5.1.
Instead, $\mathsf{K}_{\nabla}$ can be used to construct
$\tilde{\mathsf{K}}_{\nabla}$ and the Cholesky decomposition of
$\tilde{\mathsf{K}}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}$ is
calculated instead. Calculating the Cholesky decomposition of
$\tilde{\mathsf{K}}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}$ is
numerically stable since
$\kappa(\tilde{\mathsf{K}}_{\nabla}(\boldsymbol{\gamma})+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I})\leq\kappa_{\max}\,\forall\boldsymbol{\gamma}>0$,
as proven in Section 5. The following algorithm details how the Cholesky
decomposition for the gradient-enhanced covariance matrix should be
calculated.

Algorithm 1 Stable Cholesky decomposition for gradient-enhanced GP

1:Select evaluation points $\mathsf{X}$ and hyperparameters
$\boldsymbol{\gamma}$

2:Calculate $\mathsf{K}_{\nabla}$, $\mathsf{P}$, and
$\eta_{\tilde{\mathsf{K}}_{\nabla}}$ with Eqs. (8), (21), and (34),
respectively

3:$\tilde{\mathsf{K}}_{\nabla}=\mathsf{P}^{-1}\mathsf{K}_{\nabla}\mathsf{P}^{-1}$

4:$\tilde{\mathsf{L}}\tilde{\mathsf{L}}^{\top}=\tilde{\mathsf{K}}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}$

5:$\mathsf{L}=\mathsf{P}\tilde{\mathsf{L}}$

Algebraically we have the following relation:

$\displaystyle\mathsf{L}\mathsf{L}^{T}=\left(\mathsf{K}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{P}\mathsf{P}\right)^{-1}=\mathsf{P}^{-1}\left(\tilde{\mathsf{K}}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}\right)^{-1}\mathsf{P}^{-1}.$

This indicates that the addition of the nugget
$\eta_{\tilde{\mathsf{K}}_{\nabla}}$ after the preconditioning is equivalent
to adding a nugget that varies along the diagonal and scales with the squared
hyperparameter $\boldsymbol{\gamma}$. However, the condition number of
$\left(\tilde{\mathsf{K}}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{I}\right)$
is bounded below $\kappa_{\max}$ for $\boldsymbol{\gamma}>0$ and can be
several orders of magnitude smaller than the condition number for
$\left(\mathsf{K}_{\nabla}+\eta_{\tilde{\mathsf{K}}_{\nabla}}\mathsf{P}\mathsf{P}\right)$.
As such, the Cholesky decomposition of the former should be performed, as
detailed in Algorithm 1.

##  7 Results

###  7.1 Condition numbers of the Gaussian kernel covariance matrices

The condition numbers of the gradient-free and gradient-enhanced covariance
matrices are compared for the baseline method presented in Section 5.2, the
rescaling method from [15], and the preconditioning method introduced in this
paper. The covariance matrices depend only on the evaluation points in
$\mathsf{X}$, the nugget, and the hyperparameters $\boldsymbol{\gamma}$.
However, the maximization of the marginal log-likelihood also depends on the
function of interest and we use the Rosenbrock function:

$f(\boldsymbol{x})=\sum_{i=1}^{d-1}\left[10\left(x_{i+1}-x_{i}^{2}\right)^{2}+\left(1-x_{i}\right)^{2}\right].$
(36)

A Latin hypercube sampling centred around
$\boldsymbol{x}=[1,\ldots,1]^{\top}$, which is the minimum for the Rosenbrock
function, is used to select the evaluation points

$\mathsf{X}=10^{-3}\times\begin{bmatrix}\phantom{-}1&\phantom{-}9&\phantom{-}7&-9&-5&-7&-3&\phantom{-}5&\phantom{-}3&-1\\\
\phantom{-}1&-3&\phantom{-}7&\phantom{-}3&\phantom{-}5&-9&-7&\phantom{-}9&-1&-5\end{bmatrix}^{\top}+1,$
(37)

where $v_{\text{min}}=\sqrt{2}/500\approx 2.8\times 10^{-3}$, which is the
minimum Euclidean distance between evaluation points.

(a) Baseline gradient-free: $\log(\kappa(\Sigma))$

(b) Baseline gradient-enhanced: $\log(\kappa(\Sigma_{\nabla}))$

(c) Gradient-enhanced with the rescaling method:
$\log(\kappa(\Sigma_{\nabla}))$

(d) Gradient-enhanced with the preconditioning method:
$\log(\kappa(\tilde{\Sigma}_{\nabla}))$

Figure 3: The condition number for covariance matrices using the nugget
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=1.5\times 10^{-9}$ from Eq. (34) along
with the set of evaluation points $\mathsf{X}$ from Eq. (37), which has a
minimum Euclidean distance between evaluation points of
$v_{\text{min}}=2.8\times 10^{-3}$. The star marking indicates where the
marginal log-likelihood function from Eq. (17) is maximized with the use of
the Rosenbrock function from Eq. (36).

Fig. 3 plots the condition number of the covariance matrices as a function of
$\boldsymbol{\gamma}$ for the evaluation points from Eq. (37). The star marker
indicates where the marginal log-likelihood from Eq. (17) is maximized. The
nugget value for all cases is $\eta=1.5\times 10^{-9}$, which comes from Eq.
(34). Red regions in Fig. 3 indicate where the condition number is greater
than $\kappa_{\max}=10^{10}$.

Figs. 3(a) and 3(b) plot the condition number of the gradient-free and
gradient-enhanced covariance matrices, respectively, using the baseline
method, which does not precondition the covariance matrix but adds the nugget
to its diagonal. For the gradient-free case we have
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$.
However, for the gradient-enhanced case
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\geq\kappa_{\max}$ for most
values of $\boldsymbol{\gamma}$, including where the marginal log-likelihood
from Eq. (17) is maximized. Selecting the hyperparameters to satisfy the
constraint $\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\leq\kappa_{\max}$
results in a lower marginal log-likelihood. This impacts the accuracy of the
surrogate, which degrades the performance of the Bayesian optimizer, as will
be shown in Section 7.3.

Fig. 3(c) plots $\log(\kappa(\Sigma_{\nabla}))$ with the use of the rescaling
method. As stated in Section 5.3, the rescaling method only ensures that
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\leq\kappa_{\max}$ when
$\gamma_{1}=\ldots=\gamma_{d}$ and $\Sigma_{\nabla}$ is not diagonally
dominant. While there are several values of $\boldsymbol{\gamma}$ where
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\geq\kappa_{\max}$, the condition
number is below $\kappa_{\max}$ where the marginal log-likelihood is
maximized.

From Fig. 3(d) we have
$\kappa(\tilde{\mathsf{K}}_{\nabla})<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$
as a result of the preconditioning method. Consequently, there is no need for
a constraint on the condition number for the optimization of the
hyperparameters. This ensures that the hyperparameters are never constrained
by the need to bound the condition number of the covariance matrix.
Furthermore, this also provides a small reduction in the cost of the
hyperparameter optimization since the constraint and its gradient do not need
to be calculated.

The blue in Figs. 3(a) and 3(d) indicate regions where $\Sigma$ and
$\tilde{\Sigma}_{\nabla}$ are nearly equal to the identity matrix. While the
condition number is very small in these regions, the marginal log-likelihood
is as well. This makes it undesirable to select those values of
$\boldsymbol{\gamma}$.

###  7.2 Applications to other kernels

The preconditioning method can be applied to gradient-enhanced covariance
matrices that utilize kernels other than the Gaussian kernel considered thus
far. For example, the preconditioning method can be applied to the Matérn
$\frac{5}{2}$ and rational quadratic kernels rasmussen_gaussian_2006 :

$\displaystyle k_{\text{M}\frac{5}{2}}(\tilde{\boldsymbol{r}})$
$\displaystyle=\left(1+\sqrt{3}\|\tilde{\boldsymbol{r}}\|+\|\tilde{\boldsymbol{r}}\|^{2}\right)e^{-\sqrt{3}\|\tilde{\boldsymbol{r}}\|}$
(38) $\displaystyle k_{rq}(\tilde{\boldsymbol{r}})$
$\displaystyle=\left(1+\frac{\|\tilde{\boldsymbol{r}}\|^{2}}{2\alpha}\right)^{-\alpha},$
(39)

where $\alpha>0$ is a hyperparameter and $\tilde{r}=\gamma_{i}(x_{i}-y_{i})$.
The hyperparameters for the Matérn $\frac{5}{2}$ and rational quadratic
kernels from Eqs. (38) and (39), respectively, have been selected such that
the preconditioning matrix required to make their gradient-enhanced kernel
matrix a correlation matrix also comes from Eq. (20). Using Eq. (27) provides
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}d)$ and ensures that
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$.
Alternatively, Eq. (34) could be used to have
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$, but this does
not provide a provable upper bound for $\kappa(\tilde{\Sigma}_{\nabla})$ since
Eq. (34) is specific to the Gaussian kernel. The methodology used in Section
5.4 to derive a nugget
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$ that ensures
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$
could also be applied to other kernels, such as the Matérn $\frac{5}{2}$ and
rational quadratic kernels.

Fig. 4 plots the condition number of the baseline and preconditioned gradient-
enhanced covariance matrices for the Matérn $\frac{5}{2}$ and rational
quadratic kernels. The set of evaluation points in $\mathsf{X}$ and the nugget
come from Eqs. (37) and (34), respectively, which are the same as the ones
used for Fig. 3. It is clear from Figs. 4(a) and 4(c) that the condition
number for the baseline method, i.e. the non-preconditioned gradient-enhanced
covariance matrices, for both kernels is larger than $\kappa_{\max}$ for
several values of $\boldsymbol{\gamma}$, including where the marginal log-
likelihood is maximized at the star marker. However, with the preconditioning
method we have
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$
for both kernels as seen in Figs. 4(b) and 4(d). This demonstrates that the
gradient-enhanced covariance matrix constructed with various kernels suffers
from severe ill-conditioning. Fortunately, the preconditioning method can be
applied to bound the condition number of $\tilde{\Sigma}_{\nabla}$ constructed
with various kernels. The provable bound
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$
requires the nugget to be calculated from Eq. (27), but using Eq. (34) was
sufficient for the case considered in Fig. 4. A user of a gradient-enhanced GP
with a non-Gaussian kernel could start by using a nugget calculated with Eq.
(34), and switch to using Eq. (27) instead if the condition number of
$\tilde{\Sigma}_{\nabla}$ was found to exceed $\kappa_{\max}$.

(a) Baseline Matérn $\frac{5}{2}$ kernel: $\log(\kappa(\Sigma_{\nabla}))$

(b) Preconditioned Matérn $\frac{5}{2}$ kernel:
$\log(\kappa(\tilde{\Sigma}_{\nabla}))$

(c) Baseline rational quadratic kernel: $\log(\kappa(\Sigma_{\nabla}))$

(d) Preconditioned rational quadratic kernel:
$\log(\kappa(\tilde{\Sigma}_{\nabla}))$

Figure 4: The condition number of the gradient-enhanced covariance matrix with
the baseline and preconditioned methods. All cases use Eq. (37) for the set of
evaluation points $\mathsf{X}$ and Eq. (34) to calculate
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=1.5\times 10^{-9}$. The value of the
hyperparameters $\boldsymbol{\gamma}$ that maximizes the marginal log-
likelihood function from Eq. (17) with the Rosenbrock function from Eq. (36)
is indicated by the star marker.

###  7.3 Optimization

In this section, the baseline, rescaling, and preconditioning methods are
compared when used with a Bayesian optimizer to minimize the Rosenbrock
function from Eq. (36) with $d\in\\{2,5,10,15\\}$. For each test case, the
optimization is repeated five separate times for each method. The starting
points are selected with the Latin hypercube sampling from the open source
Surrogate Modeling Toolbox with a lower bound of $-10$ and an upper bound of
$10$ for the parameters, and the random state set to 1. This ensures that the
optimizer for each method starts from the same initial solution. The selected
acquisition function is the upper-confidence function

$h(\boldsymbol{x})=\mu_{f}(\boldsymbol{x})+\omega\sigma_{f}(\boldsymbol{x}),$
(40)

where $\omega\geq 0$. The parameter $\omega$ promotes exploitation when it is
small, and exploration when it is large. We use $\omega=0$ since we are
interested in local optimization for the unimodal Rosenbrock function. The
gradient-based SLSQP optimizer from the Python library SciPy is used to select
the hyperparameters by maximizing the marginal log-likelihood. The same
optimizer is used to minimize the acquisition function to select the next
point in the parameter space to evaluate the Rosenbrock function. A trust
region is used in the minimization of the acquisition function, similar to the
one used in eriksson_scalable_2019 , where a Bayesian optimizer was also used
for local minimization. However, our trust region is set to be a hypersphere
instead of a hyperrectangle.

(a) $d=2$

(b) $d=5$

(c) $d=10$

(d) $d=15$

Figure 5: Bayesian optimization of the Rosenbrock function from Eq. (36) using
the baseline, the rescaling, and preconditioning methods for
$d\in\\{2,5,10,15\\}$. The plots show the lowest evaluated optimality, i.e.
the $\ell_{2}$ norm of the gradient of Eq. (36), for each optimization run at
each iteration.

In Fig. 5 the optimality, which is the $\ell_{2}$ norm of the gradient, is
compared for the Bayesian optimizer using the baseline, rescaling, and
preconditioning methods. There are two important observations from these
optimality plots: the depth and rate of convergence of the optimality. In all
cases, the rescaling and preconditioning methods converge the optimality
several orders of magnitude deeper than the baseline method. In fact, the
rescaling and preconditioning methods converge the optimality to below
$10^{-12}$ in all test cases. Meanwhile, the deepest optimality that the
baseline method achieves is $10^{-7}$, and only for $d=2$. As the
dimensionality increases, the optimizer with the baseline method is not able
to converge the optimality as deeply and can only achieve an optimality of
$10^{-4}$ for the $d=15$ case. The optimizer with the rescaling and
preconditioning methods is thus able to converge the optimality 5 to 9
additional orders of magnitude relative to the optimizer with the baseline
method. For the baseline method, the hyperparameters $\boldsymbol{\gamma}$ are
selected by solving Eq. (29), where the marginal log-likelihood is maximized
with an upper bound on the condition number. As the optimality is converged,
the evaluation points get closer together in the parameter space and this
makes the ill-conditioning of the gradient-enhanced covariance matrix worse
marchildon_non-intrusive_2023 . Consequently, solving Eq. (29) results in
hyperparameters that provide a lower marginal log-likelihood since the upper
bound on the condition number becomes a more onerous constraint. The rescaling
and baseline methods do not suffer from this since, by construction, they
guarantee that the selected hyperparameters maximize the marginal log-
likelihood without being constrained by the condition number of the covariance
matrix.

It is clear from Fig. 5 that the optimizer utilizing the preconditioning
method achieves the fastest rate of convergence of the three methods for all
four test cases. The convergence of the optimality for the optimizer with the
rescaling method is significantly slower relative to the optimizer with the
baseline and preconditioning methods, particularly for the cases with $d=2$
and $d=5$. The slower convergence of the optimality for the optimizer using
the rescaling method was a trend that was also observed in marchildon_non-
intrusive_2023 . This trend was found to be a consequence of the rescaling
method providing a surrogate with gradients that have larger errors relative
to the baseline method.

In summary, the use of the preconditioning method with a gradient-enhanced
Bayesian optimizer enables the optimality to be converged more deeply than
with the use of the baseline method, and in fewer iterations than with the
rescaling method.

##  8 Conclusions

Method | Baseline | Rescale | Precondition
---|---|---|---
$\kappa(\Sigma_{\nabla}(\boldsymbol{\gamma}))\leq\kappa_{\max},\forall\boldsymbol{\gamma}>0$ | ✗ | $\gamma_{1}=\ldots=\gamma_{d}$ | ✓
Constraint free hyperparameter optz | ✗ | ✗ | ✓
Nodes can be collocated | ✓ | ✗ | ✓
Deep convergence: optimality $<10^{-8}$ | ✗ | ✓ | ✓
Provides a correlation matrix | ✗ | ✗ | ✓
Bounded $\kappa(\tilde{\Sigma}_{\nabla}(\boldsymbol{\gamma}))$ for other kernels | ✗ | ✗ | ✓
Table 1: Comparison of methods to address the ill-conditioning of the
covariance matrix $\Sigma_{\nabla}$. The baseline and rescaling methods are
summarized in Sections 5.2 and 5.3, respectively, and the implementation of
the preconditioning method is provided in Section 6.

A gradient-enhanced GP provides a more accurate probabilistic surrogate than
its gradient-free counterpart but the ill-conditioning of its covariance
matrix has been a hindrance to its use. A straightforward method has been
developed that overcomes this problem and ensures that the condition number of
the gradient-enhanced covariance matrix is always smaller than the user-set
threshold of $\kappa_{\max}$. The simple implementation is detailed in
Algorithm 1, which is found in Section 6. The method applies a diagonal
preconditioner along with a modest nugget that scales as
$\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}\sqrt{d})$ for the
Gaussian kernel, and $\eta_{\tilde{\mathsf{K}}_{\nabla}}=\mathcal{O}(n_{x}d)$
for other kernels. A tighter bound for $\eta_{\tilde{\mathsf{K}}_{\nabla}}$
for the non-Gaussian kernels will be considered in future work.

The benefits of using the preconditioning method relative to the baseline and
rescaling methods are summarized in Table 1. With the preconditioning method,
all of the data points can be kept and there is no minimum distance
requirement between evaluation points in the parameter space. The points can
even be collocated, unlike the rescaling method. Since the preconditioning
method ensures that
$\kappa(\Sigma(\boldsymbol{\gamma}))<\kappa_{\max}\,\forall\,\boldsymbol{\gamma}>0$,
no constraint is required when maximizing the marginal log-likelihood. This
simplifies the optimization and reduces its computational cost. The
preconditioning method also provides a correlation matrix, which makes the GP
easier to interpret.

In Section 7.3 the Rosenbrock function was optimized for $d\in\\{2,5,10,15\\}$
with a Bayesian optimizer using the baseline, rescaling and preconditioning
methods. The Bayesian optimizer with the preconditioning method converged the
optimality an additional 5-9 orders of magnitude relative to the optimizer
with the baseline method. Furthermore, the preconditioning method enabled the
Bayesian optimizer to converge the optimality more quickly than when the
rescaling method is used, particularly for the lower-dimensional problems. The
slower convergence of a Bayesian optimizer using the rescaling method was
previously identified to be the result of its surrogate having gradients with
larger errors. In conclusion, the preconditioning method bounds the condition
number of the preconditioned gradient-enhanced covariance matrix and it
enables a Bayesian optimizer to achieve deeper and faster convergence relative
to the use of either the baseline or rescaling methods.

The baseline, rescaling, and preconditioning methods are all available in the
open source python library GpGradPy, which can be found at
https://github.com/marchildon/gpgradpy/tree/paper_precon. All of the figures
in this paper can be reproduced with this library.

##  Appendix A Proofs

###  A.1 Proof for Proposition 2

The derivation of the upper bound for the sum of the absolute value of the
off-diagonal entries is the same for any of the last $n_{x}d$ rows of
$\tilde{\mathsf{K}}_{\nabla}$, which comes from Eq. (23). Without loss of
generality, we consider the $b$-th row of $\tilde{\mathsf{K}}_{\nabla}$, where
$b=pn_{x}+m$, and $p$ and $m$ can take any integer values that satisfy $1\leq
p\leq d$ and $1\leq m\leq n_{x}$:

$\displaystyle\sum_{\begin{subarray}{c}i=1\\\ i\neq
b\end{subarray}}^{n_{x}}\left(\tilde{\mathsf{K}}_{\nabla}\right)_{bi}$
$\displaystyle=\sum_{\begin{subarray}{c}i=1\\\ i\neq
m\end{subarray}}^{n_{x}}\left(\left|\tilde{x}_{mp}-\tilde{x}_{ip}\right|+\sum_{j=1}^{d}\left|\delta_{jp}-\left(\tilde{x}_{mp}-\tilde{x}_{ip}\right)\left(\tilde{x}_{mj}-\tilde{x}_{ij}\right)\right|\right)e^{-\frac{\|\tilde{\boldsymbol{x}}_{m:}-\tilde{\boldsymbol{x}}_{i:}\|_{2}^{2}}{2}}$
$\displaystyle\leq\sum_{\begin{subarray}{c}i=1\\\ i\neq
m\end{subarray}}^{n_{x}}\left(\left|\tilde{x}_{mp}-\tilde{x}_{ip}\right|+1+\sum_{\begin{subarray}{c}j=1\\\
j\neq
p\end{subarray}}^{d}\left|\left(\tilde{x}_{mp}-\tilde{x}_{ip}\right)\left(\tilde{x}_{mj}-\tilde{x}_{ij}\right)\right|\right)e^{-\frac{\|\tilde{\boldsymbol{x}}_{m:}-\tilde{\boldsymbol{x}}_{i:}\|_{2}^{2}}{2}}$
$\displaystyle=\sum_{\begin{subarray}{c}i=1\\\ i\neq
m\end{subarray}}^{n_{x}}\left(1+\left|\tilde{x}_{mp}-\tilde{x}_{ip}\right|\left(1+\sum_{\begin{subarray}{c}j=1\\\
j\neq
p\end{subarray}}^{n_{x}}\left|\tilde{x}_{mj}-\tilde{x}_{ij}\right|\right)\right)e^{-\frac{\|\tilde{\boldsymbol{x}}_{m:}-\tilde{\boldsymbol{x}}_{i:}\|_{2}^{2}}{2}}$
$\displaystyle\leq(n_{x}-1)\max_{\nu\geq 0,\,\check{\boldsymbol{w}}\geq
0}\left(1+\nu\left(1+\check{\textbf{1}}_{p}^{\top}\check{\boldsymbol{w}}\right)\right)e^{-\frac{1}{2}\left(\nu^{2}+\check{\boldsymbol{w}}^{\top}\check{\boldsymbol{w}}\right)},$
(41)

where $\nu=|\tilde{x}_{mp}-\tilde{x}_{ip}|$ and
$\check{w}_{j}=|\tilde{x}_{mj}-\tilde{x}_{ij}|$, except for $\check{w}_{p}=0$.
Similarly $\check{\textbf{1}}_{p}$ is a vector of ones of length $d$ with a
zero at its $p$-th entry. The first inequality is a result of
$\tilde{\mathsf{K}}_{\nabla}$ being a correlation matrix, as explained in
Section 4, and therefore
$|\frac{\partial^{2}\mathsf{K}}{\partial\tilde{x}_{a}\partial\tilde{y}_{a}}|\leq
1$.

An analogous approach to the one taken in Proposition 1 can be used to show
that the maximization of Eq. (41) requires
$\check{\boldsymbol{w}}=\alpha\check{\textbf{1}}_{p}$, i.e. that all but the
$p$-th entries in $\check{\boldsymbol{w}}$ are equal. Using
$\check{\boldsymbol{w}}=\alpha\check{\textbf{1}}_{p}$ with Eq. (41) gives

$g_{1}(\nu,\alpha;d)=\left(\nu+1+(d-1)\alpha\nu\right)e^{-\frac{\nu^{2}+(d-1)\alpha^{2}}{2}}.$
(42)

We thus need to prove that
$(n_{x}-1)g_{1}(\nu,\alpha;d)<u_{\text{G}}(n_{x},d)$ for $\nu,\alpha\geq 0$
and $d\in\mathbb{Z}^{+}$. The following lemma considers the case for $d=1$.

###### Lemma 2

For $d=1$ we have

$(n_{x}-1)\max_{\nu\geq 0,\alpha\geq
0}g_{1}(\nu,\alpha;d=1)=u_{\text{G}}(n_{x},d=1)=(n_{x}-1)\left(1+\sqrt{5}\right)e^{-\frac{3-\sqrt{5}}{4}},$
(43)

where $g_{1}(\nu,\alpha;d)$ comes from Eq. (42) and $u_{\text{G}}(n_{x},d)$
comes from Eq. (33).

###### Proof

For $d=1$ the parameter $\alpha$ cancels out and we thus have a scalar
function that we seek to maximize

$\displaystyle\frac{\partial g_{1}(\nu;d=1)}{\partial\nu}$
$\displaystyle=\frac{\partial\left(\left(\nu+1\right)e^{-\frac{\nu^{2}}{2}}\right)}{\partial\nu}$
$\displaystyle=-\left(\nu^{2}+\nu-1\right)e^{-\frac{\nu^{2}}{2}}=0$
$\displaystyle\nu^{*}_{d=1}$ $\displaystyle=\frac{-1+\sqrt{5}}{2},$

where only the positive root was kept since $\nu\geq 0$ and it is
straightforward to verify that this critical point maximizes
$g_{1}(\nu,\alpha;d=1)$. Eq. (43) is recovered by evaluating $g_{1}$ with
$\nu=\nu^{*}_{d=1}$ and $d=1$, which completes the proof.

To consider the cases for $d\geq 2$ we will need to find the value of $\alpha$
and $\nu$ that maximize $g_{1}(\nu,\alpha;d)$ from Eq. (42). The following
lemma considers the maximization of $g_{1}$ with respect to $\alpha$.

###### Lemma 3

For $\nu,\alpha\geq 0$ and $d\geq 2$ we have $g_{1}(\nu,\alpha;d)\leq
g_{2}(\nu;d)$, where $g_{1}$ comes from Eq. (42) and

$g_{2}(\nu;d)=\left(\frac{\nu+1+\sqrt{h_{1}(\nu;d)}}{2}\right)e^{-\frac{\nu^{2}}{2}+h_{2}(\nu;d)},$
(44)

where

$\displaystyle h_{1}(\nu;d)$ $\displaystyle=(\nu+1)^{2}+4\nu^{2}(d-1)$ (45)
$\displaystyle h_{2}(\nu;d)$
$\displaystyle=\frac{(\nu+1)\sqrt{h_{1}(\nu;d)}-(\nu+1)^{2}}{4\nu^{2}(d-1)}-\frac{1}{2}.$
(46)

###### Proof

To find the maximum of $g_{1}(\nu,\alpha;d)$ with respect to $\alpha$ we find
its derivative, set it to zero and solve for $\alpha$:

$\displaystyle\frac{\partial g_{1}(\nu,\alpha;d)}{\partial\alpha}$
$\displaystyle=-(d-1)\left((d-1)\nu\alpha^{2}+(\nu+1)\alpha-\nu\right)e^{-\frac{\nu^{2}+(d-1)\alpha^{2}}{2}}=0$
$\displaystyle\alpha^{*}$
$\displaystyle=\frac{-(\nu+1)+\sqrt{(\nu+1)^{2}+4\nu^{2}(d-1)}}{2\nu(d-1)},$

where only the positive root of the quadratic equation is kept since $\alpha$
must be positive and it is straightforward to verify that this provides the
maximum of $g_{1}(\nu,\alpha;d)$. The function $g_{2}(\nu;d)$ from Eq. (44) is
recovered by evaluating $g_{1}(\nu,\alpha^{*};d)$, which completes the proof.

Both $h_{1}(\nu;d)$ and $h_{2}(\nu;d)$ from Eqs. (45) and (46), respectively,
are non-polynomial functions that make it impractical to find a closed-form
maximum solution for $g_{2}(\nu;d)$. The two following lemmas provide upper
bounds for these non-polynomial functions.

###### Lemma 4

For $\nu\geq 0$ and $d\geq 2$ we have the bound $\sqrt{h_{1}(\nu;d)}\leq
h_{3}(\nu;d)$, where $h_{1}(\nu;d)$ comes from Eq. (45) and $h_{3}(\nu;d)$ is
the following $\mathbb{C}^{0}$ continuous piecewise polynomial:

$\displaystyle h_{3}(\nu;d)=\begin{cases}(2\sqrt{d}-1)\nu+1&\text{if
}0\leq\nu\leq 1\\\ 2\sqrt{d}\nu&\text{if }\nu>1.\end{cases}$ (47)

###### Proof

For $0\leq\nu\leq 1$ and $d\geq 2$ we start by showing that $h_{3}^{2}\geq
h_{1}$:

$\displaystyle\left((2\sqrt{d}-1)\nu+1\right)^{2}$
$\displaystyle\geq(\nu+1)^{2}+4\nu^{2}(d-1)$ $\displaystyle
4\nu(1-\nu)(\sqrt{d}-1)$ $\displaystyle\geq 0.$

Next we demonstrate that $h_{3}^{2}\geq h_{1}(\nu;d)$ for $\nu\geq 1$:

$\displaystyle\left(2\sqrt{d}\nu\right)^{2}$
$\displaystyle\geq(\nu+1)^{2}+4\nu^{2}(d-1)$
$\displaystyle\left(\nu+\frac{1}{3}\right)\left(\nu-1\right)$
$\displaystyle\geq 0.$

Finally, it is straightforward to verify that $h_{3}(\nu;d)$ is
$\mathbb{C}^{0}$ continuous:

$\lim_{\nu\rightarrow 1^{-}}h_{3}(\nu;d)=\lim_{\nu\rightarrow
1^{+}}h_{3}(\nu;d)=2\sqrt{d},$ (48)

which completes the proof.

###### Lemma 5

The maximum value for $h_{2}(\nu;d)$ from Eq. (46) for $\nu\geq 0$ and
$d\in\mathbb{Z}^{+}$ is

$\max_{\nu\geq 0}h_{2}(\nu;d)=\lim{\nu\rightarrow 0}h_{2}(\nu,d)=0.$ (49)

###### Proof

We start by proving that $h_{2}(\nu;d)$ is monotonically decreasing with
respect to $\nu$ by showing that its derivative with respect to $\nu$ is
nonpositive for $\nu\geq 0$ and $d\in\mathbb{Z}^{+}$

$\displaystyle\frac{\partial h_{2}(\nu;d)}{\partial\nu}$
$\displaystyle=-\frac{h_{1}-(\nu+1)\sqrt{h_{1}}-2\nu^{2}(d-1)}{2\nu^{3}\sqrt{h_{1}}(d-1)}.$

Since the denominator of $\frac{\partial h_{2}(\nu;d)}{\partial\nu}$ is always
nonnegative for $\nu\geq 0$ and $d\in\mathbb{Z}^{+}$, we only need to show
that its numerator is nonpositive for the same range of parameters:

$\displaystyle h_{1}-(\nu+1)\sqrt{h_{1}}-2\nu^{2}(d-1)$ $\displaystyle\geq 0$
$\displaystyle\left(h_{1}-2\nu^{2}(d-1)\right)^{2}$
$\displaystyle\geq(\nu+1)^{2}a$ $\displaystyle 4\nu^{4}(d-1)^{2}$
$\displaystyle\geq 0.$

Since $h_{2}(\nu;d)$ is monotonically decreasing with respect to $\nu$ for
$\nu\geq 0$ and $d\in\mathbb{Z}^{+}$, its maximum is at $\nu=0$. To evaluate
$h_{2}(\nu;d)$ we use a limit and apply l’Hôpital’s rule twice:

$\displaystyle\lim_{\nu\rightarrow 0}h_{2}(\nu;d)$
$\displaystyle=\lim_{\nu\rightarrow
0}\frac{(\nu+1)\sqrt{h_{1}}-(\nu+1)^{2}}{4(d-1)\nu^{2}}-\frac{1}{2}$
$\displaystyle=\lim_{\nu\rightarrow
0}\frac{2\sqrt{h_{1}}+\frac{4\nu(d-1)}{\sqrt{h_{1}}}-2(\nu+1)}{8\nu(d-1)}-\frac{1}{2}$
$\displaystyle=\lim_{\nu\rightarrow
0}\frac{\frac{4(\nu+1)(d-1)}{h_{1}^{\frac{3}{2}}}+\frac{(8d-6)\nu+2}{\sqrt{h_{1}}}-2}{8(d-1)}-\frac{1}{2}$
$\displaystyle=0,$

which completes the proof.

Thanks to Lemmas 4 and 5 it is now possible to derive a closed-form solution
for an upper bound of $g_{2}(\nu;d)$ from Eq. (44) for $\nu\geq 0$ and $d\geq
2$. This is considered in the two following lemmas that consider the case for
$0\leq\nu\leq 1$ and $\nu>1$, respectively.

###### Lemma 6

For $0\leq\nu\leq 1$ and $d\geq 2$ we have
$(n_{x}-1)g_{2}(\nu;d)<u_{\text{G}}(n_{x},d)$, where $g_{2}(\nu;d)$ and
$u_{\text{G}}(n_{x},d)$ come from Eqs. (44) and (33), respectively.

###### Proof

The function $g_{2}(\nu;d)$ from Eq. (44) contains the nonlinear functions
$h_{1}(\nu;d)$ and $h_{2}(\nu;d)$ from Eqs. (45) and (46), respectively. We
use the upper bounds provided by Lemmas 4 and 5 for these functions and
$0\leq\nu\leq 1$ to get $g_{2}(\nu;d)<g_{3}(\nu;d)$, where

$\displaystyle g_{3}(\nu;d)$
$\displaystyle=\frac{\nu+1+\left[(2\sqrt{d}-1)\nu+1\right]}{2}e^{-\frac{\nu^{2}}{2}}$
$\displaystyle=\left(\sqrt{d}\nu+1\right)e^{-\frac{\nu^{2}}{2}}.$ (50)

We now find the value of $\nu$ that maximizes $g_{3}(\nu;d)$

$\displaystyle\frac{\partial g_{3}}{\partial\nu}$
$\displaystyle=\left(\sqrt{d}-\nu\left(\sqrt{d}\nu+1\right)\right)e^{-\frac{\nu^{2}}{2}}=0$
$\displaystyle\nu_{3}^{*}$ $\displaystyle=\frac{-1+\sqrt{1+4d}}{2\sqrt{d}},$
(51)

where only the positive root was kept and it is straightforward to verify that
$0<\nu^{*}<1$ for $d\geq 2$, and that this is the maximum for the function
$g_{3}$. Using $\nu_{3}^{*}$ from Eq. (51) gives
$g_{3}(\nu_{3}^{*};d)=u_{\text{G}}(n_{x},d)$, where $u_{\text{G}}(n_{x},d)$
comes from Eq. (33). Therefore, we have for $d\geq 2$:

$\max_{0\leq\nu\leq 1}g_{2}(\nu;d)<\max_{0\leq\nu\leq
1}g_{3}(\nu;d)=g_{3}(\nu_{3}^{*};d)=u_{\text{G}}(n_{x},d),$ (52)

which completes the proof.

###### Lemma 7

We have $(n_{x}-1)g_{2}(\nu;d)<u_{\text{G}}(n_{x},d)$ for $\nu\geq 1$ and
$d\geq 2$, where $g_{2}(\nu;d)$ and $u_{\text{G}}(n_{x},d)$ come from Eqs.
(44) and (33), respectively.

###### Proof

We now consider the case for $\nu>1$ by substituting $h_{3}(\nu;d)$ from Eq.
(47) for $\nu>1$ into $g_{2}(\nu;d)$ for $\sqrt{h_{1}(\nu;d)}$ and using the
results from Lemma 5 for an upper bound on $h_{2}(\nu;d)$. We get the bound
$g_{2}(\nu;d)\leq g_{4}(\nu;d)$, where

$g_{4}(\nu;d)=\frac{\nu+1+\left[2\sqrt{d}\nu\right]}{2}e^{-\frac{\nu^{2}}{2}}.$
(53)

We now find the value of $\nu\geq 1$ that maximizes $g_{4}(\nu;d)$:

$\displaystyle\frac{\partial g_{4}}{\partial\nu}$
$\displaystyle=\frac{\left(2\sqrt{d}+1\right)-\nu\left(\left(2\sqrt{d}+1\right)\nu+1\right)}{2}e^{-\frac{\nu^{2}}{2}}=0$
$\displaystyle(2\sqrt{d}+1)\nu^{2}+\nu-(2\sqrt{d}+1)$ $\displaystyle=0$
$\displaystyle\nu_{4}^{*}$
$\displaystyle=\frac{-1+\sqrt{1+4\left(2\sqrt{d}+1\right)^{2}}}{2(2\sqrt{d}+1)},$

where only the positive root was kept and it is straightforward to show that
this provides the maximum for $g_{4}(\nu;d)$. However, we now demonstrate that
this root does not satisfy the constraint $\nu\geq 1$:

$\nu_{4}^{*}<\frac{-1+\left[1+2\left(2\sqrt{d}+1\right)\right]}{2(2\sqrt{d}+1)}=1,$

where we used the inequality $\sqrt{b_{1}+b_{2}}<\sqrt{b_{1}}+\sqrt{b_{2}}$
for $b_{1},b_{2}>0$. Since there are no roots for $\nu\geq 1$ that maximize
$g_{4}(\nu;d)$ for $d\geq 2$, it is either maximized at $\nu=1$ or
$\nu\rightarrow\infty$. For $\lim\nu\rightarrow\infty$ we have
$g_{4}(\nu;d)=0$ and for $\nu=1$ we have

$(n_{x}-1)g_{4}(\nu=1,d)=(n_{x}-1)g_{3}(\nu=1,d)<(n_{x}-1)g_{3}(\nu_{3}^{*},d)=u_{\text{G}}(n_{x},d),$
(54)

where $g_{3}=g_{4}$ for $\nu=1$ since both functions used the relation
$\sqrt{h_{1}(\nu;d)}\leq h_{3}(\nu;d)$ and it was shown in Lemma 4 that
$h_{3}(\nu)$ from Eq. (47) is $\mathbb{C}^{0}$ continuous. We thus have
$(n_{x}-1)g_{2}(\nu;d)<u_{\text{G}}(n_{x},d)$ for $\nu\geq 1$ and $d\geq 2$,
which completes the proof.

It has been proven that the function $g_{1}$ from Eq. (42), which provides an
upper bound for the sum of absolute values for the off-diagonal entries for
any of the last $n_{x}d$ rows of $\tilde{\mathsf{K}}_{\nabla}$, is smaller
than $u_{\text{G}}(n_{x},d)$ for $n_{x},d\in\mathbb{Z}^{+}$, which completes
the proof.

###  A.2 Proof for Lemma 1

We start by deriving an upper bound for the exponent for
$\eta_{\tilde{\mathsf{K}}_{\nabla}}$ in Eq. (34). To do this we take the
derivative of the exponent, which we denote as $g$, with respect to $d$:

$\displaystyle\frac{\partial g}{\partial d}$
$\displaystyle=\frac{\sqrt{1+4d}-1-2d}{4d^{2}\sqrt{1+4d}}$
$\displaystyle\leq\frac{\left[1+2\sqrt{d}\right]-1-2d}{4d^{2}\sqrt{1+4d}}$
$\displaystyle=\frac{\sqrt{d}-d}{2d^{2}\sqrt{1+4d}},$

which is always negative for $d\geq 1$. Therefore, $g(d)$ is monotonically
decreasing with respect to $d$ for $d\in\mathbb{(}Z)^{+}$ and is thus
maximized at $g(d=1)=-\frac{3-\sqrt{5}}{4}$. An upper bound for
$\eta_{\tilde{\mathsf{K}}_{\nabla}}$ from Eq. (34) is now derived

$\displaystyle\eta_{\tilde{\mathsf{K}}_{\nabla}}(n_{x},d)$
$\displaystyle=\frac{1+(n_{x}-1)\frac{1+\sqrt{1+4d}}{2}e^{-\frac{1+2d-\sqrt{1+4d}}{4d}}}{\kappa_{\max}-1}$
$\displaystyle<\frac{1+(n_{x}-1)\frac{1+\left[\sqrt{1}+\sqrt{4d}\right]}{2}e^{-\frac{3-\sqrt{5}}{4}}}{\kappa_{\max}-1}$
$\displaystyle=\frac{1+(n_{x}-1)(1+\sqrt{d})e^{-\frac{3-\sqrt{5}}{4}}}{\kappa_{\max}-1},$

where it is clear that
$\eta_{\tilde{\mathsf{K}}_{\nabla}}(n_{x},d)=\mathcal{O}(n_{x}\sqrt{d})$,
which completes the proof.

###### Acknowledgements.

The authors would like to thank the Natural Sciences and Engineering Research
Council of Canada and the Ontario Graduate Scholarship Program for their
financial support.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * (1) Ababou, R., Bagtzoglou, A.C., Wood, E.F.: On the condition number of covariance matrices in kriging, estimation, and simulation of random fields.  Mathematical Geology 26(1), 99–133 (1994).  DOI 10.1007/BF02065878
  * (2) Ameli, S., Shadden, S.C.: Noise Estimation in Gaussian Process Regression.  arXiv p. 41 (2022)
  * (3) Dalbey, K.: Efficient and robust gradient enhanced Kriging emulators.  Tech. Rep. SAND2013-7022, 1096451 (2013).  DOI 10.2172/1096451
  * (4) Davis, G.J., Morris, M.D.: Six Factors Which Affect the Condition Number of Matrices Associated with Kriging.  Mathematical Geology 29(5), 669–683 (1997).  DOI 10.1007/BF02769650
  * (5) De Roos, F., Gessner, A., Hennig, P.: High-Dimensional Gaussian Process Inference with Derivatives.  In: 38th International Conference on Machine Learning, pp. 2535–2545 (2021)
  * (6) Eriksson, D., Dong, K., Lee, E., Bindel, D., Wilson, A.G.: Scaling Gaussian Process Regression with Derivatives.  In: 32nd Conference on Neural Information Processing Systems. Montreal, Canada (2018)
  * (7) Eriksson, D., Pearce, M., Gardner, J., Turner, R.D., Poloczek, M.: Scalable Global Optimization via Local Bayesian Optimization.  In: 33rd Conference on Neural Information Processing Systems, p. 12. Vancouver, Canada (2019)
  * (8) Han, Z.H., Görtz, S., Zimmermann, R.: Improving variable-fidelity surrogate modeling via gradient-enhanced kriging and a generalized hybrid bridge function.  Aerospace Science and Technology 25(1), 177–189 (2013).  DOI 10.1016/j.ast.2012.01.006
  * (9) He, X., Chien, P.: On the Instability Issue of Gradient-Enhanced Gaussian Process Emulators for Computer Experiments.  SIAM/ASA Journal on Uncertainty Quantification 6(2), 627–644 (2018).  DOI 10.1137/16M1088247
  * (10) Higham, N.J.: Cholesky factorization.  Wiley Interdisciplinary Reviews: Computational Statistics 1(2), 251–254 (2009).  DOI 10.1002/wics.18
  * (11) Hung, T.H., Chien, P.: A Random Fourier Feature Method for Emulating Computer Models With Gradient Information.  Technometrics 63(4), 500–509 (2021).  DOI 10.1080/00401706.2020.1852973
  * (12) Kostinski, A.B., Koivunen, A.C.: On the condition number of Gaussian sample-covariance matrices.  IEEE Transactions on Geoscience and Remote Sensing 38(1), 329–332 (2000).  DOI 10.1109/36.823928
  * (13) Laurent, L., Le Riche, R., Soulier, B., Boucard, P.A.: An Overview of Gradient-Enhanced Metamodels with Applications.  Archives of Computational Methods in Engineering 26(1), 61–106 (2019).  DOI 10.1007/s11831-017-9226-3
  * (14) March, A., Willcox, K., Wang, Q.: Gradient-based multifidelity optimisation for aircraft design using Bayesian model calibration.  The Aeronautical Journal 115(1174), 729–738 (2011).  DOI 10.1017/S0001924000006473
  * (15) Marchildon, A.L., Zingg, D.W.: A Non-intrusive Solution to the Ill-Conditioning Problem of the Gradient-Enhanced Gaussian Covariance Matrix for Gaussian Processes.  Journal of Scientific Computing 95(3), 65 (2023).  DOI 10.1007/s10915-023-02190-w
  * (16) Mohammadi, H., Riche, R.L., Durrande, N., Touboul, E., Bay, X.: An analytic comparison of regularization methods for Gaussian Processes (2017)
  * (17) Ollar, J., Mortished, C., Jones, R., Sienz, J., Toropov, V.: Gradient based hyper-parameter optimisation for well conditioned kriging metamodels.  Structural and Multidisciplinary Optimization 55(6), 2029–2044 (2017).  DOI 10.1007/s00158-016-1626-8
  * (18) Osborne, M.A., Garnett, R., Roberts, S.J.: Gaussian Processes for Global Optimization.  In: 3rd International Conference on Learning and Intelligent Optimization. Trento, Italy (2009)
  * (19) Rasmussen, C.E., Williams, C.K.I.: Gaussian Processes for Machine Learning.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, Mass (2006)
  * (20) Schulz, E., Speekenbrink, M., Krause, A.: A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions.  Journal of Mathematical Psychology 85, 1–16 (2018).  DOI 10.1016/j.jmp.2018.03.001
  * (21) Shahriari, B., Swersky, K., Wang, Z., Adams, R.P., de Freitas, N.: Taking the Human Out of the Loop: A Review of Bayesian Optimization.  Proceedings of the IEEE 104(1), 148–175 (2016).  DOI 10.1109/JPROC.2015.2494218
  * (22) Toal, D.J., Bressloff, N.W., Keane, A.J., Holden, C.M.: The development of a hybridized particle swarm for kriging hyperparameter tuning.  Engineering Optimization 43(6), 675–699 (2011).  DOI 10.1080/0305215X.2010.508524
  * (23) Toal, D.J.J., Bressloff, N.W., Keane, A.J.: Kriging Hyperparameter Tuning Strategies.  AIAA Journal 46(5), 1240–1252 (2008).  DOI 10.2514/1.34822
  * (24) Toal, D.J.J., Forrester, A.I.J., Bressloff, N.W., Keane, A.J., Holden, C.: An adjoint for likelihood maximization.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 465(2111), 3267–3287 (2009).  DOI 10.1098/rspa.2009.0096
  * (25) Ulaganathan, S., Couckuyt, I., Dhaene, T., Degroote, J., Laermans, E.: Performance study of gradient-enhanced Kriging.  Engineering with Computers 32(1), 15–34 (2016).  DOI 10.1007/s00366-015-0397-y
  * (26) Won, J.H., Kim, S.J.: Maximum Likelihood Covariance Estimation with a Condition Number Constraint.  In: 2006 Fortieth Asilomar Conference on Signals, Systems and Computers, pp. 1445–1449. IEEE, Pacific Grove, CA, USA (2006).  DOI 10.1109/ACSSC.2006.354997
  * (27) Wu, A., Aoi, M.C., Pillow, J.W.: Exploiting gradients and Hessians in Bayesian optimization and Bayesian quadrature.  arXiv:1704.00060 [stat] (2018)
  * (28) Wu, J., Poloczek, M., Wilson, A.G., Frazier, P.: Bayesian Optimization with Gradients.  In: 31st Conference on Neural Information Processing Systems. Long Beach, CA, USA (2017)
  * (29) Zimmermann, R.: On the Maximum Likelihood Training of Gradient-Enhanced Spatial Gaussian Processes.  SIAM Journal on Scientific Computing 35(6), A2554–A2574 (2013).  DOI 10.1137/13092229X
  * (30) Zimmermann, R.: On the condition number anomaly of Gaussian correlation matrices.  Linear Algebra and its Applications 466, 512–526 (2015).  DOI 10.1016/j.laa.2014.10.038
  * (31) Zingg, D.W., Nemec, M., Pulliam, T.H.: A comparative evaluation of genetic and gradient-based algorithms applied to aerodynamic optimization.  European Journal of Computational Mechanics 17(1-2), 103–126 (2008).  DOI 10.3166/remn.17.103-126