UQL AB U SER M ANUAL

K RIGING (G AUSSIAN P ROCESS M ODELING )

C. Lataniotis, D. Wicaksono, S. Marelli, B. Sudret

C HAIR OF R ISK , S AFETY AND U NCERTAINTY Q UANTIFICATION Risk, Safety &

S TEFANO -F RANSCINI -P LATZ 5 Uncertainty Quantification

CH-8093 Z ÜRICH
How to cite UQL AB

S. Marelli, and B. Sudret, UQLab: A framework for uncertainty quantification in Matlab, Proc. 2nd Int. Conf. on
Vulnerability, Risk Analysis and Management (ICVRAM2014), Liverpool, United Kingdom, 2014, 2554-2563.

How to cite this manual

C. Lataniotis, D. Wicaksono, S. Marelli, B. Sudret, UQLab user manual – Kriging (Gaussian process modeling),
Report UQLab-V2.1-105, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2024

BibTeX entry

@TechReport{UQdoc_21_105,
author = {Lataniotis, C. and Wicaksono, D. and Marelli, S. and Sudret, B.},
title = {{UQLab user manual -- Kriging (Gaussian process modeling)}},
institution = {Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich,
Switzerland},
year = {2024},
note = {Report UQLab-V2.1-105}
}
Document Data Sheet

Document Ref. UQL AB-V2.1-105

Title: UQL AB user manual – Kriging (Gaussian process modeling)

Authors: C. Lataniotis, D. Wicaksono, S. Marelli, B. Sudret

Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich,
Switzerland
Date: 15/04/2024

Doc. Version Date Comments

V2.1 15/04/2024 UQL AB V2.1 release
• Added GP resampling.
V2.0 01/02/2022 UQL AB V2.0 release
V1.4 01/02/2021 UQL AB V1.4 release
V1.3 19/09/2019 UQL AB V1.3 release
• Added a new introduction section
• Added documentation related to the Gaussian process re-
gression (GPR) feature.
• Updated figures and default options
V1.2 22/02/2019 UQL AB V1.2 release
• Added an automatic calculation of validation error if a val-
idation set is provided
V1.1 05/07/2018 UQL AB V1.1 release
• Minor updates in the reference list
V1.0 01/05/2017 UQL AB V1.0 release
• Minor consistency fixes in the reference list
• New usage sections regarding scaling and sampling from a
Gaussian process posterior
V0.9 01/07/2015 Initial release
 Abstract

Kriging is a stochastic modeling algorithm that has many applications in most fields of en-
gineering and applied mathematics. The UQLab metamodeling tool provides an efficient,
modular, and easy to use Kriging module that allows one to obtain efficient Kriging pre-
dictors with minimal effort. For experienced users, numerous more advanced configuration
options can be easily set.

This guide is an extensive manual of the Kriging metamodeling module and divided into
three parts:

• A short introduction to the main concepts and techniques behind Kriging as interpola-
tion and regression models; with a selection of relevant references in the literature

• A detailed example-based user guide, with explanations of most of the available options
and methods;

• A comprehensive reference list of all the available options and functionalities of the
module.

Keywords: UQL AB, Kriging, Gaussian process metamodeling, Gaussian process regression
Contents

1 Theory 1

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Kriging basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1 Prediction with noise-free responses (interpolation) . . . . . . . . . . . 2

1.2.2 Prediction with noisy responses (regression) . . . . . . . . . . . . . . . 3

1.2.3 Kriging metamodel ingredients . . . . . . . . . . . . . . . . . . . . . . 5

1.3 Trend types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.3.1 Commonly used trends . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4 Correlation functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.4.1 Correlation families . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1.4.2 Correlation function types . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.4.3 Isotropic correlation functions . . . . . . . . . . . . . . . . . . . . . . . 11

1.4.4 Nugget and numerical stability . . . . . . . . . . . . . . . . . . . . . . 12

1.5 Estimation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.5.1 Maximum-likelihood estimation . . . . . . . . . . . . . . . . . . . . . . 12

1.5.2 Cross-validation estimation . . . . . . . . . . . . . . . . . . . . . . . . 14

1.6 Optimization methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

1.7 A posteriori error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.7.1 Leave-one-out cross-validation error . . . . . . . . . . . . . . . . . . . 17

1.7.2 Validation error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.8 Gaussian process trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2 Usage 20

2.1 Reference problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

 2.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.3 Kriging metamodel calculation: noise-free case . . . . . . . . . . . . . . . . . 21

2.3.1 Accessing the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.4 Kriging metamodel setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.4.1 Specification and generation of experimental design . . . . . . . . . . 25

2.4.2 Trend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4.3 Correlation function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4.4 Estimation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.4.5 Optimization methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

2.5 Kriging metamodels with noise (regression) . . . . . . . . . . . . . . . . . . . 35

2.5.1 Unknown homoscedastic noise . . . . . . . . . . . . . . . . . . . . . . 35

2.5.2 Known noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

2.6 Kriging metamodels of vector-valued models . . . . . . . . . . . . . . . . . . . 38

2.6.1 Accessing the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.7 Using a validation set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

2.8 Using Kriging predictor as a model . . . . . . . . . . . . . . . . . . . . . . . . 39

2.9 Specifying manually a Kriging predictor (predictor-only mode) . . . . . . . . . 40

2.10 Performing Kriging on an auxiliary space (scaling) . . . . . . . . . . . . . . . 41

2.11 Drawing sample paths from a Gaussian process posterior . . . . . . . . . . . . 42

2.12 Using Kriging with constant inputs . . . . . . . . . . . . . . . . . . . . . . . . 43

3 Reference List 45

3.1 Create a Kriging metamodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.1.1 Experimental design options . . . . . . . . . . . . . . . . . . . . . . . . 49

3.1.2 Trend options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.1.3 Correlation function options . . . . . . . . . . . . . . . . . . . . . . . . 51

3.1.4 Estimation method options . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1.5 Hyperparameters optimization options . . . . . . . . . . . . . . . . . . 52

3.1.6 Regression options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.1.7 Custom Kriging options . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.1.8 Validation Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.2 Accessing the Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3.2.1 Internal fields (advanced) . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.3 Kriging predictor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

3.4 Printing and visualizing a Kriging metamodel . . . . . . . . . . . . . . . . . . 64

3.4.1 Printing the results: uq_print . . . . . . . . . . . . . . . . . . . . . . 64

3.4.2 Visualize the results: uq_display . . . . . . . . . . . . . . . . . . . . 66

References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Chapter 1

Theory

1.1 Introduction

In modern applied sciences and engineering, computational models have become more and
more complex and hence expensive to evaluate. In this context, metamodeling (or surro-
gate modeling) can reduce the associated computational costs by substituting an expensive
computational model with a metamodel, a functional approximation of the original model
that is much faster to evaluate. It is constructed by learning the approximation from a rela-
tively small set of input parameters and their corresponding model responses generated from
running the original expensive model. Because it is faster to evaluate, a metamodel allows
for more sophisticated analyses, including, e.g., sensitivity and reliability analysis, Bayesian
model calibration, etc.

In its original form, Kriging (also known as Gaussian process modeling) is a statistical interpo-
lation method that capitalizes on Gaussian processes to interpolate a wide range of complex
functions. It was first developed as a spatial interpolation tool in geostatistics by Krige dating
back to the 1950s (Krige, 1951) and formalized by Matheron in the 1960s (Matheron, 1963).
Kriging was later introduced in the context of metamodeling and computer experiments in
the work of Sacks et al. (1989) in which Kriging was used to represent an input/output
mapping of an expensive computational model. In the mid-2000s, Gaussian processes en-
joyed renewed interest due to their applications in machine learning regression and classi-
fication, thanks to the introduction of GP-regression, which supports noisy data(Rasmussen
and Williams, 2006). The reader is referred to Santner et al. (2003) for an in-depth introduc-
tion of Kriging as a metamodeling tool and to Rasmussen and Williams (2006) for a similar
introduction to Kriging as a regression and classification tool.

This user manual presents the Kriging metamodeling tool of UQL AB that supports Kriging
for both interpolation and regression. This chapter briefly presents the basics of Kriging,
including its formulation, main ingredients, as well as the estimation of its parameters from
data. Section 2 then details the usage of Kriging within UQlab. Finally, Section 3 provides a
comprehensive list of the available commands and options.

1
UQL AB user manual

1.2 Kriging basics

Kriging (or Gaussian process modeling) is a stochastic algorithm which assumes that the
model output M(x) is a realization of a Gaussian process indexed by x ∈ DX ⊂ RM . A
Kriging metamodel is described by the following equation (Santner et al., 2003)

MK (x) = β T f (x) + σ 2 Z (x, ω) . (1.1)

The first term in Eq. (1.1), β T f (x), is the mean value of the Gaussian process (i.e., its trend);
it consists of P arbitrary functions {fj ; j = 1, . . . , P } and the corresponding coefficients
{βj ; j = 1, . . . , P }. The second term in Eq. (1.1) consists of the (constant) variance of the
Gaussian process σ 2 and a zero-mean, unit-variance, stationary Gaussian process Z(x, ω).
The underlying probability space is represented by ω and is defined in terms of a correlation
function R (i.e., correlation family) and its hyperparameters θ. The correlation function
R = R(x, x0 ; θ), in turn, describes the correlation between two sample points in the output
space, that depends on x, x0 , and the hyperparameters θ.

1.2.1 Prediction with noise-free responses (interpolation)

In the context of metamodeling, it is of interest to predict MK (x) for a new point x, given
the (observed) experimental design X = x(1) , . . . , x(N ) and the corresponding noise-free


model responses Y = y (1) = M x(1) , . . . , y (N ) = M x(N ) . A Kriging metamodel (or



Kriging predictor) provides the prediction based on the Gaussian properties of the process.

The Gaussian assumption states that the vector formed by the prediction at x, i.e., Yb (x) and
the true model responses Y, has a joint Gaussian distribution defined by
 T
r T (x)
    
Yb (x) f (x) β 2 1
∼ NN +1 ,σ (1.2)
Y Fβ r (x) R

where:

• F is the observation (design) matrix of the Kriging metamodel trend. It reads

 
Fij = fj x(i) , i = 1, . . . , N ; j = 0, . . . , P. (1.3)

• r(x) is the vector of cross-correlations between the prediction point x and each one of
the observations whose elements read as
 
ri = R x, x(i) ; θ , i = 1, . . . , N. (1.4)

• R is the correlation matrix with elements:

 
Rij = R x(i) , x(j) ; θ , i, j = 1, . . . , N. (1.5)

UQL AB-V2.1-105 -2-

 Kriging (Gaussian process modeling)

Consequently, the mean and variance of the Gaussian random variate Yb (x) conditional on
the observed data X and Y can be computed as follows (Santner et al., 2003; Dubourg, 2011)
 
µYb (x) = f (x)T βb + r (x)T R−1 Y − F βb , (1.6)

−1 
σY2b (x) = σ 2 1 − r T (x) R−1 r (x) + uT (x) F T R−1 F u (x) , (1.7)

where

−1
βb = F T R−1 F F T R−1 Y (1.8)

is the generalized least-squares estimate of β and

u (x) = F T R−1 r (x) − f (x) . (1.9)

Eqs. (1.6) and (1.7) are referred to as the mean and variance of the Kriging predictor, respec-
tively. An important property of this particular predictor is that the variance of the prediction
at an experimental design point x ∈ X collapses to zero. In other words, the Kriging predictor
is an interpolant with respect to the experimental design points.

Another useful corollary of the Gaussian assumption is that

 
Yb (x) ∼ N µYb (x) , σY2b (x) (1.10)

and therefore
t − µYb (x)
   
P Y (x) ≤ t = Φ
b , (1.11)
σYb (x)
where Φ(·) denotes the Gaussian cumulative density function. Note that this is the probability
that Yb (x) is smaller than t at a given x1 . Based on Eq. (1.11), the confidence intervals on the
predictor can be calculated by
h  α  α i
Yb (x) ∈ µYb (x) − Φ−1 1 − σYb (x) , µYb (x) + Φ−1 1 − σYb (x) (1.12)
2 2
with probability 1 − α.

1.2.2 Prediction with noisy responses (regression)

There are cases in which the observed responses are noisy, that is

y = M(x) + ε, (1.13)

where it is often assumed that the additive noise ε follows a zero-mean Gaussian distribution

ε ∼ N (0, Σn ) , (1.14)

where Σn is the covariance matrix of the noise term.

1
This shall not be confused with a probability of failure computed in the framework of rare event simulation.
For details, see UQL AB User Manual – Structural Reliability.

UQL AB-V2.1-105 -3-

UQL AB user manual

Depending on the properties of Σn , three classes of noise can be identified:

• Σn = σn2 I (where I is an identity matrix) corresponds to homogeneous (homoscedastic)

noise, in which the variance σn2 is constant and the same for all observed responses. In
other words, the noise is an independent and identically distributed (iid) Gaussian.

• Σn = diag σn2 corresponds to independent heterogeneous (heteroscedastic) noise, in

which the noise variances σn2 differ for each observed response but are not correlated.

• Σn corresponds to general heteroscedastic noise, in which the noise variance can differ
for each observed response and correlations between observations are possible.

The joint Gaussian distribution formed by the prediction at x, i.e., Yb (x) and the observed,
albeit noisy, model responses Y now becomes:
 T
σ2 σ 2 r T (x)
    
Yb (x) f (x) β
∼ NN +1 , . (1.15)
Y Fβ σ 2 r (x) σ 2 R + Σn

Because in general the noise variance cannot be factored out from the covariance matrix (i.e.,
due to heteroscedasticity), the covariance matrix C = σ 2 R + Σn is used in the subsequent
formulation. Consequently, the Kriging mean and variance predictors (see Eq (1.6) and (1.7))
become  
µYb (x) = f (x)T βb + c(x)T C −1 Y − F βb , (1.16)

σY2b (x) = σ 2 − cT (x) C −1 c(x) + uTc (x)(F T C −1 F )−1 uc (x)

(1.17)

where c = σ 2 r (x) is the cross-covariance vector,

