Additive Kernels for Gaussian Process Modeling

keyword

Kriging, Computer Experiment, Additive Models, GAM, Maximum Likelihood Estimation, Relaxed Optimization, Sensitivity Analysis

Introduction

The study of numerical simulators often deals with calculation intensive computer codes. This cost implies that the number of evaluations of the numerical simulator is limited and thus many methods such as uncertainty propagation, sensitivity analysis, or global optimization are unaffordable. A well known approach to circumvent time limitations is to replace the numerical simulator by a mathematical approximation called metamodel (or response surface or surrogate model) based on the responses of the simulator for a limited number of inputs called the Design of Experiments (DoE). There is a large number of metamodels types and among the most popular we can cite regression, splines, neural networks… In this article, we focus on a particular type of metamodel: the Kriging method, more recently referred to as Gaussian Process modeling . Originally presented in spatial statistics as an optimal Linear Unbiased Predictor (LUP) of random processes, Kriging has become very popular in machine learning, where its interpretation is usually restricted to the convenient framework of Gaussian Processes (GP). Beyond the LUP —which then elegantly coincides with a conditional expectation—, the latter GP interpretation allows indeed the explicit derivation of conditional probability distributions for the response values at any point or set of points in the input space.

The first part of this paper focuses on the case of additive Gaussian Processes, their associated kernels and the properties of associated additive kriging models. The second part deals with a Relaxed Likelihood Maximization (RLM) procedure for the estimation of kernel parameters for Additive Kriging models. Finally, the proposed algorithm is compared with existing methods on a well known test function: the Sobol’s g-function . It is shown within the latter example that Additive Kriging with RLM outperforms standard Kriging and produce similar performances as GAM. Due to its approximation performance and its built-in probabilistic framework both demonstrated later in this article, the proposed Additive Kriging model appears as a serious and promising challenger among additive models.

Towards Additive Kriging

The proof of this property is given in appendix for $d=2$. For $d=n$ the proof follows the same pattern but the notations are more cumbersome. Note that the class of additive processes is not actually limited to processes with additive kernels. For example, let us consider $Z_{1}$ and $Z_{2}$ two correlated Gaussian processes on $(\Omega,\mathcal{F},P)$ such that the couple $(Z_{1},Z_{2})$ is Gaussian. Then $Z_{1}(x_{1})+Z_{2}(x_{2})$ is also a Gaussian process with additive paths but its kernel is not additive. However, in the next section, the term additive process will always refer to GP with additive kernels.

2 Invertibility of covariance matrices

A first approach is to remove some points in order to avoid any linear combination, which is furthermore in accordance with the aim of parsimonious evaluations for costly simulators. Algebraic methods may be used for determining the combination of points leading to a linear relationship between the values of the random process but this procedure is out of the scope of this paper.

3 Additive Kriging

Equations 3 are valid for any s.p.d kernel, thus they can be applied with additive kernels. In this case, the additivity of the kernel implies the additivity of the kriging mean: for example in dimension 2, for $K(\mathbf{x},\mathbf{y})=K_{1}(x_{1},y_{1})+K_{2}(x_{2},y_{2})$ we have

Another interesting property concerns the variance: $v$ can be null at points that do not belong to the DoE. Let us consider a two dimensional example where the DoE is composed of the 3 points represented on the left pannel of figure 1: $\mathcal{X}=\{\mathbf{x}^{(1)}\ \mathbf{x}^{(2)}\ \mathbf{x}^{(3)}\}$. Direct calculations presented in Appendix B shows that the prediction variance at the point $\mathbf{x}^{(4)}$ is equal to 0. This particularity follows from the fact that the value of the additive process are known almost surely at the point $x^{(4)}$ based on the observations at $\mathcal{X}$. In the next section, we illustrate the potential of Additive Kriging on an example and propose an algorithm for parameter estimation.

4 Illustration and further consideration on a 2D example

We present here a first basic example of an additive kriging model. We consider $D=^{2}$, and a set of 5 points in $D$ where the value of the observations $\mathbf{F}$ are arbitrarily chosen. Figure 2 shows the obtained kriging model. We can see on this figure the properties we mentioned above: the kriging mean is an additive function and the prediction variance can be null for points that do no belong to the DoE.

The effect of any variable can be isolated and represented so as the metamodel can be split in a sum of univariate sub-models. Moreover, we can get confidence intervals for each univariate model. As the expression of the first univariate model is

the effect of the direction 2 can be seen as an observation noise. We thus get an expression for the prediction variance of the first main effect

The expression of $m_{i}^{*}(x_{i})$ is straightforward whereas $v_{i}^{*}(x_{i})$ requires more calculations which are given in Appendix C.