−1
βb = F T C −1 F F T C −1 Y, (1.18)

is the generalized least-square estimate of the regression coefficients β,

b and

uc (x) = F T C −1 c(x) − f (x). (1.19)

In the special case of homoscedastic noise, C = σ 2 R + σn2 I, thus with

2
σtotal = σ 2 + σn2 (1.20)

and
σn2
τ= 2 (1.21)
σtotal
Eq. (1.15) can be written as

(1 − τ ) reT (x)
   T   
Yb (x) f (x) β 2
∼ NN +1 , σtotal (1.22)
Y Fβ re (x) Re

where
re = (1 − τ ) r (1.23)

UQL AB-V2.1-105 -4-

 Kriging (Gaussian process modeling)

and
e = (1 − τ ) R + τ I.
R (1.24)

In this case, the mean and variance of the Gaussian random variate Yb (x) can be computed
as follows (Rasmussen and Williams, 2006)
 
µYb (x) = f (x)T βb + re(x)T R
e −1 Y − F βb (1.25)
 
σY2b (x) = σtotal
2 e −1 re(x) + uT (x)(F T R
(1 − τ ) − reT (x) R e −1 F )−1 u(x) (1.26)

where βb and u in the above follow Eqs. (1.8) and (1.9), but with R and r are replaced by R
e
and re, respectively.

Unlike in the noise-free case in Section 1.2.1, the variance of the prediction at an experi-
mental design point x ∈ X does not collapse to zero and the Kriging predictor becomes a
regression model. Other corollaries following the Gaussian assumptions, however, still apply.

1.2.3 Kriging metamodel ingredients

In practice, in order to obtain a Kriging metamodel, the following steps are needed:

• Select a functional basis of the Kriging trend. This is further discussed in Section 1.3.

• Select a correlation function R(x, x0 ; θ). This is further discussed in Section 1.4.

• If the hyperparameters θ are unknown, they need to be estimated from the available
data. This involves setting up an optimization problem (see Section 1.5) and solving it
(see Section 1.6).

• In the case of regression, then either the noise variance σn2 (if it is unknown) or the
Gaussian process variance σ 2 (if σn2 is known) need to be estimated.

• Using the optimal value of θ, calculate the rest of the unknown Kriging parameters
(e.g., β) and, if necessary, σ 2 .

Predictions at new points can then be made in terms of the mean and variance of Yb (x).

Two examples of one-dimensional Kriging metamodels are shown in Figure 1. They are
constructed using two sets of experimental design and the corresponding responses; both sets
come from the same underlying model. In the case of noise-free responses (Figure 1(a)), the
Kriging predictor returns the mean and the variance of a Gaussian process that interpolates
the design points. In the case of noisy responses (Figure 1(b)), the predictor serves as a
regression model that fits the noisy responses. Using Eq. (1.12), the 95% confidence intervals
of both predictors are also calculated and shown in the plot.

UQL AB-V2.1-105 -5-

UQL AB user manual

(a) noise-free responses (interpolation) (b) noisy responses (regression)

Figure 1: Examples of one-dimensional Kriging metamodels with two different responses.

1.3 Trend types

The trend refers to the mean of a Kriging metamodel, that is, the term β T f (x) in Eq. (1.1).
While using a non-zero trend is optional, it is commonly preferred2 . In the literature, a
different naming is given to the Kriging metamodel depending on the type of trend that is
used (Dubourg, 2011), namely:

• Simple Kriging
In simple Kriging, the trend is known
P
X
β T f (x) = fj (x),
j=1

where fj ’s are arbitrary but fully specified functions. Note that no estimation on β is
taken place and the coefficients β are all 1’s
• Ordinary Kriging
In ordinary Kriging, the trend has a constant yet unknown value

β T f (x) = β0 f0 (x) = β0

where by convention f0 (x) = 1.

• Universal Kriging
Universal Kriging, the most general and flexible formulation, assumes that the trend is

2
Note that the mean of the Kriging predictor given in Eq. (1.6), (1.16), or (1.25) is not confined to be zero
even though the trend is zero.

UQL AB-V2.1-105 -6-

 Kriging (Gaussian process modeling)

a linear combination of prescribed arbitrary functions (e.g., polynomials)

P
X
T
β f (x) = βj fj (x).
j=0

According to the above, the simple and ordinary Kriging are special cases of universal Kriging.
In the current version of UQL AB, the functions fj ’s can either be constants, polynomials, or
any other user-defined functions.

1.3.1 Commonly used trends

The most commonly used trends are based on polynomial basis and summarized in Table 1.

Table 1: Trend options combinations

Trend Formula
constant (ordinary Kriging) β0
PM
linear β0 + i=1 βi xi
PM
quadratic β0 + i=1 βi xi + ΣM M
i=1 Σj=1 βij xi xj

Based on the Table 1, a multivariate polynomial trends can be written general form as follows

X M
Y
T
β f (x) = βα fα (x) , fα (x) = xαi i (1.27)
α∈AM,P i=1

where α = {α1 , . . . , αM } is a vector of indices and AM,P = α ∈ NM : |α| ≤ P denotes the



set of indices that corresponds to all polynomials in the M input variables up to degree P .

The Kriging module in UQL AB offers the possibility to use a multivariate polynomial of arbi-
trary degree P (Eq. (1.27)) as well as different types of basis functions. See Section 2.4.2 for
more information and Section 3.1.2 for a list of all the available options regarding the trend.

1.4 Correlation functions

The correlation function (also called the kernel or covariance3 function in the literature) is a
crucial ingredient for a Kriging predictor, since it contains the assumptions about the approx-
imation function. Generally speaking, it describes the “similarity” between observations and
new points, i.e., how similar such points are, depending on the distance between the input
points. In this sense, input points that are close to each other are expected to have similar
outputs or responses.

An arbitrary function of (x, x0 ) is, in general, not a valid correlation function. The necessary
conditions for a function R (x, x0 ) to be a correlation function are:

3
The covariance function is to the correlation function multiplied by the Gaussian process variance σ 2 .

UQL AB-V2.1-105 -7-

UQL AB user manual

• The correlation matrix4 whose elements are computed as Rij = R x(i) , x(j) , (xi , xj ) ∈

X × X is positive semi-definite for any choice of number of sample points N in experi-

mental design X .

• The function is symmetric, i.e., R(x, x0 ) = R(x0 , x), ∀x0 , x ∈ DX .

In general, a correlation function can be expressed in the form R(x, x0 ; θ) where θ is a vector
that contains a set of parameters. Typically θ ∈ RM , yet this is not necessarily true in the
general case, since more than one parameter might correspond to each input dimension. For
notational clarity, however, it is assumed that one element of θ is used per input dimension
in the subsequent discussion.

1.4.1 Correlation families

All the correlation families described below are stationary correlation functions that depend
only on the relative position of its two inputs. Moreover, these families are one-dimensional
and defined for a pair of one-dimensional input x, x0 ∈ R and parametrized by θ ∈ R>0 , often
referred to as the characteristic length scale or simply the scale parameter.

For each of the available correlation families in UQL AB, the function is plotted in Figures 2-6
with different scale values, along with the corresponding sample paths (i.e., realizations) of
a zero-mean, unit-variance Gaussian process having this correlation function.

• Linear:
|x − x0 |
 
0


R x, x ; θ = max 0, 1 − . (1.28)
θ

Figure 2: The linear correlation function (Eq. (1.28)) and sample paths drawn from the
corresponding zero-mean unit-variance Gaussian process for various scale parameters.

• Exponential:
|x − x0 |
 
0


R x, x ; θ = exp − . (1.29)
θ
The resulting sample paths are C 0 , i.e., continuous but non-differentiable.
4
It is also a Gram matrix.

UQL AB-V2.1-105 -8-

 Kriging (Gaussian process modeling)

Figure 3: The exponential correlation function (Eq. (1.29)) and sample paths drawn from
the corresponding zero-mean unit-variance Gaussian process for various scale parameters.

• Gaussian (or squared exponential):

" 2 #
|x − x0 |

1
R x, x0 ; θ = exp −


. (1.30)
2 θ

The resulting sample paths of the corresponding process are infinitely differentiable.

Figure 4: The Gaussian correlation function (Eq. (1.30)) and sample paths of the correspond-
ing zero-mean unit-variance Gaussian process for various scale parameters.

• Matérn: The general form of the Matérn correlation function is given by

√ |x − x0 | v √ |x − x0 |
   
0

1
R x, x ; θ, v = v−1 2 v Kv 2 v (1.31)
2 Γ(v) θ θ

where θ is the correlation function scale parameter, v ≥ 1/2 is the shape parameter, Γ
is the Euler’s Gamma function, and Kν is the modified Bessel function of the second
kind5 .

An interesting feature of this correlation function family is that the sample paths of the
corresponding Gaussian process are dv − 1e times differentiable, where d·e denotes the
ceiling function. For v = 1/2, Matérn kernel coincides with the exponential correlation
5
It is also also known as the Bessel function of the third kind. For details, see Abramovitz and Stegun (1965).

UQL AB-V2.1-105 -9-

UQL AB user manual

function which generates C 0 sample paths. For v → ∞, it tends towards the Gaussian
correlation function which generates C ∞ sample paths.

The Matérn functions can be computed using simplified formulas when v = p + 1/2,
where p is a non-negative integer (Rasmussen and Williams, 2006). The Matérn func-
tions with v = 3/2 and 5/2 (abbreviated as Matérn-3/2 and Matérn-5/2, respectively)
are the most commonly used and available in UQL AB. Their formulas are given below.

For v = 3/2
√ |x − x0 | √ |x − x0 |
   
R (h; θ, v = 3/2) = 1 + 3 exp − 3 , (1.32)
θ θ

Figure 5: The Matérn 3/2 correlation function (Eq. (1.32)) and sample paths drawn from the
corresponding zero-mean unit-variance Gaussian process for various scale parameters.

and for v = 5/2

2 !
√ |x − x0 | 5 |x − x0 | √ |x − x0 |
  
0
R(x, x ; θ, v = 5/2) = 1+ 5 + exp − 5 . (1.33)
θ 3 θ θ

Figure 6: The Matérn 5/2 correlation function (Eq. (1.33)) and sample paths drawn from the
corresponding zero-mean unit-variance Gaussian process for various scale parameters.

UQL AB-V2.1-105 - 10 -
 Kriging (Gaussian process modeling)

1.4.2 Correlation function types

When the input dimension M is greater than one, multi-dimensional correlation function
can be constructed from one-dimensional correlation families using one of the following
constructions:

• Ellipsoidal correlation functions (Rasmussen and Williams, 2006), calculated as follows

"M  #0.5
X xi − x0 2
0 i


R x, x ; θ = R (h) , h = . (1.34)
θi
i=1

• Separable correlation functions (Sacks et al., 1989; Dubourg, 2011), calculated as fol-
lows
M
Y
0
R xi , x0i , θi .



R x, x ; θ = (1.35)
i=1

The function R(·) (resp, R(·, ·; ·)) that appears on the right-hand side of Eq. (1.34) (resp.
Eq. (1.35)) corresponds to one-dimensional correlation functions described in Section 1.4.1.

Note: When using one of the correlation families in Section 1.4.1 to build a multi-
0|
dimensional ellipsoidal correlation function, the term |x−x
θ inside the function is
replaced with h.

1.4.3 Isotropic correlation functions

The multi-dimensional correlation functions given in Section 1.4.2 above correspond to the
anisotropic case, in which there is a unique scale parameter for each input dimension. On the
other hand, a multi-dimensional correlation function is called isotropic when a single scale
parameter θ is associated with all the input dimensions. Generally speaking, a correlation
function is called isotropic when it has the same behavior over all dimensions.

In that sense, for each type of correlation function, the isotropy is defined as follows:

• Isotropic ellipsoidal correlation function given as

"M #0.5
0

1 X 0 2


R x, x ; θ = R (h) ; h = xi − xi . (1.36)
θ
i=1

• Isotropic separable correlation function given as

M
Y
0
R xi , x0i ; θ .



R x, x ; θ = (1.37)
i=1

UQL AB-V2.1-105 - 11 -
UQL AB user manual

1.4.4 Nugget and numerical stability

Regardless on how the correlation matrix R is calculated, it is often the case that it needs to
be inverted at various stages of the Kriging process (see, for example, Eq. (1.6) to (1.9)). This
inversion is well known to suffer from numerical instabilities, especially when the distances
between the design points xi are small and the responses are noise-free. To circumvent this
limitation, a nugget ν can be introduced. Nugget is a set of values that are added to the main
diagonal of R, such that
Rii = 1 + νi . (1.38)

1.5 Estimation methods

To obtain a Kriging metamodel of Eq. (1.1), usually the hyperparameters θ are unknown and
thus must be estimated. The estimation is achieved by solving an optimization problem that
differs depending on the selected estimation method, and whether the responses are noisy.
The available estimation methods in UQL AB are discussed in the following subsections.

1.5.1 Maximum-likelihood estimation

The idea behind the maximum-likelihood (ML) estimation method is to find the set of Krig-
ing parameters β, σ 2 , θ, and if apply, σn2 , such that the likelihood of the observations Y =

T
M x(1) , . . . , M x(N )


is maximized. Since Y is assumed to follow a multivariate Gaus-
sian distribution (recall basic Kriging assumptions), the likelihood function L β, σ 2 , θ; Y

reads
(det C)−1/2
 
2

1 T −1
L β, σ , θ; Y = exp − (Y − F β) C (Y − F β) (1.39)
(2π)N/2 2
where the covariance matrix C sums up the covariance matrix of the Gaussian process σ 2 R
and covariance matrix of the noise of the responses Σn as follows

C = σ 2 R + Σn . (1.40)

Depending on whether such noise is present as well as the nature of the noise, different
optimization problems can be set up as described in the sequel.

1.5.1.1 Noise-free Kriging

For models with noise-free responses, the covariance matrix C reduces to σ 2 R, and the
likelihood function can be written as

(det R)−1/2
 
2

1 T −1
L β, σ , θ; Y = exp − 2 (Y − F β) R (Y − F β) . (1.41)
(2πσ 2 )N/2 2σ

By maximizing the quantity described in Eq. (1.41), the following analytical estimates of

UQL AB-V2.1-105 - 12 -
 Kriging (Gaussian process modeling)

β and σ 2 that are strictly functions of θ are obtained (for proof and more details, refer to
Dubourg (2011); Santner et al. (2003))

−1
βb = β (θ) = F T R−1 F F T R−1 Y, (1.42)

1
b2 = σ 2 (θ) =
σ (Y − F β)T R−1 (Y − F β) . (1.43)
N
The hyperparameters θ, in turn, are obtained from solving the following optimization prob-
lem6
θb = arg min [− log L (θ; Y)] . (1.44)
θ∈Dθ

Based on Eqs. (1.41) to (1.43), the optimization problem of Eq. (1.44) can be written as
1
log (det R) + N log 2πσ 2 + N .


θb = arg min (1.45)
θ∈Dθ 2

1.5.1.2 Kriging with unknown homogeneous noise (homoscedastic)

In the case of noisy responses with an unknown homogeneous (homoscedastic) noise variance,
the covariance matrix C reduces to σ 2 R, e where σ 2 = σ 2 + σ 2 , R
total
e = (1 − τ ) R + τ I, and
total n
σn2
τ= 2
σtotal
(Section 1.2.2). The corresponding likelihood then reads

 −1/2
det R
e 
1

T e −1
L β, σ 2 , θ, τ ; Y =



N/2 exp − 2 (Y − F β) R (Y − F β) , (1.46)
2
2πσtotal 2σtotal

where σ 2 and R in Eq. (1.41) have been replaced by σtotal

2 and R,
e respectively.
2
Similarly, the estimates for β and σtotal are written as
 −1
e −1 F
βb = β (θ) = F T R e −1 Y,
FTR (1.47)

1
σ 2
btotal 2
= σtotal (θ) = (Y − F β)T R
e −1 (Y − F β) . (1.48)
N
Based on Eqs. (1.46) to (1.48), the optimization problem can be formulated as follows
1h  
2

i
θ,
b τb = arg min log det R
e + N log 2πb
σtotal +N . (1.49)
θ∈Dθ ,τ ∈(0,1) 2

Notice that compared to Eq. (1.45), there is an additional parameter τ to be optimized and
that it is bounded in (0, 1).

6
The logarithm of the likelihood is usually taken in the optimization to avoid numerical underflow problem.

UQL AB-V2.1-105 - 13 -
UQL AB user manual

1.5.1.3 Kriging with known noise (homo- and heteroschedastic)

Finally, in the case of noisy responses with known noise variance, the covariance matrix is
given by Eq. (1.40).

An analytical form for the estimate of β as function of θ is nevertheless available (Rasmussen

and Williams, 2006) and reads

−1 T −1
βb = β θ, σ 2 , σn2 = F T C −1 F


F C Y. (1.50)

There is, however, no analytical expression for the Gaussian process variance σ 2 . and it
must be simultaneously optimized with the rest of the hyperparameters θ. The optimization
problem can be formulated as follows:
 
2 1  T
−1