The benefits of using $m_{i}^{*}$ and $v_{i}^{*}$ and then to define the sub-models up to a constant can be seen on the right panel of figure 3. At the end, the probabilistic framework gives an insight on the error of the metamodel but also of each sub-model.

Parameter estimation

Obtaining the optimal parameters $\psi_{i}^{*}$ relies on the succesful use of a non-convex global optimization routine. This can be severely hindered for large values of $d$ since the search space of kernel parameters becomes high dimensional. One way to cope with this issue is to separate the variables and split the optimization into several low-dimensional subproblems.

2 The Relaxed Likelihood Maximization algorithm

The aim of the Relaxed Likelihood Maximization (RLM) algorithm is to treat separately the optimization in each direction. In this way, RLM can be seen as a cyclic relaxation optimization procedure with initial values of the parameters $\sigma^{2}_{i}$ set to zero. As we will see, the main originality here is to consider a kriging model with an observation noise variance $\tau^{2}$ that fluctuates during the optimization. This parameter account for the metamodel error (if the function is not additive for example) but also for the inaccuracy of the intermediate values of $\sigma_{i}$ and $\theta_{i}$.

The first step of the algorithm is to estimate the parameters of the kernel $K_{1}$. The simplification of the method is to consider that all the variations of $Y$ in the other directions can be summed up as a white noise. Under this hypothesis, $l$ depends on $\psi_{1}$ and $\tau$:

Then, the couple $\{\psi_{1}^{*},\tau^{*}\}$ that maximizes $l(\psi_{1},\tau)$ can be obtained by numerical optimization.

The second step of the algorithm consists in estimating $\psi_{2}$, with $\psi_{1}$ fixed to $\psi_{1}^{*}$:

This operation can be repeated for all the directions until the estimation of $\psi_{d}$. However, even if all the parameters $\psi_{i}$ have been estimated, it is fruitful to re-estimate them such that the estimation of the parameter $\psi_{i}$ can benefit of the values $\psi_{j}^{*}$ for $j>i$. Thus, the algorithm is composed of a cycle of estimations that treat each direction one after each other:

$\tau$ is a parameter tuning the fidelity of the metamodel since for $\tau=0$ the kriging mean interpolates the data. In practice, this parameter is decreasing at almost each new estimation. Depending on the observations and on the DoE, $\tau$ converges either to a constant or to zero (cf. the g-function example and figure 6). When zero is not reached, $\tau^{2}$ should correspond to the part of the variance that cannot be explained by the additive model. Thus, the comparison between $\tau^{2}$ and the $\sigma_{i}^{2}$ allows us to quantify the degree of additivity of the objective function according to the metamodel.

This procedure of estimation is not meant to be applied for kernels that are not additive. The method developed by Welch for tensor product kernels in has similarities since it corresponds to a sequential estimation of the parameters. One interesting feature of Welch’s algorithm is to choose at each step the best search direction for the parameters. The RLM algorithm could easily be adapted in a similar way to improve the quality of the results but the corresponding adapted version would be much more time consuming.

Comparison between the optimization’s methods

The aim of this section is to compare the RLM algorithm to the Usual Likelihood Maximization (ULM). The test functions that are considered are paths of an additive GP $Y$ with Gaussian additive kernel $K$. For this example, the parameters of $K$ are fixed to $\sigma_{i}=1$, $\theta_{i}=0.2$ for $i\in 1\dots d$ but those values are supposed to be unknown.

Here, $2\times d+1$ parameters have to be estimated: $d$ for the variances, $d$ for the range and 1 for the noise variance $\tau^{2}$. For ULM, they are estimated simultaneously, whereas the RLM is a 3-dimensional optimization at each step. In both cases, we use the L-BFGS-B method of the function optim with the R software. To judge the effectiveness of the algorithms, we compare here the best value found for the log-likelihood $l$ to the computational budget (the number of call to $l$) required for the optimization. As the optim function returns the number $nc$ of call to $l$ and the best value $bv$ at the end of each optimization, we obtain for the MLE on one path of $Y$ one value of $nc$ and $bv$ for ULM and $nb\_iteration\times d$ values of $nv$ and $bv$ for RLM since there is one optimization at each step of each iteration.