θ, σ
b b = arg min log (det C) + N (log(2π) + N + Y − F β C
b Y − Fβ .
b
θ∈Dθ ,σ 2 ∈Dσ2 2
(1.51)
Note that, if the known noise is independent (i.e., a vector σn2 ), the noise covariance matrix
Σn is simplified to diag σn2 ; whereas, if the known noise variance is homoscedastic (i.e., a

scalar σn2 ), the noise covariance matrix Σn is simplified to σn2 I, where I is the identity matrix.

1.5.2 Cross-validation estimation

The general principle of a cross-validation (CV) method known as the K-fold cross-validation
is to split the whole data set of observations D = x(i) , y (i) , i = 1, . . . , N into K mutually


exclusive and collectively exhaustive subsets {Dk , k = 1, . . . , K} such that

K
[
Di ∩ Dj = ∅ , ∀(i, j) ∈ {1, . . . , K}2 and Dk = D. (1.52)
k=1

The k-th subset of cross-validated predictions is obtained by estimating the model using all
the subsets, except for the k-th one, and using the model to predict that specific k-th subset
that was left apart. The leave-one-out (LOO) CV procedure corresponds to the special case
that the number of subsets is equal to the number of observations (i.e., K = N ).

The cross-validation error of the k-th set is computed as

X   2
CV,k = y (i) − µYb x(i) ; β, σ 2 , θ, σn2 , D \ Dk
(x(i) ,y(i) )∈Dk
X   2 (1.53)
= y (i) − µYb ,\Dk x(i) ; β, σ 2 , θ, σn2
(x(i) ,y(i) )∈Dk

where µYb ,\Dk denotes the mean of the Kriging predictor (see, for instance, Eq. (1.6) for the
noise-free prediction) on the k-th CV subset x(i) , y (i) ∈ Dk conditioned on the other K − 1


UQL AB-V2.1-105 - 14 -
 Kriging (Gaussian process modeling)

subsets. In the case of LOO CV, the number of elements in Dk is equal to one.

The overall cross-validation error then reads

K
2 1 X
, θ, σn2 ; Y


CV β, σ = CV,k . (1.54)
N
k=1

The idea behind CV estimation method is to find the set of Kriging parameters β, σ 2 , θ, and,
if apply, σn2 , that minimizes the cross-validation error. As in the case of maximum-likelihood
estimation, depending on whether the noise in the response is present and the nature of the
noise, different optimization problems based on the CV estimation can also be set up.

1.5.2.1 Noise-free Kriging

For models with noise-free responses, the optimization problem in a CV estimation reduces
to (Santner et al., 2003; Bachoc, 2013)

θb = arg min CV (θ; Y) (1.55)

θ∈Dθ

Notice that in this case the CV error of Eq. (1.54) can be written only as a function of θ be-
cause of the analytical formula for β (Eq. (1.8)), the non-dependence of σ 2 from the predictor
mean (Eq. (1.6)), and the absence of σn2 .
b the estimate of β is computed as in Eq. (1.8), while the estimate of σ 2 is computed
Using θ,
using the following equation (Cressie, 1993; Bachoc, 2013)
  2
K y (i) − µYb ,\Dk x(i) ; θb
b = 1
X X
b2 = σ 2 (θ)
σ   (1.56)
N c2b x(i) ; θb
k=1 i∈Dk Y ,\Dk

where c2b denotes the normalized variance of the Kriging predictor on the k-th CV subset
Y ,\Dk
 (i) (i)
x ,y ∈ Dk conditional on all points of the other K − 1 subsets. In other words,

σ 2b
Y ,\Dk
c2Yb ,\D = (1.57)
k σ2
where σ 2b denotes the variance of the Kriging predictor (Eq. (1.7)) on the same k-th CV
Y ,\Dk
subset conditional on all points of the other K − 1 subsets.

1.5.2.2 Kriging with unknown homogeneous noise (homoscedastic)

In the case of noisy responses with unknown homogeneous noise variance, an additional pa-
rameter τ is to be optimized along with the hyperparameters θ as follows

θ,
b τb = arg min ECV (θ, τ ; Y) . (1.58)
θ∈Dθ ,τ ∈(0,1)

UQL AB-V2.1-105 - 15 -
UQL AB user manual

in which µYb ,\Dk x(i) ; θ, τ 2 given in Eq. (1.25) is used as the Kriging predictor.

2
Using θb and τb, the estimate of β is computed as in Eq. (1.47), while the estimate of σtotal is
computed using the following equation
  2
K X y (i) − µYb ,\Dk x(i) ; θ,b τb
b τb = 1
  X
2 2
σ
btotal = σtotal θ,   (1.59)
N c2b x(i) ; θ,
b τb
k=1 i∈Dk Y ,\Dk

where c2b above is analogous to the one appeared in Eq. (1.56) but based on Eq. (1.26),
Y ,\Dk
instead of Eq. (1.7).

1.5.2.3 Kriging with known noise (homo- and heteroschedastic)

Similarly to the case of the ML estimation (Section 1.5.1.3), models with known responses
noise variance also require that the Gaussian process variance σ 2 is optimized with rest of
the hyperparameters θ. The optimization problem now reads

b2 = arg min CV θ, σ 2 ; Y, Σn

θ,
bσ (1.60)
θ∈Dθ ,σ 2 ∈Dσ
2

where Σn is the known noise covariance matrix. Here, the Kriging predictor used in the
computation of the CV error uses the formulation of Eq. (1.16).

If the noise is independent, the noise covariance matrix Σn is simplified to diag σn2 , where

σn2 is the vector of noise variances; while if the noise is homoscedastic (i.e., a scalar σn2 ), the
noise covariance matrix Σn is simplified to σn2 I, where I is the identity matrix.

1.6 Optimization methods

To solve the optimization problems described in Section 1.5, there are trade-offs between
choosing some local (usually gradient-based) or choosing global (e.g., evolutionary algo-
rithms) methods. Local methods tend to converge faster and require fewer objective function
evaluations, but may perform poorly due to the existence of flat regions and multiple lo-
cal minimas, especially for increasing input dimension. Alternatively, the so-called hybrid
methods combine the features of both methods.

The currently available optimization methods are briefly discussed below:

• Interior point with L-BFGS Hessian approximation

An interior point gradient-based method is used to approximate the Hessian matrix us-
ing a limited-memory variant of the quasi-Newton Broyden-Fletcher-Goldfarb-Shanno
(L-BFGS) method (Byrd et al., 1999; Nocedal, 1980).
• Genetic Algorithm (GA)
GA is a well-established global optimization method (Goldberg, 1989) that has been
applied in many fields for the past decades. A hybrid approach can also be used, in

UQL AB-V2.1-105 - 16 -
 Kriging (Gaussian process modeling)

which the final solution of the genetic algorithm is used as a starting point of the
gradient-based method previously mentioned.
• Covariance matrix adaptation–evolution strategy (CMA-ES)
CMA-ES is a derandomized stochastic search algorithm introduced by Hansen and Os-
termeier (2001). It proceeds by adapting the covariance matrix of a normal distribution
such that directions that have improved the objective function in the recent past itera-
tions are more likely to be sampled again. Similar to GA, a hybrid approach can also
be used with CMAES.

An example of the objective function landscape for a one-dimensional problem is given in

−1
Figure 7. The original model is M(x) = 1 + 25x2 (i.e., the Runge function). The ex-
perimental design consists of 8 points as shown in Figure 7(a). The evaluations of objective
functions as a function of θ based on the ML and CV estimations (both are to be minimized)
are shown in Figures 7(c) and 7(b), respectively.

1.7 A posteriori error estimation

After the Kriging metamodel is set up, its predictive accuracy on new set of data can be
assessed either by using leave-one-out (LOO) cross-validation error or validation error.

1.7.1 Leave-one-out cross-validation error

The leave-one-out (LOO) cross-validation (CV) error is calculated on the initial experimental
design X , and its corresponding responses Y = M (X ) as follows

PN  2 
1  i=1 M(xi ) − µYb ,(−i) (xi ) 
LOO = (1.61)
N Var [Y]
 

where µYb ,(−i) (xi ) denotes the Kriging metamodel that is obtained using all the points of the
experimental design X , except xi .

1.7.2 Validation error

The validation
nerror is calculated
 as the relative
o generalization error on an independent set of
(i) (i)
data Dval = xval , yval , i = 1, . . . , Nval as follows
 2 
N val  
(i) (i)
MK
P
M(xval − xval
Nval − 1 
 i=1

val = (1.62)

 N 
Nval  P val 
(i)
 2  
M xval − µ̂Yval
i=1

UQL AB-V2.1-105 - 17 -
UQL AB user manual

(a) The Experimental Design, drawn from the Runge function.

(b) Maximum likelihood (ML) estimation objective function landscape.

(c) Cross-validation (CV) estimation objective function landscape.

Figure 7: Examples of one-dimensional objective function landscapes as functions of scale

parameter θ for ML and CV estimation methods using various Gaussian process correlation
families.

 
1 PNval (i)
where µ̂Yval = Nval i=1 M xval is the sample mean of the validation set responses. This
error measure is useful to compare the performance of different metamodels evaluated on
the same validation set.

UQL AB-V2.1-105 - 18 -
 Kriging (Gaussian process modeling)

1.8 Gaussian process trajectories

As mentioned in Section 1.2, a Kriging model is a Gaussian process defined by its mean
and covariance functions. From a Bayesian persperctive, they are the posterior mean and
covariance obtained by conditionning a prior Gaussian process, defined by its trend and auto-
correlation function, over existing observations. The posterior variance is given by Eq. (1.17)
(Note that this is simplified into Eq. (1.7) for interpolation and Eq. (1.26) for regression
with homoscedastic noise). Generalizing the expression given in Eq. (1.17), the posterior
covariance reads:

Σ x, x0 = C x, x0 − cT (x) C −1 c(x0 ) + uTc (x)(F T C −1 F )−1 uc x0 ,




(1.63)

where C (x, x0 ) = σ 2 R (x, x0 ; θ) is the prior covariance matrix.

For efficiency, this conditional Gaussian process is discretized using the Karhunen-Loève ex-
pansion (Ghanem and Spanos, 1991) and expansion optimal linear estimation (EOLE, Li and
Der Kiureghian (1993)). The two methods rely on the spectral decomposition of the covari-
ance matrix of the Gaussian process (See Section 1.2.1 for interpolation and Section 1.2.2
for regression). The reader is referred to the UQL AB User Manual – Random Fields for more
details on these two methods.

This feature is only available for models with scalar outputs. For multi-output problems, it is
recommended to build a single Kriging model for each output.

UQL AB-V2.1-105 - 19 -
Chapter 2

Usage

In this chapter, a reference problem is set up to showcase how each of the Kriging ingredients
described in Section 1.2.3 can be used in UQL AB.

2.1 Reference problem

In this manual, Kriging is used to build a surrogate of a model given a set of observed inputs
and model responses. To this end, the following one-dimensional function is used as the true
model
M (x) = x sin (x) , x ∼ U (0, 15) . (2.1)

2.2 Problem setup

The Kriging module creates a MODEL object. The basic options common to any Kriging meta-
model read
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';

Recalling (and slightly expanding) Eq. (1.1), a Kriging metamodel reads

P
X     
MK (x) = βj θb fj (x) + σ 2 Z x ; R θb , (2.2)
j=1

where θb is obtained by solving an optimization problem

θb = arg min J(θ) (2.3)

θ∈Dθ

In practice, the objective function J(θ) differs depending on the choice of estimation method
(i.e., maximum-likelihood or cross-validation). The main ingredients that need to be set up
to obtain a Kriging metamodel are the following:

20
 Kriging (Gaussian process modeling)

• An experimental design, X , and the corresponding model responses, Y. If they are not
available, they can be generated by defining a full computational (true) model and a
probabilistic input model1 .

• A trend specification, i.e., selection of the type of trend component functions fj (x) , j =
1, . . . , P .

• An appropriate correlation function R(x, x0 ; θ).

• A hyperparameter estimation method that defines the objective function J(θ) in Eq. (2.3);
this choice also defines the approach to estimate the Gaussian process variance σ 2 .

• An optimization method to solve the problem in Eq. (2.3).

Note that the minimal amount of information that a user needs to supply is the experimental
design and the corresponding model responses. The user can either select and tune each of
the ingredients, or leave them to their default values.

2.3 Kriging metamodel calculation: noise-free case

Below a minimal configuration example to create a Kriging metamodel in UQL AB is given as

was discussed in Section 2.2.
% S t a r t t h e UQLab framework
uqlab
rng(100,'twister') % For r e p r o d u c i b l e results

% Create experimental design

X = transpose(0:2:14);
Y = X.*sin(X);

% D e f i n e a K r i g i n g metamodel
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;

% C r e a t e t h e K r i g i n g metamodel
myKriging = uq_createModel(MetaOpts);

Once the metamodel is created, a report of the Kriging results can be printed on screen by:
uq_print(myKriging)

which then shows:

%--------------- Kriging metamodel ---------------%
1
Note that this step belongs to the general context of metamodeling, i.e., it is not specific to Kriging (see, for
instance, UQL AB User Manual – Polynomial Chaos Expansions)

UQL AB-V2.1-105 - 21 -
UQL AB user manual

Object Name: Model 1

Input Dimension: 1
Output Dimension: 1

Experimental Design
Sampling: User
X size: [8x1]
Y size: [8x1]

Trend
Type: ordinary
Degree: 0
Beta: [ 31.66776 ]

Gaussian Process (GP)

Corr. type: ellipsoidal
Corr. isotropy: anisotropic
Corr. family: matern-5_2
sigmaˆ2: 1.18220e+05
Estimation method: Cross-validation

Hyperparameters
theta: [ 2.90596 ]
Optim. method: Hybrid Genetic Alg.

GP Regression
Mode: interpolation

Error estimates
Leave-one-out: 5.55516e-01
%--------------------------------------------------%

It can be observed that the default values regarding the trend, correlation function, estima-
tion, and optimization method have been assigned. A visual representation of the metamodel
can be obtained by:
uq_display(myKriging)

The figure produced by uq_display is shown in Figure 8.

Note: Note that uq_display for Kriging MODEL objects can only be used for model one-
and two-dimensional inputs.

2.3.1 Accessing the results

The Kriging MODEL object myKriging created in Section 2.3 uses all the default configura-
tion options (as listed in the output of uq_print(myKriging)). This object can be directly
accessed to obtain various results and information of the metamodel.

UQL AB-V2.1-105 - 22 -
 Kriging (Gaussian process modeling)

Figure 8: The output of uq_display of a one-dimensional Kriging MODEL object.

2.3.1.1 Kriging parameters

The values of the basic parameters of the Kriging metamodel, namely θ, σ 2 , and β are con-
tained in the structure myKriging.Kriging:
struct with fields:

beta: 31.6678
sigmaSQ: 1.1822e+05
theta: 2.9060

2.3.1.2 Model evaluations

The experimental design and the corresponding model responses are stored in the
myKriging.ExpDesign structure:
struct with fields:

Sampling: 'User'
X: [8x1 double]
Y: [8x1 double]
NSamples: 8
U: [8x1 double]

UQL AB-V2.1-105 - 23 -
UQL AB user manual

The fields contain the following information:

• ExpDesign.Sampling: the source of the experimental design

• ExpDesign.NSamples: the size of the experimental design
• ExpDesign.X: the experimental design X
• ExpDesign.Y: the corresponding full model responses, Y = M(X )
• ExpDesign.U: the scaled experimental design (refer to Section 2.10 for details)

2.3.1.3 A posteriori error estimates

The Leave-One-Out (LOO) cross-validation error of the metamodel is stored in myKriging.Error:

struct with fields:

LOO: 0.5555

This error is calculated according tofollowing Eq. (1.61). Refer to Section 1.7.1 for details.

2.4 Kriging metamodel setup

In the following subsections, the various configuration options of each of the ingredients of
a Kriging metamodel are presented. These ingredients are summarized below:

• MetaOpts.ExpDesign contains the options regarding the generation or specification

of the experimental design. The way to use each option is discussed in Section 2.4.1
and list of all available options can be found in Table 5.
PP
• MetaOpts.Trend contains the options regarding the term j=1 βj (θ)fj (x) in Eq. (2.2),
see Section 1.3 for a theoretical introduction. The way to use each option is discussed
in Section 2.4.2 and the list of all available options can be found in Table 6.

• MetaOpts.Corr contains the options regarding the correlation function that is used in
order to compute R in Eq. (2.2), see Section 1.4 for a theoretical introduction. The
way to use each option is discussed in Section 2.4.3 and the list of all available options
can be found in Table 8.

• MetaOpts.EstimMethod refers to the method that is used for estimating the hyperpa-
rameters θ (maximum-likelihood or cross-validation). The choice of estimation method
corresponds to different ways of calculating σ 2 (θ) as well as J(θ). See Section 1.5 for
a theoretical introduction of each estimation method.

• MetaOpts.Optim contains the options related to the method for solving the optimiza-
tion problem in Eq. (2.3). See Section 1.6 for a brief description of the available opti-
mization methods. The way to use each option is discussed in Section 2.4.5 and the list
of all available options can be found in Table 10.

UQL AB-V2.1-105 - 24 -
 Kriging (Gaussian process modeling)

2.4.1 Specification and generation of experimental design

An experimental design for constructing a Kriging metamodel can either be specified using
existing data or generated from a MODEL and an INPUT objects.

2.4.1.1 Using existing data

Two alternatives are offered depending on where the experimental design data is stored:

• If data is stored in the M ATLAB workspace, e.g., in the variables X, Y:

MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;

• If data is stored in a .mat file, e.g., mydata.mat:

MetaOpts.ExpDesign.Sampling = 'Data';
MetaOpts.ExpDesign.DataFile = 'mydata.mat';

Note: The experimental design that is contained in a .mat file must be saved with
variable names X and Y, respectively.

The input and output dimensions of the problem M, Nout are automatically inferred from the
number of columns of X and Y, respectively.

2.4.1.2 Using INPUT and MODEL objects

In order to use INPUT object, it first needs to be created. For the reference problem in Eq. (2.1)
this reads:
InputOpts.Marginals.Type = 'Uniform';
InputOpts.Marginals.Parameters = [0 15];
myInput = uq_createInput(InputOpts);

The input and output dimensions of the problem are automatically inferred form the con-
figuration of the INPUT object myInput and MODEL object myModel. For details about the
configuration options available for INPUT and MODEL objects, please refer to the UQL AB User
Manual – the INPUT module and the UQL AB User Manual – the MODEL module, respectively.

Then, a MODEL object is created as follows:

ModelOpts.mString = 'X.*sin(X)';
myModel = uq_createModel(ModelOpts);

Now, the INPUT and MODEL objects are specified in to the Kriging metamodel configuration:
MetaOpts.Input = myInput;
MetaOpts.FullModel = myModel;

UQL AB-V2.1-105 - 25 -
UQL AB user manual

Finally, the size of the experimental design is specified, e.g., for N = 8 sample points:
MetaOpts.ExpDesign.NSamples = 8;

Note that, by default, Latin Hypercube sampling (LHS) is used in order to obtain X . However,
other sampling strategies can be selected through the option MetaOpts.ExpDesign.Sampling,
see Table 5 in Section 3.1.1 for a list of the available sampling strategies.

2.4.2 Trend

2.4.2.1 Standard options

Some typical configuration options are the following:

• When considering a constant (ordinary Kriging), linear, or quadratic trend, the field
Metaopts.Trend.Type must contain the appropriate string value, that is, 'ordinary',
'linear', and 'quadratic' respectively.

• For a polynomial trend of an arbitrary degree (Eq. (1.27)), set the following option:
MetaOpts.Trend.Type = 'polynomial';
MetaOpts.Trend.Degree = q;

where q is an integer equal to the polynomial degree.

• For a constant trend with fixed value (i.e., simple Kriging), set the following option:
MetaOpts.Trend.Type = 'simple';
MetaOpts.Trend.CustomF = V;

where V is a real number. If the option .Trend.customF is not specified, the default
value 0 is set.

A summary of the available trend options in UQL AB along with the corresponding formula
and the additionally required options is given in Table 2.

Table 2: Trend options combinations

Trend.Type Formula Trend.Degree Trend.CustomF
'simple' f (x) (no trend estimation) Not required Not required
'ordinary' β0 Not required Not required
PM
'linear' β0 + i=1 βi xi Not required Not required
PM PM PM
'quadratic' β0 + i=1 βi xi + i=1 j=1 βij xi xj Not required Not required
'polynomial' see Eq. (1.27) Required Not required
PP
'custom' i=1 βi fi (x) Not required Required

In Figure 9, the mean and variance of various Kriging predictors having different trend con-

UQL AB-V2.1-105 - 26 -
 Kriging (Gaussian process modeling)

Figure 9: The output of the example script uq_Example_Kriging_03_TrendTypes. The

mean and variance of various Kriging predictors are plotted, using different trend configura-
tions. In the top figure, black line indicates the true model prediction.

figurations are plotted as in the example script uq_Example_Kriging_03_TrendTypes.

Note: The field MetaOpts.Trend.Type determines the trend type of a Krig-

ing metamodel. If the user does not define a trend type, the default
Metaopts.Trend.Type = 'ordinary' (i.e., ordinary Kriging) is used.

2.4.2.2 Advanced options

A user can define arbitrary custom trends via advanced options. First, the trend type is
specified as follows:
Metaopts.Trend.Type = 'custom';

Then the trend can either be defined as a string or a function handle. If, for example, the
trend is a sine function sin(x) then the options is set as follows:
MetaOpts.Trend.CustomF = 'sin';

or:
MetaOpts.Trend.CustomF = @sin;

The function can be any of the built-in M ATLAB (as in the example above), UQL AB, or user-

UQL AB-V2.1-105 - 27 -
UQL AB user manual

defined functions (as long as the function is available in the M ATLAB path).

In general, the dimension of the information matrix F is N × P . Depending on the trend

type (Metaopts.Trend.Type) the value of P varies:

• If MetaOpts.Trend.Type is 'simple', then P equals to 1.

• If MetaOpts.Trend.Type is either 'ordinary', 'linear', or 'quadratic', then P

(M +2)(M +1)
equals to 1, M + 1, 2 , respectively.

• If MetaOpts.Trend.Type='polynomial', then depending on the polynomial degree

q, defined in Metaopts.Trend.Degree, P equals to Mq+q .

• If MetaOpts.Trend.Type is 'custom', then P depends on how Metaopts.Trend.CustomF

is specified:

– If MetaOpts.Trend.CustomF is a string or double, then P equals to 1.

– If MetaOpts.Trend.CustomF is a cell array of strings or function handles, then
P equals the length of the cell array. In other words, each element of the array
corresponds to a single column of F .

2.4.3 Correlation function

The key ingredients for specifying a correlation function is done through the standard op-
tions, while urther flexibility such as using a custom correlation function can be accessed
through the advanced options.

2.4.3.1 Standard options

The three key ingredients for specifying a correlation function are:

• Type that specifies the type of the autocorrelation function as described in Section 1.4.2.
By default, the type is set to 'ellipsoidal'. To use instead a separable correlation
function, the following option is set:
MetaOpts.Corr.Type = 'separable';

Note: If the user does not specify the correlation type, it is, by default, set to be
'ellipsoidal'.

• Family that specifies the one-dimensional correlation function family R(·) of the Gaus-
sian process as described in Section 1.4. Its inputs depend on the type of correlation
function:

– For ellipsoidal correlation functions, the correlation family is a function of the form
R(h) where h > 0 corresponds to the normalized Euclidean distance between x

UQL AB-V2.1-105 - 28 -
 Kriging (Gaussian process modeling)

and x0 (see Eq. (1.34)). Each coordinate difference between the rows of x and x0
is normalized by a positive scalar θ.
– For separable correlation functions, the correlation family is a function of the form
R(x, x0 ; θ) where θ is a positive scalar.

To use a built-in correlation family, its name must be set, see Table 8 in Section 3.1 for
the name of each built-in correlation family and Section 1.4.1 for a brief description of
each family. For example, to use the exponential family, the following option is set:
MetaOpts.Corr.Family = 'exponential';

Note: By default, the correlation family is set to be 'matern-5_2'.

• Isotropic a flag (with true or false value) that specifies whether the correlation func-
tion is isotropic or not. By default, the correlation function is considered anisotropic.
To change the default, the following option is set:
MetaOpts.Corr.Isotropic = true;

Note: By default, the correlation function is considered anisotropic, that is,

MetaOpts.Corr.Isotropic = false.

In Figure 10, the mean and variance of various Kriging predictors having different correlation
families (as in the example script uq_Example_Kriging_01_1D) are plotted.

2.4.3.2 Advanced options

Through the advanced options, users can specify custom correlation function families, custom
routines to compute the whole correlation matrix R, as well as a nugget term to stabilize the
inversion of R.

• User-defined correlation family: Apart from the built-in one-dimensional correlation

families, users can also defined one themselves. Depending on the selected type of
correlation function (i.e., separable or ellipsoidal), UQL AB expects different inputs for
the user-defined correlation family:

– Separable case: UQL AB expects that a user-defined correlation family to be of the

form R (x, x0 ; θ) where x, x0 and θ are scalars and a positive scalar, respectively.
For example, to specify the following correlation family:
"   #
x − x 0 1.5
R x, x0 ; θ = exp −


,
θ

the following options are set:

UQL AB-V2.1-105 - 29 -
UQL AB user manual

Figure 10: The output of the example script uq_Example_Kriging_01_1D. The mean and
variance of various Kriging predictors are plotted, using different correlation families. In the
top figure, black line indicates the true model prediction.

MetaOpts.Corr.Type = 'separable';
MetaOpts.Corr.Family = @(x1,x2,th) exp(-((x1-x2)/th).ˆ1.5);

Notice that for the separable case, a user-defined correlation function can also be
non-stationary.

Note: UQL AB expects that the user-defined correlation family to be vector-

ized, i.e., it should be able to handle inputs x1 and x2 that are vectors,
either with the same or different length.

– Ellipsoidal case: UQL AB expects that an ellipsoidal user-defined correlation fam-

ily to be of the form R(h) where h is a positive scalar. UQL AB calculates the value
of h following Eq. (1.34) and then sends it to the user-defined function.
For example, to specify the following correlation family:

R (h) = exp −h1.5 ,

the following options are set:

MetaOpts.Corr.Type = 'ellipsoidal';
MetaOpts.Corr.Family = @(h) exp(-h.ˆ1.5);

UQL AB-V2.1-105 - 30 -
 Kriging (Gaussian process modeling)

Note: UQL AB expects that the user-defined correlation family to be vector-

ized, i.e., it should be able to handle input h that is vector.

• User-specified routine to calculate the correlation matrix (R): All the aforemen-
tioned options regarding the correlation function are handled by the built-in function
uq_eval_Kernel. Users may instead specify a fully custom correlation function that
computes the correlation matrix R directly as long as the following rules are followed:

– The function should return the correlation matrix R, accept x1 , x2 , θ, and a struc-
ture of options (in that order) such as in the following function declaration:
function R = my_eval_R(x1, x2, theta, Options)

Note: The order of the input parameters in the custom function matters and
must be specified as is.

– Inputs x1 and x2 are matrices with arbitrary number of rows and M number of
columns (where M is the input dimension).

– When the custom correlation function my_eval_R is called, for instance, during
the calculation of or prediction with a Kriging metamodel, the structure Options
will contain all options that were set inside MetaOpts.Corr.

– The input theta of the function my_eval_R corresponds to the hyperparameters

θ and it is expected to be a vector of arbitrary length. There is no restriction
regarding its length but its dimension should be reflected in the choice of either
initial value and/or optimization bounds (see Section 2.4.5 for more information
about the optimization options).

This custom routine can then be specified in the Kriging metamodel options as:
MetaOpts.Corr.Handle = @my_eval_R;

• Nugget: A nugget is a small numerical value added to the diagonal of the R matrix to
improve the stability of its numerical inversion. To add a constant nugget, say ν = 10−4 ,
the following option is set:
MetaOpts.Corr.Nugget = 1e-4;

If the nugget is specified as a vector ν of length N , each of its elements is added to the
corresponding diagonal element of R, that is, Rii = 1 + νi , {i = 1, . . . , N }.

Note: By default, the nugget value in UQL AB is set to ν = 10−10 .

UQL AB-V2.1-105 - 31 -
UQL AB user manual

Figure 11: The output of the example script uq_Example_Kriging_02_VariousMethods.

The mean and variance of various Kriging predictors are plotted, using two different op-
timization methods. On the left-hand side, the maximum likelihood estimation method is
used and on the right-hand side the LOO cross-validation method is used.

2.4.4 Estimation methods

The estimation method determines the objective function that is minimized in Eq. (2.3) to
estimate the hyperparameters. If required, this choice also determines how the (constant)
Gaussian process variance σ 2 is calculated. For details, see Section 1.5.

The maximum likelihood and cross-validation methods are available for estimating the Krig-
ing parameters and it suffices to select either 'ML' or 'CV' in the field Metaopts.EstimMethod,
for the maximum-likelihood and cross-validation estimation methods, respectively.

An example that illustrates different Kriging predictors created using different estimation is
given in the script uq_Example_Kriging_02_VariousMethods (Figure 11).

2.4.4.1 Advanced Options

If the cross-validation (CV) method is selected as the estimation method, then by default,
UQL AB uses the N -fold CV (or equivalently, the leave-one-out) approach. Other K-fold CV
approaches (K ≤ N ) can be used by changing the number of sample points left out from the
estimation via the MetaOpts.CV.LeaveKOut option as follows:
MetaOpts.CV.LeaveKOut = k;

where k is an integer. The number of folds (subsets) in K-fold CV is computed as follows

K = dN/ke (2.4)

UQL AB-V2.1-105 - 32 -
 Kriging (Gaussian process modeling)

As an example, for N = 100 and k = 3, the number of folds in K-fold CV equals to 34.

Note: By default, the estimation method is set to be the N -fold (leave-one-

out) cross-validation method, that is, MetaOpts.EstimMethod = 'CV' and
Metaopts.CV.LeaveKOut = 1.

Note: If K-fold CV is used, UQL AB randomly permutes the experimental design before
creating the folds. Furthermore, the MetaOpts.CV.LeaveKOut must be smaller
than the size of the experimental design N , such that at least two folds are avail-
able for cross-validation. Otherwise UQL AB will throw an error.

2.4.5 Optimization methods

Various configuration options associated with the method to solve the optimization prob-
lem in Eq. (2.3) are available. The available optimization methods can be divided into two
categories:

• Gradient-based methods
Currently, only BFGS method is available as a local gradient-based method. To use this
method, the following option is set:
MetaOpts.Optim.Method = 'BFGS';

As in any gradient-based methods, an initial value of the hyperparameters (denoted

here as θ0 ) is required. By default, θ0 = 1 in each dimension. However, a different
initial value can be selected, e.g., θ0 = 0.5 (in each dimension):
MetaOpts.Optim.InitialValue = 0.5;

Furthermore, different initial value per dimension can also be specified using an M × 1
vector instead of a scalar.

Note: Note that, by default, the hyperparameters are defined in to the scaled
(auxiliary) space U as described in Section 2.10.

Note: The BFGS method requires a separate license for the Optimization Toolbox
of M ATLAB.

• Global and hybrid methods

The Genetic Algorithm (GA), the covariance matrix adaptation–evolution strategy (CMA-
ES), and their hybrid counterparts (HGA and HCMAES, respectively) are available as
optimization methods in UQL AB. To use one of these methods, e.g., GA, the following
option is set:
MetaOpts.Optim.Method = 'GA';

UQL AB-V2.1-105 - 33 -
UQL AB user manual

As in any global methods, the bounds of the optimization variable(s) needs to be set.
By default, the bounds (lower, upper) [0.001, 10]T are used. To specify different bounds,
for example, [0.01, 1]T , the following options is set:
MetaOpts.Optim.Bounds = [0.01; 1];

Note: Note that, by default, the hyperparameters are defined in to the scaled
(auxiliary) space U as described in Section 2.10.

Additional options which are specific to each method may exist. For example, when
using the GA method, one might need to set a different value of stall generations, e.g.,
20. This, in turn, can by accomplished by setting:
MetaOpts.Optim.GA.nStall = 20;

Note: GA requires a separate license for the M ATLAB Global Optimization Toolbox
and all hybrid methods require a separate license for the M ATLAB Optimiza-
tion Toolbox. CMAES, in contrast, is built into UQL AB and thus requires no
additional license.

Moreover, regardless of the method, users can manually specify the following options:

• The number of iterations or generations (its actual meaning depends on the optimiza-
tion method), for example:
MetaOpts.Optim.MaxIter = 100;

• The convergence tolerance, for example:

MetaOpts.Optim.Tol = 1e-5;

• The verbosity of the optimization process, e.g., to show only the final result:
MetaOpts.Optim.Display = 'final';

For a list of all the available options for each method, refer to Table 10 in Section 3.1.5.

Note: If the user does not specify any method and if the required toolboxes are avail-
able, the Hybrid Genetic Algorithm (HGA) is used by default. If none of the
toolboxes are available, the default is the built-in covariance matrix adaptation–
evolution strategy (CMA-ES).

If no optimization of the hyperparameters is required, due to, for example, a prescribed value
theta_user is used, the following option is set:
MetaOpts.Optim.Method = 'none';
MetaOpts.Optim.InitialValue = theta_user;

UQL AB-V2.1-105 - 34 -
 Kriging (Gaussian process modeling)

A comparison of the resulting Kriging predictors using different optimization (and estima-
tion) methods can be found in the example script uq_Example_Kriging_02_VariousMethods
(see Figure 11).

2.5 Kriging metamodels with noise (regression)

To demonstrate the Kriging metamodel with noisy responses, we consider once more the
reference problem described in Section 2.1 and add to the observed response an additive
noise
y = M (x) + ε (2.5)

The following subsections demonstrate how to create several Kriging metamodels in UQL AB
for different cases of noisy responses (see Section 1.2.2).

2.5.1 Unknown homoscedastic noise

Assuming that the noise is homoscedastic (i.e. its variance is constant), UQL AB can estimate
the unknown noise variance from the observed data during the calculation of the Kriging
model. To estimate the noise, the following option is set:
MetaOpts.Regression.SigmaNSQ = 'auto';

The example below demonstrates how to create a Kriging model with an unknown ho-
moscedastic noise in the response. First, an illustrative noisy data set is created:
% Create experimental design with noisy responses
rng(100,'twister') % f o r r e p r o d u c i b l e r e s u l t s
X = linspace(0,14,100)';
noiseVar = 3.0;
Y = X.*sin(X) + sqrt(noiseVar) * randn(size(X,1),1);

Then, a Kriging metamodel as in Section 2.3 is specified, now with the additional option:
% D e f i n e a K r i g i n g metamodel
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
% Estimate the noise
MetaOpts.Regression.SigmaNSQ = 'auto';

Finally, the metamodel is created:

myKriging = uq_createModel(MetaOpts);

The report produced by uq_print now displays an additional information about the noise
estimation process. Notice that the estimate noise variance is now reported.

UQL AB-V2.1-105 - 35 -
UQL AB user manual

%---------------- Kriging metamodel ----------------%

Object Name: Model 1
Input Dimension: 1
Output Dimension: 1

Experimental Design
Sampling: User
X size: [100x1]
Y size: [100x1]

Trend
Type: ordinary
Degree: 0
Beta: [ 4.18461 ]

Gaussian Process (GP)

Corr. type: ellipsoidal
Corr. isotropy: anisotropic
Corr. family: matern-5_2
sigmaˆ2: 6.43875e+01
Estimation method: Cross-validation

Hyperparameters
theta: [ 0.82604 ]
Optim. method: Hybrid Genetic Alg.

GP Regression
Mode: regression
Est. noise: true
sigmaNˆ2: 2.89456e+00

Error estimates
Leave-one-out: 9.03089e-02
%---------------------------------------------------%

A visual representation of the Kriging metamodel can be obtained via the uq_display func-
tion whose result is shown in Figure 12.

2.5.2 Known noise

If the noise variance is known a priori, it can be used directly in the regression model without
being estimated. To specify a known noise variance, the following options are set:
MetaOpts.Regression.SigmaNSQ = constSigmaNSQ;

where constSigmaNSQ is the known noise variance. The specification constSigmaNSQ de-
pends on the nature of the noise as follows:

• Homoscedastic noise: constSigmaNSQ is a scalar. This is the case in which the noise in
the response is constant everywhere.
• Independent heteroscedatic noise: constSigmaNSQ is a N × 1 column vector. This is the
case in which the noises are independent, but may differ at each observation point.

UQL AB-V2.1-105 - 36 -
 Kriging (Gaussian process modeling)

Figure 12: The output of uq_display of a Kriging MODEL object having a one-dimensional
input, with noisy responses.

• Correlated heteroscedastic noise: constSigmaNSQ is a N × N matrix. This is the case in

which the noises are neither constant nor independent. In this case, constSigmaNSQ is
the covariance matrix of the noise.

The case of known homoscedastic noise is illustrated using the previous example with smaller
set of data as follows:
% Create a h y p o t h e t i c a l experimental design with noisy responses
rng(100,'twister')
X = linspace(0,14,50)';
noiseVar = 3;
Y = X.*sin(X) + sqrt(noiseVar) * randn(size(X,1),1);
% D e f i n e a K r i g i n g metamodel
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
% Impose known h o m o s c e d a s t i c n o i s e v a r i a n c e
MetaOpts.Regression.SigmaNSQ = noiseVar;
% C r e a t e t h e K r i g i n g metamodel
myKriging = uq_createModel(MetaOpts);

UQL AB-V2.1-105 - 37 -
UQL AB user manual

Finally, below is another example of GP regression for an independent heteroscedastic noise

case in which the noise variances at individual observation locations are given as vector:
% Create a h y p o t h e t i c a l experimental design with noisy responses
X = linspace(0,14,15)';
% Define a vector of noise variance at i n d i v i d u a l data l o c a t i o n s
noiseVar = [0.3 0.4 4 0.25 0.16 ...
0.133 0.5 0.9 5 0.1600 ...
0.571 0.02 0.0225 0.02 0.8]';
% Create noisy responses data s e t with the noise
Y = X.*sin(X) + sqrt(noiseVar) .* randn(size(X,1),1);
% D e f i n e a K r i g i n g metamodel
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
% Impose t h e known h e t e r o s c e d a s t i c n o i s e v a r i a n c e
MetaOpts.Regression.SigmaNSQ = noiseVar;
% C r e a t e t h e K r i g i n g metamodel
myKrigingRegression = uq_createModel(MetaOpts);

The results of both cases are visualized in Figure 13. Notice that in Figure 13(b), the noise
variance associated for each experimental design (observation) points differs. A Kriging ex-
ample for the different cases of noise variance specification is provided in uq_Example_Kriging_08_Regressio

(a) known homoscedastic noise (b) known heteroscedastic noise

Figure 13: Example of one-dimensional GP regression models predictions with known ho-
moscedastic and independent heteroscedastic noise variances.

2.6 Kriging metamodels of vector-valued models

All the examples presented so far in this chapter dealt with scalar-valued models. If the model
(or when manually specified, the experimental design) is having vector-valued (multiple)

UQL AB-V2.1-105 - 38 -
 Kriging (Gaussian process modeling)

outputs, UQL AB performs an independent Kriging metamodel for each output component
with shared (common) experimental design and trend basis functions. No additional config-
uration is needed to enable this behavior. A Kriging example for model with multiple outputs
can be found in the UQL AB example script uq_Example_Kriging_07_MultipleOutputs.

2.6.1 Accessing the results

Calculating a Kriging metamodel on a model with multiple outputs will result in an output
structure with multiple elements (i.e., a structure array). As an example, a model with nine
outputs will produce the following output structure myKriging.Kriging:
1x9 struct array with fields:

beta
sigmaSQ
theta

Each element of the Kriging structure is functionally identical to its scalar counterpart in
Section 2.3.1. Similarly, the myKriging.Error structure becomes a structure array:
1x9 struct array with fields:

LOO

2.7 Using a validation set

If an independent validation set is provided (see Table 19 in Section 3.1.8), UQL AB auto-
matically computes the validation error according to Eq. (1.62). To provide a validation set
stored in, for example, XVal and YVal, the following option is set:
MetaOpts.ValidationSet.X = XVal;
MetaOpts.ValidationSet.Y = YVal;

The value of the validation error is stored in myKriging.Error.Val (see Table 22) and will
also be displayed when invoking uq_print(myKriging).

2.8 Using Kriging predictor as a model

Regardless of the configuration options that are used to construct a Kriging metamodel (as
described in Section 2.4), the resulting Kriging metamodel can be used to predict new points
according to, for instance in the case of noise-free responses, Eqs. (1.6)-(1.7). Indeed, after a
Kriging MODEL object is created in UQL AB, it can be used just like an ordinary UQL AB MODEL
object (for details, see the UQL AB User Manual – the MODEL module).

Consider the example in Section 2.1. After obtaining a Kriging metamodel, users can now
evaluate the mean of the Kriging predictor on point x = 1 as follows:

UQL AB-V2.1-105 - 39 -
UQL AB user manual

>> x = 1;
>> YmuKG = uq_evalModel(x)

YmuKG =

1.7544

The variance of the predictor can be also evaluated by using a two-component output:
>> [YmuKG,YvarKG] = uq_evalModel(x)

YmuKG =

1.7544

YvarKG =

3.0040

As most functions within UQL AB, model evaluations are vectorized. Therefore, evaluating
multiple points at a time is much faster than repeatedly evaluating one point at a time. For
example:
X = transpose(0:0.1:15)';
% E v a l u a t e t h e mean and v a r i a n c e o f t h e K r i g i n g p r e d i c t o r on p o i n t s x
[YmuKG,YvarKG] = uq_evalModel(X);

2.9 Specifying manually a Kriging predictor (predictor-only mode)

The Kriging module in UQL AB can also be used to build custom Kriging-based models which
may be used as predictors (as in Section 2.8). This allows, for instance, to import a Kriging
metamodel calculated with another software into the UQL AB framework, or to create an
ad-hoc Kriging model.

Below is an example of how to create a custom Kriging metamodel that is similar to the one
that was obtained in Section 2.3.
% Create the experimental design
X = transpose(0:2:15);
% C a l c u l a t e t h e r e s p o n s e o f t h e model
Y = X.*sin(X);

% D e f i n e a custom K r i g i n g o b j e c t
MetaOpts.Type = 'Metamodel';
MetaOpts.MetaType = 'Kriging';
MetaOpts.ExpDesign.Sampling = 'User';
MetaOpts.ExpDesign.X = X;
MetaOpts.ExpDesign.Y = Y;
% Select ordinary Kriging
MetaOpts.Kriging.Trend.Type = 'ordinary';

UQL AB-V2.1-105 - 40 -
 Kriging (Gaussian process modeling)

% S e t t h e v a l u e s o f b e t a , sigma ˆ2 and t h e t a
MetaOpts.Kriging.beta = 69.84;
MetaOpts.Kriging.sigmaSQ = 2.566e5;
MetaOpts.Kriging.theta = 9.999;
% Select e l l i p s o i d a l , Matern −3/2 c o r r e l a t i o n function
MetaOpts.Kriging.Corr.Type = 'ellipsoidal';
MetaOpts.Kriging.Corr.Family = 'matern-3_2';

% C r e a t e t h e metamodel
myCustomKriging = uq_createModel(Metaopts);

% E v a l u a t e t h e metamodel on some new p o i n t s

Xnew = transpose(1:2:13);
Ynew = uq_evalModel(Xnew);

If the metamodel has more than one output, it is sufficient to specify the same information
for each output by adding an index i to the MetaOpts.Kriging(i) structure array. Further-
more, to define a custom regression model, a value to an additional field MetaOpts.Kriging.sigmaNSQ
must be specified.

Note: By default, the nugget value is set to ν = 10−10 when creating a custom Kriging
metamodel (see Section 2.4.3.2 on nugget).

2.10 Performing Kriging on an auxiliary space (scaling)

The Kriging module offers various options for scaling the experimental design before calcu-
lating the metamodel. That is, instead of X , UQL AB may use the scaled experimental design
U. Depending on the value of the option Metaopts.Scaling and whether a probabilistic
input model has been defined (in Metaopts.Input), the transformation X 7→ U varies. All
the possible cases are summarized in Table 3.

Table 3: Available experimental design transformations

Input
No INPUT defined INPUT defined
Scaling
true U = (X − µX ) /σX U = (X − µX ) /σX
false U =X U =X
INPUT object - U = τ (X )

By setting MetaOpts.Scaling = true, the experimental design is scaled so that it has a zero
mean and a unit variance (i.e., standardized). If no input model has been specified, then the
mean and variance (µX , σX ) are empirically estimated from the available data specified in
Metaopts.ExpDesign.X. If, however, an input model has been specified via an INPUT object,
say myInput, as follows:
MetaOpts.Input = myInput;

UQL AB-V2.1-105 - 41 -
UQL AB user manual

then (µX , σX ) are estimated from the moments of the marginal distribution of each compo-
nent in X.

If an additional INPUT object, say myScaledInput, is set to the MetaOpts.Scaling option:

MetaOpts.Scaling = myScaledInput;

then U = τ (X ) where τ denotes the generalized isoprobabilistic transformation between the

two probability spaces (for more information, refer to the UQL AB User Manual – the INPUT
module).

Note: By default, Metaopts.Scaling = true, that is the experimental design is

scaled.

2.11 Drawing sample paths from a Gaussian process posterior

To draw sample paths of a Kriging model using spectral decomposition of the Kriging covari-
ance matrix as described in Section 1.8, the user needs to enable this feature by using the
following command:
MetaOpts.isTrajectory = true;

All options of the random field module can be accessed through the field .GRF. For instance,
the user may define the discretization domain by using the field .GRF.Domain as follows:
MetaOpts.GRF.Domain = [0; 15];

Note: In general, default values from the random field module are considered, except
for the following:

– Metaopts.GRF.DiscScheme = 'EOLE'
– MetaOpts.GRF.EOLE.Sampling = 'LHS'
– MetaOpts.GRF.EOLE.SPC = 20 (for M = 1)
– MetaOpts.GRF.EOLE.NSamples = 2000 (for M > 1)
– MetaOpts.GRF.EnergyRatio = 0.995
– MetaOpts.GRF.Domain = [minX_ED-0.1*R maxX_ED+0.1*R]

where minX_ED = min (X ), maxX_ED = max (X ) and R is the range of the exper-
imental design, i.e., R = max (X ) − min (X )
The same options are used when the user defines 'KL' as discretization scheme.

To generate trajectories, the uq_evalModel command needs to be called with an additional

argument that represents the desired number of trajectories or replications:
R = 100;
Ytraj = uq_evalModel(X,R);

UQL AB-V2.1-105 - 42 -
 Kriging (Gaussian process modeling)

For an input matrix X of size N × M , the output matrix Ytraj is now a three-dimensional
array of size N × Nout × R, where the third dimension corresponds to the R trajectories.

Given the myKriging object that was calculated in Section 2.3, a set of 15 sample paths is
plotted in Figure 14 together with the mean of the Kriging predictor.

Figure 14: The mean of the Kriging predictor together with 15 sample paths of the posterior
Gaussian process.

2.12 Using Kriging with constant inputs

In some analyses, users may need to assign a constant value to one or more input vari-
ables. When this is the case, the Kriging metamodel is calculated by internally removing the
constant parameters from the full inputs. This process is transparent to users as the model
is still evaluated using the full set of parameters (including those which are set constant).
UQL AB will automatically and appropriately account for the set of input parameters which
are declared constant.

To set a parameter to constant, e.g., with constValue as its value, the following options are
specified (for details, see UQL AB User Manual – the INPUT module):
InputOpts.Marginals.Type = 'Constant';
InputOpts.Marginals.Parameters = constValue;

UQL AB-V2.1-105 - 43 -
UQL AB user manual

Furthermore, when the standard deviation of an input parameter is set to zero, UQL AB auto-
matically sets the marginal of this parameter to the type Constant. For example, the follow-
ing uniformly distributed variable whose upper and lower bounds are identical (therefore its
standard deviation is zero) is automatically set to a constant with value 1:
InputOpts.Marginals.Type = 'Uniform';
InputOpts.Marginals.Parameters = [1 1];

UQL AB-V2.1-105 - 44 -
Chapter 3

Reference List

How to read the reference list

Structures play an important role throughout the UQL AB syntax. They offer a natural way
to semantically group configuration options and output quantities. Due to the complexity of
the algorithms implemented, it is not uncommon to employ nested structures to fine-tune
the inputs and outputs. Throughout this reference guide, a table-based description of the
configuration structures is adopted.

The simplest case is given when a field of the structure is a simple value or array of values:

Table X: Input
.Name String A description of the field is put here

which corresponds to the following syntax:

Input.Name = 'My Input';

The columns, from left to right, correspond to the name, the data type and a brief description
of each field. At the beginning of each row a symbol is given to inform as to whether the
corresponding field is mandatory, optional, mutually exclusive, etc. The comprehensive list
of symbols is given in the following table:

Mandatory
 Optional
⊕ Mandatory, mutually exclusive (only one of
the fields can be set)

Optional, mutually exclusive (one of them
can be set, if at least one of the group is set,
otherwise none is necessary)

When one of the fields of a structure is a nested structure, a link to a table that describes the

45
UQL AB user manual

available options is provided, as in the case of the Options field in the following example:

Table X: Input
.Name String Description
 .Options Table Y Description of the Options
structure

Table Y: Input.Options
.Field1 String Description of Field1
 .Field2 Double Description of Field2

In some cases, an option value gives the possibility to define further options related to that
value. The general syntax would be:
Input.Option1 = 'VALUE1' ;
Input.VALUE1.Val1Opt1 = ...;
Input.VALUE1.Val1Opt2 = ...;

This is illustrated as follows:

Table X: Input
.Option1 String Short description
'VALUE1' Description of 'VALUE1'
'VALUE2' Description of 'VALUE2'

.VALUE1 Table Y Options for 'VALUE1'

.VALUE2 Table Z Options for 'VALUE2'

Table Y: Input.VALUE1
 .Val1Opt1 String Description
 .Val1Opt2 Double Description

Table Z: Input.VALUE2
 .Val2Opt1 String Description
 .Val2Opt2 Double Description

Note: In the sequel, double and doubles mean a real number represented in double
precision and a set of such real numbers, respectively.

UQL AB-V2.1-105 - 46 -
 Kriging (Gaussian process modeling)

3.1 Create a Kriging metamodel

Syntax

myKriging = uq_createModel(MetaOpts)

Input

The struct variable MetaOpts contains the configuration information for a Kriging meta-
model. The detailed list of available options is reported in Table 4.

Table 4: MetaOpts
.Type 'Metamodel' Select the metamodeling tool
.MetaType 'Kriging' Select Kriging
 .Name String Unique identifier for the metamodel
 .Display String Level of information displayed
default: 'standard' during calculation
'quiet' Minimum display level, displays
nothing or very few information.
'standard' Default display level, shows the most
important information.
'verbose' Maximum display level, shows all the
information on runtime, like updates
on iterations, etc.
 .ExpDesign Table 5 Option to specify an experimental
design (Section 2.4.1)
 .Input String or INPUT object Option to specify an INPUT object
that describes the inputs of the
metamodel
String Name of the INPUT object
INPUT object Variable that refers to the INPUT
object
 .FullModel String or MODEL object Option to specify the MODEL object
that describes the full computational
model of the metamodel. The object
is used to compute the model
responses.
String Name of the MODEL object
MODEL object Variable that refers to the MODEL
object
 .Trend Table 6 Options to specify the Kriging trend
(Sections 1.3 and 2.4.2)

UQL AB-V2.1-105 - 47 -
UQL AB user manual

 .Corr Table 8 Options to specify the correlation

function (Section 2.4.3)
 .EstimMethod String Select the method to estimate
default: 'CV' Kriging parameters (Section 2.4.4)
'CV' Cross-validation estimation
(Section 1.5.2)
'ML' Maximum-likelihood estimation
(Section 1.5.1)
 .CV Table 9 Options relevant to the
cross-validation estimation method.
It is only taken into account if the
selected method is cross-validation.
 .Optim Table 10 Options related to the optimization
in the estimation of Kriging
parameters (Sections 1.6 and 2.4.5)
 .Regression Table 16 Options to specify the noise in the
model responses for a Gaussian
process regression model
(Sections 1.2.2 and 2.5)
 .Kriging Table 18 Field to specify the parameters of a
custom Kriging metamodel
(Section 2.9)
 .KeepCache Logical Flag to keep (cached) some
default: true important matrices. Keeping the
cached matrices may speed up the
Kriging predictor calculation for
some problems. If set to false,
some cached matrices are removed
after the metamodel has been
created.
 .Scaling Logical or INPUT object Conduct Kriging calculation in an
default: true auxiliary space, either by
standardizing the experimental
design or using an INPUT object. This
affects the value of various Kriging
parameters, notably the correlation
function hyperparameters
(Section 2.10).
Logical Flag to scale the experimental
design. If set to true, the
experimental design will be scaled
(specifically, standardized) before
calculating the Kriging metamodel.
INPUT object INPUT object that defines the
auxiliary probabilistic space
 .ValidationSet Table 19 Independent validation data set
(Section 3.1.8)

UQL AB-V2.1-105 - 48 -
 Kriging (Gaussian process modeling)

 .isTrajectory Logical Enable/Disable the generation of

default: false sample paths with the Kriging model
(Sections 1.8 and 2.11)
true Create the underlying random field
discreitzation needed for sampling
trajectories
false Do not create the underlyging
random field discretization. Note
that in this case, it will not be
possible to generate sample paths
once the Kriging model is built.
 .GRF Structure Options for the Gaussian process
discretization (See UQL AB User
Manual – Random Fields). By
default:
Metaopts.GRF.DiscScheme =
'EOLE'
MetaOpts.GRF.EOLE.SPC = 10
MetaOpts.GRF.EnergyRatio =
0.995 MetaOpts.GRF.Domain =
[minX_ED maxX_ED] where
minX_ED = min (X ) and maxX_ED
= max (X ), with X being the
experimental design.

3.1.1 Experimental design options

Table 5: MetaOpts.ExpDesign
.Sampling String Type of sampling
default: 'LHS'
'MC' Monte Carlo sampling
'LHS' Latin Hypercube sampling
'Sobol' Sobol’ sequence sampling
'Halton' Halton sequence sampling
'Data' A sample directly obtained from an
external .mat file, defined in
.ExpDesign.DataFile.
'User' User-defined experimental design,
defined in .ExpDesign.X and
.ExpDesign.Y.
 .DataFile String Name of the .mat file that contains
the experimental design. It only
applies, and is required, if the type of
sampling is 'Data'.
 .NSamples Integer Number of sample points to draw. It
is required for all types of sampling,
except 'Data' and 'Code'.

UQL AB-V2.1-105 - 49 -
UQL AB user manual

 .X N × M Double User-defined experimental design

(model inputs) X . It only applies,
and is required, if the type of
sampling is 'User'.
 .Y N × Nout Double User-defined model responses Y. It
only applies, and is required, if the
type of sampling is 'User'.

3.1.2 Trend options

Table 6: MetaOpts.Trend
 .Type String Type of the Kriging trend
default: 'ordinary' (Sections 1.3 and 2.4.2)
'simple' Known trend function (i.e., simple
Kriging), defined in
.Trend.CustomF. No estimation of
the coefficients takes place.
'ordinary' Unknown constant trend (i.e.,
ordinary Kriging)
'linear' Linear polynomial trend
'quadratic' Quadratic polynomial trend
'polynomial' Polynomial trend of an arbitrary
degree, defined in .Trend.Degree.
'custom' User-defined trend (Section 2.4.2.2)
 .Degree Integer Degree of the polynomial trend. It
default: 0 only applies if the Kriging trend type
is 'polynomial'.
 .CustomF Double, String, or Cell Field to set up custom trend function
array
Double Constant trend in simple Kriging
String Function name that returns f (x), an
arbitrary trend function
Cell arrays of strings Respective names of functions
fi (x) , i = 1, . . . , P , such that each
function evaluated on X returns the
i’th column of F
Cell arrays of function Respective handles of functions
handles fi (x) , i = 1, . . . , P , such that each
function evaluated on X returns the
i’th column of F
 .TruncOptions Table 7 Polynomial basis truncation options

Table 7: MetaOpts.Trend.TruncOptions

UQL AB-V2.1-105 - 50 -
 Kriging (Gaussian process modeling)

 .qNorm Double Parameter for the q-norm (or

default: 1 hyperbolic) truncation scheme. For
details, see the UQL AB User Manual
– Polynomial Chaos Expansions.
 .MaxInteraction Integer Maximum interaction of the input
variables. This is used to limit the
basis terms to up to the given
MaxInteraction. For details, see
the UQL AB User Manual –
Polynomial Chaos Expansions.
 .Custom Double vector User-defined basis specified using
multi-indices

3.1.3 Correlation function options

Table 8: MetaOpts.Corr
 .Family String or function handle One-dimensional correlation family
default: 'Matern-5_2' (Section 1.4.1)
'Linear' Linear correlation (Eq. (1.28))
'Exponential' Exponential correlation (Eq. (1.29))
'Gaussian' Gaussian correlation (Eq. (1.30))
'Matern-3_2' Matérn-3/2 correlation (Eq. (1.32))
'Matern-5_2' Matérn-5/2 correlation (Eq. (1.33))
Function handle User-defined correlation family
 .Type String Type of the correlation function
default:
'Ellipsoidal'
'Ellipsoidal' Ellipsoidal correlation (Eq. (1.34))
'Separable' Separable correlation (Eq. (1.35))
 .Isotropic Logical Flag to determine whether the
default: false correlation function is isotropic or
anisotropic (Section 1.4.3)
 .Nugget Double scalar or vector Set of values that are added to the
default: 10−10 main diagonal of the correlation
matrix R (Section 1.4.4)
• If scalar, add this quantity to
the diagonal elements of R.
• If vector, add each element to
the corresponding diagonal
element of R.

UQL AB-V2.1-105 - 51 -
UQL AB user manual

 .Handle Function handle Handle of the function that is used to

default: calculate the correlation matrix R
@uq_eval_Kernel (Section 2.4.3.2)

3.1.4 Estimation method options

Table 9: MetaOpts.CV
 .LeaveKOut Integer (< N ) Number of sample points left out in
default: 1 cross-validation (CV) (Section 2.4.4)
• The number of sample points
in Dk (Eq. (1.53)). It can be
any integer between 1 and N
(the size of experimental
design).
• If the total number of sample
points is not divisible by the
number of sample points left
out, the last subset contains
the remainder.
• The default is 1, also known as
leave-one-out (LOO) CV.

3.1.5 Hyperparameters optimization options

Table 10: MetaOpts.Optim

 .InitialValue Scalar or M × 1 Double Initial estimate of the correlation
default: 1.0 parameters
Scalar Double Initial estimate of the correlation
parameter. If M > 1, the same value
is assigned for all input dimensions.
M × 1 Double Initial estimates of the correlation
parameters for each of the M input
dimensions
 .Bounds 2 × M or 2 × 1 Double Bound of admissible values of the
default: [10−3 , 10]T correlation parameters, in the form
[lower_bound; upper_bound].
2 × 1 Double Bound of the correlation parameter.
If M > 1, the same value is assigned
for all input dimensions.
2 × M Double Bounds of the correlation parameters
for each of the M input dimensions
 .Display String Option to determine the verbosity of
default: 'none' the optimization process

UQL AB-V2.1-105 - 52 -
 Kriging (Gaussian process modeling)

'none' Nothing will be printed to the

command window
'final' Only the final result of the
optimization process will be printed
'iter' State of the optimization process will
be printed after each iteration
 .MaxIter Integer Maximum number of iterations or
default: 20 generations
 .Tol Double Covergence tolerance of the
default: 10−4 optimization

 .Method String Optimization method (Section 2.4.5)

default: 'HGA'
'none' No optimization. The value
.Optim.InitialValue is used.
'LBFGS' Gradient-based optimization
(L-BFGS) using the fmincon
function of M ATLAB
'GA' Genetic Algorithm (GA) optimization
using the ga function of M ATLAB
'HGA' Hybrid Genetic Algorithm (HGA)
optimization using the functions ga
and fmincon of M ATLAB
'CMAES' Covariance Matrix
Adaptation-Evolution Strategy
(CMA-ES) optimization using the
function uq_cmaes of UQL IB
'HCMAES' Hybrid Covariance Matrix
Adaptation-Evolution Strategy
(HCMA-ES) optimization using the
functions uq_cmaes of UQL IB and
fmincon of M ATLAB

.BFGS Table 11 Options relevant to the L-BFGS
optimization method

.GA Table 12 Options relevant to the GA
optimization method

.HGA Table 13 Options relevant to the HGA
optimization method

.CMAES Table 14 Options relevant to the CMA-ES
optimization method

.HCMAES Table 15 Options relevant to the HCMA-ES
optimization method

UQL AB-V2.1-105 - 53 -
UQL AB user manual

Note: The convergence tolerance defined under Metaopts.Optim.Tol may slightly dif-
fer depending on the optimization method. For M ATLAB optimization functions
(fmincon and ga), check the definition of the option 'TolFun' of the respective
function.

Note: Some of the above optimization methods require special optimization toolboxes
licenses in M ATLAB. An error is issued when the user specifies an algorithm for
which the corresponding toolbox is not available. The default optimizer, i.e.,
'HGA', is then replaced by:

– 'GA', if the Optimization Toolbox is not available;

– 'HCMAES', if the Global Optimization Toolbox is not available;
– 'CMAES', if both the Global Optimization and the Optimization toolboxes
are not available.

Table 11: MetaOpts.Optim.BFGS

 .nLM Integer (≥ 1) Limited memory size of the L-BFGS
default: 5 method (.nLM ≥ 1)

Table 12: MetaOpts.Optim.GA

 .nPop Integer Population size of each generation
default: 30
 .nStall Integer Maximum number of stall
default: 5 generations

Table 13: MetaOpts:Optim.HGA

 .nPop Integer Population size of each generation
default: 30
 .nStall Integer Maximum number of stall
default: 5 generations
 .nLM Integer (≥ 1) Limited memory size of the L-BFGS
default: 5 method (.nLM ≥ 1). The method is
executed starting, as the initial value,
with the optimal value obtained from
the GA optimization.

Table 14: MetaOpts.Optim.CMAES

 .nPop Integer The population size of each
default: 30 generation
 .nStall Integer The maximum number of stall
default: 5 generations

Table 15: MetaOpts.Optim.HCMAES

UQL AB-V2.1-105 - 54 -
 Kriging (Gaussian process modeling)

 .nPop Integer Population size of each generation

default: 30
 .nStall Integer Maximum number of stall
default: 5 generations
 .nLM Integer Limited memory size of the L-BFGS
default: 5 method (.nLM ≥ 1). The method is
executed starting, as the initial value,
with the optimal value obtained from
the CMA-ES optimization.

3.1.6 Regression options

Table 16: MetaOpts.Regression

 .SigmaNSQ 'none', 'auto', Logi- Specification of the output noise
cal, or Double variance in a regression model
default: 'none' (Section 2.5)
'none' or false Flag to ignore the noise variance.
The resulting Kriging model will be
in interpolation mode.
'auto' or true Flag to estimate the homogeneous
(homoscedastic) noise variance
Scalar Double Homoscedastic noise variance σn2
N × 1 Double Independent heterogeneous
(heteroscedastic) noise variance σn2
N × N Double Noise covariance matrix Σn
 .SigmaSQ Table 17 Initial value and bound for the
optimization of the Gaussian process
variance σ 2 in the case of a
regression model with known noise
variance

Note: If the metamodel has more than one outputs, then MetaOpts.Regression is ex-
pected to be a structure array with length equal to the number of outputs and that
the parameters in Table 16 are defined for each MetaOpts.Regression(idx)
(where idx ranges from 1 to the total number of outputs). If only one regression
option is specified, then the option applies to all output by default.

Table 17: MetaOpts.Optim.Regression.SigmaSQ

 .InitialValue Double Initial value of the Gaussian process
default: 0.5Var [Y] variance
 .Bound 2 × 1 Double Lower and upper bounds for the
default: Gaussian process variance
[0.1 × Var [Y] ; 10 × Var [Y]]

UQL AB-V2.1-105 - 55 -
UQL AB user manual

where Var [Y] = N1 i (yi − ȳ)2 and ȳ = N1 i yi are the sample variance
P P
and mean of the observed responses, respectively.

3.1.7 Custom Kriging options

To create a user-defined Kriging model as described in Section 2.9, a MetaOpts.Kriging

options must be set according to Table 18.

Table 18: MetaOpts.Kriging

.Trend Table 6 Kriging trend definition
.beta Double Values of the trend coefficients β in
Eq. (2.2)
.sigmaSQ Double Value of σ 2 in Eq. (2.2)
.theta Double Values of θ in Eq. (2.2)
.Corr Table 8 Correlation function definition (See
below Note)
 .SigmaNSQ Scalar, Vector, or Matrix Values of the output noise variance
Double in a regression model (Sections 1.2.2
and 1.2.2)
Scalar Double Homoscedastic noise variance σn2
N × 1 Double Independent heteroscedastic noise
variance σn2
N × N Double Noise covariance matrix Σn

Note: In the current version of UQL AB, only the default correlation function handle
uq_eval_Kernel is supported. All other options follow Table 8.

Note: If the metamodel has more than one outputs, then MetaOpts.Kriging is ex-
pected to be a structure array with length equal to the number of outputs and the
parameters in Table 18 are defined for each MetaOpts.Kriging(idx) (where
idx ranges from 1 to the total number of outputs).

Note: To create a custom Kriging model, the fields Metaopts.ExpDesign.X and

Metaopts.ExpDesign.Y are also required (see Table 5).

3.1.8 Validation Set

If an independent validation set is provided, UQL AB automatically calculates the validation

error of the created Kriging model (Section 1.7.2). The options are set according to Table 19.

Table 19: MetaOpts.ValidationSet

UQL AB-V2.1-105 - 56 -
 Kriging (Gaussian process modeling)

.X N × M Double User-specified validation set

experimental design Xval
.Y N × Nout Double User-specified validation set response
Yval

UQL AB-V2.1-105 - 57 -
UQL AB user manual

3.2 Accessing the Results

Table 20: myKriging = uq_createModel(...)

.Name String Name of the Kriging metamodel
.Kriging Table 21 Kriging results
.Error Table 22 Error metrics of the metamodeling result
.ExpDesign Table 23 Options and values related to experimental
design
.Options Table 4 Options that were defined in MetaOpts
variable (Section 3.1)
.Internal Table 24 Internal fields

Table 21: myKriging.Kriging

.beta P × 1 Double Values of β(θ) in Eq. (2.2)
.sigmaSQ Double Value of σ 2 (θ) in Eq. (2.2)
.theta Scalar or Values of θ in Eq. (2.2)
1 × M Double
.sigmaNSQ Scalar, Vector, Value of noise variance in Gaussian process
or Matrix regression (Sections 1.2.2 and 1.2.2)
Double
Scalar Double Homoscedastic noise variance σn2
N × 1 Double Independent heteroscedastic noise variance σn2
N × N Double Noise covariance matrix Σn

Table 22: myKriging.Error

.LOO Double Leave-One-Out (LOO) error, calculated using
Eq. (1.61)
.Val Double Validation error (see Eq. (1.62) and
Section 2.7). Only available if a validation set
is provided (see Table 19).

Note: In general, the fields myKriging.Kriging and myKriging.Error are structure

arrays with length equal to the number of outputs of the Kriging model. To access
the results that correspond to the respective output, one should use, for example,
myKriging.Internal.Kriging(k).field where k = 1, ..., Nout .

Table 23: myKriging.ExpDesign

.NSamples Double Number of sample points
.Sampling String Sampling method

UQL AB-V2.1-105 - 58 -
 Kriging (Gaussian process modeling)

.X N × M Double Values in the experimental design

.U N × M Double Experimental design values in the auxiliary
(scaled) space. For details, see Section 2.10
.Y N × Nout Observed model responses (output) Y that
Double corresponds to the experimental design X

3.2.1 Internal fields (advanced)

Note: The internal fields of the Kriging metamodel object are not intended to be ac-
cessed or changed in most typical usage scenarios of the module. Note that some
internal fields are not documented in this user manual.

Table 24: myKriging.Internal

.Kriging Table 25 Internal fields with data related to the Kriging
model
.Runtime Structure Internal fields with variables that are used
during the calculation of the kriging
metamodel
.Error Table 26 Internal fields related to the metamodeling
error
.ExpDesign Table 27 Internal fields related to the experimental
design
.Scaling Logical or Scaling of the experimental design (see
INPUT object Section 2.10)
.Regression Table 28 Internal fields related to a Gaussian process
regression model (e.g., noise of the responses).
It is only available in the case of a regression
model

Table 25: myKriging.Internal.Kriging

.Trend Table 29 Internal fields related to the trend
.GP Table 30 Internal fields related to the Gaussian process
.Optim Table 31 Internal fields related to the optimization of
the hyperparameters
.sigmaNSQ Double Internal field that contains the value of noise
variance in Gaussian process regression (refer
to the entry .sigmaNSQ in Table 21)
.Cached Structure Internal fields related to cached variables

UQL AB-V2.1-105 - 59 -
UQL AB user manual

Note: The fields myKriging.Internal.Kriging and myKriging.Internal.Error

are, in general, structure arrays with length equal to the number of outputs of
the metamodel. To access the results that correspond to the respective output,
one should use, for example, myKriging.Internal.Kriging(k).field where
k = 1, ..., Nout .

Table 26: myKriging.Internal.Error

.LOOmean N × 1 Double Error between the observed data and the mean
of the Kriging predictor from leave-one-out
(LOO) cross-validation (CV) (i.e., from the
prediction made on a point conditioned on all
the other points in the observed data)
.LOOsd N × 1 Double Standard deviation of the Kriging predictor
from LOO CV
.varY 1 × Nout Double The variance of the observed model responses
(output) Y. This quantity can be multiplied
with the relative LOO error in
myKriging.Error.LOO to get the absolute
value of the LOO error (Eq. (1.61)).

Table 27: myKriging.Internal.ExpDesign

.muX 1 × M Double Mean of the inputs used in the scaling of the
experimental design (Section 2.10).
• If an experimental design is specified in
MetaOpts.ExpDesign, then this
corresponds to the empirical mean of
each input dimension.
• If INPUT object is specified in
MetaOpts.Input, then this
corresponds to the first moment of the
corresponding marginal of each input
dimension.

.stdX 1 × M Double Standard deviation of the inputs used in the

scaling of the experimental design
(Section 2.10).
• If an experimental design is specified in
MetaOpts.ExpDesign, then this
corresponds to the empirical standard
deviation of each input dimension.
• If INPUT object is specified in
MetaOpts.Input, then this
corresponds to the second moment of
the corresponding marginal of each
input dimension.

UQL AB-V2.1-105 - 60 -
 Kriging (Gaussian process modeling)

.varY 1 × Nout Double Variance of the observed model responses

(output) Y. This quantity can be multiplied
with the relative LOO error in
myKriging.Error.LOO to get the absolute
value of the LOO error (Eq. (1.61)).

Table 28: myKriging.Internal.Regression

.IsRegression Logical Flag that indicates the current Kriging model is
a Gaussian process regression model
.EstimNoise Logical Flag that indicates that the noise variance is
estimated
.Tau Structure Structure that contains the initial value and
bounds used in the optimization of the τ
parameter (Eq. (1.21))
.SigmaSQ Structure Structure that contains the initial value and
bounds used in the optimization of the σ 2
(Gaussian process variance)
.SigmaNSQ Scalar, Vector, Estimated or specified noise variance.
or Matrix
Double
.IsHomoscedastic Logical Flag that indicates the Gaussian process
regression model is homoscedastic, with a
constant noise variance.

Table 29: myKriging.Internal.Kriging.Trend

.F N × P Double Observation matrix, i.e., the basis functions of
the Kriging trend fj ’s evaluated on the
experimental design (Eq. (2.2)).
.beta p × 1 Double Values of β(θ)
b in Eq. (2.2)

Table 30: myKriging.Internal.Kriging.GP

.R N × N Double Correlation matrix of the Gaussian process on
the experimental design points (i.e.,
R = R(X , X ))
.sigmaSQ Double Value of σ 2 in Eq. (2.2)

Table 31: myKriging.Internal.Kriging.Optim

.Theta Scalar or Optimum value of θ
1 × M Double
.Tau Double Ratio of σn2 to σ 2 + σn2 as defined in Eq. (1.21).
Only available in a regression with unknown
homoscedastic noise variance.

UQL AB-V2.1-105 - 61 -
UQL AB user manual

.SigmaSQ Double Gaussian process variance, directly optimized,

in the case of known noise variance. Only
available in a regression with known noise
variance (Sections 1.5.2.3 and 1.5.1.3)
.ObjFun Double Objective function value at the optimum θ
.InitialObjFun Double Objective function value at the initial estimate
of θ
.nEval Double Number of objective function evaluations
during the optimization process
.nIter Double Number of iterations (or generations) during
the optimization process. For hybrid methods
(i.e., 'HGA', 'HCMAES'), only the number of
generations is taken into account.

UQL AB-V2.1-105 - 62 -
 Kriging (Gaussian process modeling)

3.3 Kriging predictor

Syntax

Y_mu = uq_evalModel(X)
Y_mu = uq_evalModel(myKriging,X)
[Y_mu,Y_sigma2] = uq_evalModel(...)
[Y_mu,Y_sigma2,Y_cov] = uq_evalModel(...)
Y_traj = uq_evalModel(X,R)
%Y t r a j = u q e v a l M o d e l ( X , ' evalTraj ' , true )
%Y t r a j = u q e v a l M o d e l ( X , ' randomSeed ' , s e e d )

Description

Y_mu = uq_evalModel(X) returns the mean of the Kriging predictor (N × Nout ) on the
points of X (N × M ) using the most recently created Kriging model.

Y_mu = uq_evalModel(myKriging,X) returns the mean of the Kriging predictor (N ×

Nout )on the points of N × M of X using the Kriging metamodel object myKriging.

[Y_mu,Y_sigma2] = uq_evalModel(...) additionally returns the variance of the Kriging

predictor (N × Nout ).

[Y_mu,Y_sigma2,Y_cov] = uq_evalModel(...) additionally returns the covariance ma-

trices of the Kriging predictor (N × N × Nout ), i.e., one covariance matrix per output
dimension. Note that the diagonal elements of Y_cov (in each output dimension) are
equal to the corresponding elements in Y_sigma2, i.e., Y_sigma2 = diag(Y_cov).

Y_traj = uq_evalModel(X,R) returns R trajectories of the currently active Kriging model.

Y_traj is a matrix of size N ×Nout ×R. R must be specified as a 1×1 double, containing
an integer value.

Note: By default, the most recently created model or metamodel is the current active
model.

UQL AB-V2.1-105 - 63 -
UQL AB user manual

3.4 Printing and visualizing a Kriging metamodel

UQL AB offers two commands to conveniently print reports containing contextually relevant
information for a given Kriging metamodel object.

3.4.1 Printing the results: uq_print

Syntax

uq_print(myKriging)
uq_print(myKriging,outIdx)
uq_print(myKriging, outIdx, 'Option1', 'Option2', ...)

Description

uq_print(myKriging) prints a report of the configuration options and results of the Krig-
ing metamodel object myKriging. If the model has multiple outputs, only the report
on the information about the first output dimension are printed.

uq_print(myKriging,outIdx) prints a report on the configuration options and results of

the Kriging metamodel object myKriging for the output dimensions specified in the
array outIdx.

uq_print(myKriging, outIdx, 'Option1', 'Option2', ...) prints a report only on

selected configuration options or results, the selection of which, are specified by the
string options 'Option1', 'Option2', etc.See Table 32 for the available options.

Table 32: uq_print options

'beta' Prints out the regression coefficients β
'theta' Prints out the final hyperparameters θ (either the optimal value
or the user-specified value if they are manually specified)
'F' Prints out the observation matrix F
'optim' Prints out hyperparameters optimization information (i.e.,
optimization methods and optimized hyperparameters)
'trend' Prints out trend-related information (e.g., the type, coefficients)
'GP' Prints out the Gaussian process-related information (e.g., the
correlation family, type)
'R' Prints out the correlation matrix R
'regression' Prints out the Gaussian process regression-related information

UQL AB-V2.1-105 - 64 -
 Kriging (Gaussian process modeling)

Examples

uq_print(myKriging,[1 3]) prints the Kriging metamodeling results for output dimen-
sions 1 and 3.

uq_print(myKriging, 3, 'theta', 'R') only prints the hyperparameters θ and corre-

lation matrix R for output dimension 3.

UQL AB-V2.1-105 - 65 -
UQL AB user manual

3.4.2 Visualize the results: uq_display

Syntax

uq_display(myKriging)
uq_display(myKriging, outIdx)
uq_display(myKriging, outIdx, 'R')

Description

uq_display(myKriging) creates a visualization of the Kriging metamodel object myKriging.

It plots the Kriging model predictions (i.e., the mean and standard deviation) against
the input. If the model has multiple outputs, only the prediction of the first output
dimension is plotted.

uq_display(myKriging,outIdx) plots the Kriging model predictions against the input for
the model output dimensions specified in the array outIdx.

uq_display(myKriging, outIdx, 'R') creates a display of the correlation matrix R of

the metamodel object myKriging for the model output dimensions specified in the
array outIdx.

Note: Creating plots of Kriging predictions against the model inputs using uq_display
is only available for Kriging model in 1- and 2-dimension. In the case of 1-
dimensional model, the function creates a line plot; while in the case of 2-
dimensional, the function creates a contour plot. Visualizing the correlation ma-
trix, however, can be used for arbitrary input dimensions.

Examples

uq_display(myKriging, [1 3]) creates two plots, one plot for each selected output di-
mension, of Kriging predictions against the input for the Kriging metamodel object
myKriging.

uq_display(myKriging, 1:3, 'R') displays the correlation matrices of the metamodel

object myKriging for output dimensions 1, 2, and 3 in three separate figures.

UQL AB-V2.1-105 - 66 -
References

Abramovitz, M. and Stegun, I. A. (1965). Handbook of mathematical functions. New York:

Dover Publications Inc. 9

Bachoc, F. (2013). Cross-validation and maximum likelihood estimations of hyper-

parameters of Gaussian processes with model misspecification. Computational Statistics
and Data Analysis, 66:55–69. 15

Byrd, R. H., Hribar, M. E., and Nocedal, J. (1999). An interior point algorithm for large scale
nonlinear programming. SIAM Journal on Optimization, 9(4):877–900. 16

Cressie, N. A. C. (1993). Statistics for spatial data. John Wiley & Sons Inc. 15

Dubourg, V. (2011). Adaptive surrogate models for reliability analysis and reliability-based
design optimization. PhD thesis, Université Blaise Pascal, Clermont-Ferrand, France. 3, 6,
11, 13

Ghanem, R. and Spanos, P.-D. (1991). Spectral stochastic finite-element formulation for
reliability analysis. J. Eng. Mech., 117(10):2351–2372. 19

Goldberg, D. E. (1989). Genetic Algorithms in Search, Optimization, and Machine Learning.

Addison-Wesley Professional. 16

Hansen, N. and Ostermeier, A. (2001). Completely derandomized self-adaptation in evolu-

tion strategies. Evolutionary Computation, 9(2):159–195. 17

Krige, D. G. (1951). A statistical approach to some mine valuation and allied problems on
the Witwatersrand. Master’s thesis, University of the Witwatersrand, South Africa. 1

Li, C. and Der Kiureghian, A. (1993). Optimal discretization of random fields. J. Eng. Mech.,
119(6):1136–1154. 19

Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(2):1246–1266. 1

Nocedal, J. (1980). Updating Quasi-Newton Matrices with Limited Storage. Mathematics of

Computation, 35(151):773–782. 16

Rasmussen, C. and Williams, C. (2006). Gaussian processes for machine learning. Adaptive
computation and machine learning. MIT Press, Cambridge, Massachusetts. 1, 5, 10, 11, 14

67
UQL AB user manual

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

Santner, T., Williams, B., and Notz, W. (2003). The design and analysis of computer experi-
ments. Springer series in Statistics. Springer. 1, 2, 3, 13, 15

UQL AB-V2.1-105 - 68 -