The panel (a) of figure 4 presents the results for the two optimizations on a path of a GP for $d=5$. On this example, we can see that ULM needs 1500 calls to the log-likelihood before convergence whereas RLM requires much more calls before convergence. However, the result of the two methods are similar for 1500 calls but the result of RLM after 5000 calls is substantially improved. In order to get more robust results we simulate 20 paths of $Y$ and we observe the global distribution of the variables $nv$ and $bv$. Furthermore, we study the evolution of the algorithm performances when the dimension increases choosing various values for the parameter $d$ from 3 to 18 with a Latin Hypercube (LH) Design with maximin criteria containing $10\times d$ points. We observe on the panels (b), (c) and (d) of figure 4 that optimization with the RLM requires more calls to the function $l$, but this method appears to be more efficient and robust than ULM. Those results are stressed by figure 5 where the final best value of RLM and ULM are compared. This figure also shows that the advantage of using RLM comes bigger when $d$ is getting larger.

Application to the g-function of Sobol

In order to illustrate the methodology and to compare it to existing algorithms, an analytical test case is considered. The function to approximate is the g-function of Sobol defined over $^{d}$ by

This popular function in the literature is obviously not additive. However, depending on the coefficients $a_{k}$, $g$ can be very close to an additive function. As a rule, the g-function is all the more additive as the $a_{k}$ are large. One main advantage for our study is that the Sobol sensitivity indices can be obtained analytically so we can quantify the degree of additivity of the test function. For $i=1,\dots,d$ the indice $S_{i}$ associated to the variables $x_{i}$ is

Here we limit ourselves to the case $d=4$ and following we choose $a_{k}=k$ for $k\in\{1,...,4\}$. For this combination of parameters, the sum of the first order Sobol indices is 0.95 so the g-function is almost additive. The considered DoE are LH maximin designs based on 40 points. To asses the quality of the obtained metamodels, the predictivity coefficient $Q_{2}$ is computed on a test sample of $n=1000$ points uniformly distributed over $^{4}$. Its expression is:

where $\mathbf{y}$ is the vector of the values at the test points, $\hat{\mathbf{y}}$ is the vector of predicted values and $\bar{\mathbf{y}}$ is the mean of $\mathbf{y}$.

We run on this example 5 iterations of the RLM algorithm with kernel Matèrn 3/2. The evolution of the estimated observation noise $\tau^{2}$ is represented on figure 6. On this figure, it appears that the observation noise is decreasing as the estimation of the parameters is improved. Here, the convergence of the algorithm is reached at iteration 4. The overall quality of the constructed metamodel is high since $Q_{2}=0.91$ and the final value for $\tau^{2}$ is 0.01.

As previously the expression of the univariate sub-metamodels is

The univariate functions obtained are presented on figure 7. The confidence intervals are not represented here in order to enhance the readability of the graphics and the represented values are centered to ensure that the observations and the univariate functions are comparable.

As the value of $Q_{2}$ is likely to fluctuate with the DoE and the optimization performances, we compare here the proposed RLM algorithm with other methods for 20 different LHS. The other methods used for the test are (a) additive kriging model with ULM, (b) kriging with usual tensor-product kernel, (c) the GAM algorithm. The results for classical kriging and GAM are obtained with the DiceKrigingAs for RLM and ULM, DiceKriging also use the BFGS algorithm for the likelihood maximization and the GAM packages for R available on the CRAN . As the value of the $a_{k}$ are the same as in where Marrel et al. presents a specific algorithm for sequential parameter estimation in non-additive kriging models, the results of this paper are also presented as method (d). The mean and the standard deviation of the obtained $Q_{2}$ are gathered in table 1.

Concluding remarks

The proposed methodology seems to be a good challenger for additive modeling. On the example of the GP paths, the RLM appears to be more efficient than usual likelihood maximization and well suited for high dimensional modeling. On the second example, additive models benefit of the important additive component of the g-function and outperform non additive models even if the function is not purely additive. The predictivity of the RLM is equivalent to that of GAM but its robustness is higher for this example.

One main difference between RLM and GAM backfitting is that RLM takes into account the estimated parameters into the covariance structure whereas GAM subtracts from the observation the predicted value for all the sub-models obtained in the other directions.

At the end, the proposed methodology is fully compatible with Kriging-based methods and its versatile applications. For example, one can choose a well suited kernel for the function to approximate or use additive kriging for high-dimensional optimization strategies relying on the expecting improvement criteria.

References

Appendix A: Proof of proposition 1 for d=2𝑑2d=2

Appendix B: Calculation of the prediction variance

Let consider a DoE composed of the 3 points $\{\mathbf{x}^{(1)}\ \mathbf{x}^{(2)}\ \mathbf{x}^{(3)}\}$ represented on the left pannel of figure 1. We want here to show that although $\mathbf{x}^{(4)})$ does not belongs to the DoE we have $v(\mathbf{x}^{(4)})=0$.