UQ[PY]LAB USER MANUAL - PYTHON

KRIGING (GAUSSIAN PROCESS MODELING)

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

CHAIR OF RISK, SAFETY AND UNCERTAINTY QUANTIFICATION

STEFANO-FRANSCINI-PLATZ 5

CH-8093 Z

URICH

Risk, Safety &

Uncertainty Quantification

How to cite UQ[PY]LAB

C. Lataniotis, S. Marelli, B. Sudret, Uncertainty Quantiﬁcation in the cloud with UQCloud, Proceedings of the 4th

International Conference on Uncertainty Quantiﬁcation in Computational Sciences and Engineering (UNCECOMP

2021), Athens, Greece, June 27–30, 2021.

How to cite this manual

C. Lataniotis, D. Wicaksono, S. Marelli, B. Sudret, UQ[py]Lab user manual - Python – Kriging (Gaussian process

modeling), Report UQ[py]Lab -V1.0-105, Chair of Risk, Safety and Uncertainty Quantiﬁcation, ETH Zurich,

Switzerland, 2024

BIBT

X entry

@TechReport{UQdoc_10_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 UQ[py]Lab - V1.0-105}

}

List of contributors:

Name Contribution

A. Giannoukou, A. Hlobilov

a Translation from the UQLab manual

Document Data Sheet

Document Ref. UQ[PY]LAB-V1.0-105

Title: UQ[PY]LAB user manual - Python – Kriging (Gaussian process mod-

eling)

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

Chair of Risk, Safety and Uncertainty Quantiﬁcation, ETH Zurich,

Switzerland

Date: 27/05/2024

Doc. Version Date Comments

V0.8 28/09/2022 Initial release

V1.0 27/05/2024 Minor update for UQ[PY]LAB release 1.0

Abstract

Kriging is a stochastic modeling algorithm that has many applications in most ﬁelds of engi-

neering and applied mathematics. The UQ[PY]LAB metamodeling tool provides an efﬁcient,

modular, and easy to use Kriging module that allows one to obtain efﬁcient Kriging predictors

with minimal effort. For experienced users, numerous more advanced conﬁguration 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: UQCloud, UQ[PY]LAB, Kriging, Gaussian process metamodeling, Gaussian pro-

cess 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

2 Usage 19

2.1 Reference problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

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

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

2.4 Kriging metamodel setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.4.1 Speciﬁcation and generation of experimental design . . . . . . . . . . 24

2.4.2 Trend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.4.3 Correlation function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.4.4 Estimation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.4.5 Optimization methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

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

2.5.1 Unknown homoscedastic noise . . . . . . . . . . . . . . . . . . . . . . 33

2.5.2 Known noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

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

2.6.1 Accessing the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2.7 Using a validation set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

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

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

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

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

2.12 Using Kriging with constant inputs . . . . . . . . . . . . . . . . . . . . . . . . 41

3 Reference List 42

3.1 Create a Kriging metamodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.1.1 Experimental design options . . . . . . . . . . . . . . . . . . . . . . . . 46

3.1.2 Trend options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.1.3 Correlation function options . . . . . . . . . . . . . . . . . . . . . . . . 48

3.1.4 Estimation method options . . . . . . . . . . . . . . . . . . . . . . . . 49

3.1.5 Hyperparameters optimization options . . . . . . . . . . . . . . . . . . 49

3.1.6 Regression options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.1.7 Custom Kriging options . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.1.8 Validation Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.2 Accessing the Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.2.1 Internal ﬁelds (advanced) . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.3 Kriging predictor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

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

3.4.1 Printing the results: uq.print . . . . . . . . . . . . . . . . . . . . . . . 61

3.4.2 Visualize the results: uq.display . . . . . . . . . . . . . . . . . . . . . 63

References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

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 ﬁrst 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-

ﬁcation, 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 classiﬁcation tool.

This user manual presents the Kriging metamodeling tool of UQ[PY]LAB that supports Krig-

ing for both interpolation and regression. This chapter brieﬂy 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 UQ[PY]LAB. Finally, Section 3

provides a comprehensive list of the available commands and options.

UQ[PY]LAB user manual - Python

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 ∈ D

⊂ R

. A

Kriging metamodel is described by the following equation (Santner et al., 2003)

(x) = β

f (x) + σ

Z (x, ω) . (1.1)

The ﬁrst term in Eq. (1.1), β

f(x), is the mean value of the Gaussian process (i.e., its trend);

it consists of P arbitrary functions {f

; j = 1, . . . , P } and the corresponding coefﬁcients

{β

; j = 1, . . . , P }. The second term in Eq. (1.1) consists of the (constant) variance of the

Gaussian process σ

and a zero-mean, unit-variance, stationary Gaussian process Z(x, ω).

The underlying probability space is represented by ω and is deﬁned in terms of a correlation

function R (i.e., correlation family) and its hyperparameters θ. The correlation function

R = R(x, x

′

; θ), in turn, describes the correlation between two sample points in the output

space, that depends on x, x

′

, and the hyperparameters θ.

1.2.1 Prediction with noise-free responses (interpolation)

In the context of metamodeling, it is of interest to predict M

(x) for a new point x, given

the (observed) experimental design X =



(1)

, . . . , x

(N)



and the corresponding noise-free

model responses Y =



(1)

= M



(1)



, . . . , y

(N)

= M



(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.,

Y (x) and

the true model responses Y, has a joint Gaussian distribution deﬁned by



Y (x)



∼ N

N+1



(x) β

F β



, σ



1 r

(x)

r (x) R



(1.2)

where:

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

= f



(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

= R



x, x

(i)

; θ



, i = 1, . . . , N. (1.4)

• R is the correlation matrix with elements:

= R



(i)

, x

(j)

; θ



, i, j = 1, . . . , N. (1.5)

UQ[PY]LAB-V1.0-105 - 2 -

Kriging (Gaussian process modeling)

Consequently, the mean and variance of the Gaussian random variate

Y (x) conditional on

the observed data X and Y can be computed as follows (Santner et al., 2003; Dubourg, 2011)

(x) = f (x)

β + r (x)

−1



Y − F



, (1.6)

(x) = σ



1 − r

(x) R

−1

r (x) + u

(x)



−1



−1

u (x)



, (1.7)

where

β =



−1



−1

Y (1.8)

is the generalized least-squares estimate of β and

u (x) = F

−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

Y (x) ∼ N



(x) , σ

(x)



(1.10)

and therefore



Y (x) ≤ t



= Φ



t − µ

(x)



, (1.11)

where Φ(·) denotes the Gaussian cumulative density function. Note that this is the probability

that

Y (x) is smaller than t at a given x

. Based on Eq. (1.11), the conﬁdence intervals on the

predictor can be calculated by

Y (x) ∈

(x) − Φ

−1



1 −



(x) , µ

(x) + Φ

−1



1 −



(x)

(1.12)

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, Σ

) , (1.14)

where Σ

is the covariance matrix of the noise term.

This shall not be confused with a probability of failure computed in the framework of rare event simulation.

For details, see UQ[PY]LAB User Manual – Structural Reliability.

UQ[PY]LAB-V1.0-105 - 3 -

UQ[PY]LAB user manual - Python

Depending on the properties of Σ

, three classes of noise can be identiﬁed:

• Σ

= σ

I (where I is an identity matrix) corresponds to homogeneous (homoscedastic)

noise, in which the variance σ

is constant and the same for all observed responses. In

other words, the noise is an independent and identically distributed (iid) Gaussian.

• Σ

= diag





corresponds to independent heterogeneous (heteroscedastic) noise, in

which the noise variances σ

differ for each observed response but are not correlated.

• Σ

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.,

Y (x) and the observed,

albeit noisy, model responses Y now becomes:



Y (x)



∼ N

N+1



(x) β

F β





(x)

r (x) σ

R + Σ



. (1.15)

Because in general the noise variance cannot be factored out from the covariance matrix (i.e.,

due to heteroscedasticity), the covariance matrix C = σ

R + Σ

is used in the subsequent

formulation. Consequently, the Kriging mean and variance predictors (see Eq (1.6) and (1.7))

become

(x) = f (x)

β + c(x)

−1



Y − F



, (1.16)

(x) =



− c

(x) C

−1

c(x) + u

(x)(F

−1

F )

−1

(x)



(1.17)

where c = σ

r (x) is the cross-covariance vector,

β =



−1



−1

Y, (1.18)

is the generalized least-square estimate of the regression coefﬁcients

β, and

(x) = F

−1

c(x) − f (x). (1.19)

In the special case of homoscedastic noise, C = σ

R + σ

I, thus with

total

= σ

+ σ

(1.20)

and

τ =

total

(1.21)

Eq. (1.15) can be written as



Y (x)



∼ N

N+1



(x) β

F β



, σ

total



(1 − τ )

(x)

r (x)



(1.22)

where

r = (1 − τ) r (1.23)

UQ[PY]LAB-V1.0-105 - 4 -

Kriging (Gaussian process modeling)

and

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

In this case, the mean and variance of the Gaussian random variate

Y (x) can be computed

as follows (Rasmussen and Williams, 2006)

(x) = f (x)

β +

r(x)

−1



Y − F



(1.25)

(x) = σ

total



1 −

(x)

−1

r(x) + u

(x)(F

−1

F )

−1

u(x)



(1.26)

where

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

and

r, 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, x

′

; θ). 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 σ

(if it is unknown) or the

Gaussian process variance σ

(if σ

is known) need to be estimated.

• Using the optimal value of θ, calculate the rest of the unknown Kriging parameters

(e.g., β) and, if necessary, σ

Predictions at new points can then be made in terms of the mean and variance of

Y (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 ﬁts the noisy responses. Using Eq. (1.12), the 95% conﬁdence intervals

of both predictors are also calculated and shown in the plot.

UQ[PY]LAB-V1.0-105 - 5 -

UQ[PY]LAB user manual - Python

(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 β

f(x) in Eq. (1.1).

While using a non-zero trend is optional, it is commonly preferred

. 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

f(x) =

j=1

(x),

where f

’s are arbitrary but fully speciﬁed functions. Note that no estimation on β is

taken place and the coefﬁcients β are all 1’s

• Ordinary Kriging

In ordinary Kriging, the trend has a constant yet unknown value

f(x) = β

(x) = β

where by convention f

(x) = 1.

• Universal Kriging

Universal Kriging, the most general and ﬂexible formulation, assumes that the trend is

Note that the mean of the Kriging predictor given in Eq. (1.6), (1.16), or (1.25) is not conﬁned to be zero

even though the trend is zero.

UQ[PY]LAB-V1.0-105 - 6 -

Kriging (Gaussian process modeling)

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

f(x) =

j=0

(x).

According to the above, the simple and ordinary Kriging are special cases of universal Kriging.

In the current version of UQ[PY]LAB, the functions f

’s can either be constants, polynomials,

or any other user-deﬁned 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) β

linear β

i=1

quadratic β

i=1

+ Σ

i=1

j=1

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

f(x) =

α∈A

M,P

(x) , f

(x) =

i=1

(1.27)

where α = {α

, . . . , α

} is a vector of indices and A

M,P



α ∈ N

: |α| ≤ P



denotes the

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

The Kriging module in UQ[PY]LAB offers the possibility to use a multivariate polynomial of

arbitrary 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 covariance

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, x

′

) is, in general, not a valid correlation function. The necessary

conditions for a function R (x, x

′

) to be a correlation function are:

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

UQ[PY]LAB-V1.0-105 - 7 -

UQ[PY]LAB user manual - Python

• The correlation matrix

whose elements are computed as R

= R



(i)

, x

(j)



, (x

, x

) ∈

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

mental design X.

• The function is symmetric, i.e., R(x, x

′

) = R(x

′

, x), ∀x

′

, x ∈ D

In general, a correlation function can be expressed in the form R(x, x

′

; θ) where θ is a vector

that contains a set of parameters. Typically θ ∈ R

, 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 deﬁned for a pair of one-dimensional input x, x

′

∈ R and parametrized by θ ∈ R

, often

referred to as the characteristic length scale or simply the scale parameter.

For each of the available correlation families in UQ[PY]LAB, the function is plotted in Fig-

ures 2-6 with different scale values, along with the corresponding sample paths (i.e., realiza-

tions) of a zero-mean, unit-variance Gaussian process having this correlation function.

• Linear:



x, x

′

; θ



= max



0, 1 −

|x − x

′



. (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, x

′

; θ



= exp



−

|x − x

′



. (1.29)

The resulting sample paths are C

, i.e., continuous but non-differentiable.

It is also a Gram matrix.

UQ[PY]LAB-V1.0-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):



x, x

′

; θ



= exp

−



|x − x

′



. (1.30)

The resulting sample paths of the corresponding process are inﬁnitely 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

ern: The general form of the Mat

ern correlation function is given by



x, x

′

; θ, v



v−1

Γ(v)



√

|x − x

′





√

|x − x

′



(1.31)

where θ is the correlation function scale parameter, v ≥ 1/2 is the shape parameter, Γ

is the Euler’s Gamma function, and K

is the modiﬁed Bessel function of the second

kind

An interesting feature of this correlation function family is that the sample paths of the

corresponding Gaussian process are ⌈v − 1⌉ times differentiable, where ⌈·⌉ denotes the

ceiling function. For v = 1/2, Mat

ern kernel coincides with the exponential correlation

It is also also known as the Bessel function of the third kind. For details, see Abramovitz and Stegun (1965).

UQ[PY]LAB-V1.0-105 - 9 -

UQ[PY]LAB user manual - Python

function which generates C

sample paths. For v → ∞, it tends towards the Gaussian

correlation function which generates C

∞

sample paths.

The Mat

ern functions can be computed using simpliﬁed formulas when v = p + 1/2,

where p is a non-negative integer (Rasmussen and Williams, 2006). The Mat

ern func-

tions with v = 3/2 and 5/2 (abbreviated as Mat

ern-3/2 and Mat

ern-5/2, respectively)

are the most commonly used and available in UQ[PY]LAB. Their formulas are given

below.

For v = 3/2

R (h; θ, v = 3/2) =



1 +

√

|x − x

′



exp



−

√

|x − x

′



, (1.32)

Figure 5: The Mat

ern 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

R(x, x

′

; θ, v = 5/2) =

1 +

√

|x − x

′



|x − x

′



exp



−

√

|x − x

′



. (1.33)

Figure 6: The Mat

ern 5/2 correlation function (Eq. (1.33)) and sample paths drawn from the

corresponding zero-mean unit-variance Gaussian process for various scale parameters.

UQ[PY]LAB-V1.0-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



x, x

′

; θ



= R (h) , h =

i=1



− x

′



0.5

. (1.34)

• Separable correlation functions (Sacks et al., 1989; Dubourg, 2011), calculated as fol-

lows



x, x

′

; θ



i=1



, x

′

, θ



. (1.35)

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-

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 deﬁned as follows:

• Isotropic ellipsoidal correlation function given as



x, x

′

; θ



= R (h) ; h =

i=1



− x

′



0.5

. (1.36)

• Isotropic separable correlation function given as



x, x

′

; θ



i=1



, x

′

; θ



. (1.37)

UQ[PY]LAB-V1.0-105 - 11 -

UQ[PY]LAB user manual - Python

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 x

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

= 1 + ν

. (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 UQ[PY]LAB are discussed in the following subsections.

1.5.1 Maximum-likelihood estimation

The idea behind the maximum-likelihood (ML) estimation method is to ﬁnd the set of Krig-

ing parameters β, σ

, θ, and if apply, σ

, such that the likelihood of the observations Y =





(1)



, . . . , M



(N)



is maximized. Since Y is assumed to follow a multivariate Gaus-

sian distribution (recall basic Kriging assumptions), the likelihood function L



β, σ

, θ; Y



reads



β, σ

, θ; Y



(det C)

−1/2

(2π)

N/2

exp



−

(Y − F β)

−1

(Y − F β)



(1.39)

where the covariance matrix C sums up the covariance matrix of the Gaussian process σ

and covariance matrix of the noise of the responses Σ

as follows

C = σ

R + Σ

. (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 σ

R, and the

likelihood function can be written as



β, σ

, θ; Y



(det R)

−1/2

(2πσ

)

N/2

exp



−

2σ

(Y − F β)

−1

(Y − F β)



. (1.41)

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

UQ[PY]LAB-V1.0-105 - 12 -

Kriging (Gaussian process modeling)

β and σ

that are strictly functions of θ are obtained (for proof and more details, refer to

Dubourg (2011); Santner et al. (2003))

β = β (θ) =



−1



−1

Y, (1.42)

bσ

= σ

(θ) =

(Y − F β)

−1

(Y − F β) . (1.43)

The hyperparameters θ, in turn, are obtained from solving the following optimization prob-

lem

θ = arg min

θ∈D

[−log L(θ; Y)] . (1.44)

Based on Eqs. (1.41) to (1.43), the optimization problem of Eq. (1.44) can be written as

θ = arg min

θ∈D



log (det R) + N log



2πσ



+ N



. (1.45)

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 σ

total

R, where σ

total

= σ

+ σ

R = (1 − τ ) R + τI, and

τ =

total

(Section 1.2.2). The corresponding likelihood then reads



β, σ

, θ, τ; Y





det



−1/2



2πσ

total



N/2

exp



−

2σ

total

(Y − F β)

−1

(Y − F β)



, (1.46)

where σ

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

total

and

R, respectively.

Similarly, the estimates for β and σ

total

are written as

β = β (θ) =



−1



−1

Y, (1.47)

bσ

total

= σ

total

(θ) =

(Y − F β)

−1

(Y − F β) . (1.48)

Based on Eqs. (1.46) to (1.48), the optimization problem can be formulated as follows

θ, bτ = arg min

θ∈D

,τ∈(0,1)

log



det



+ N log



2πbσ

total



+ N

. (1.49)

Notice that compared to Eq. (1.45), there is an additional parameter τ to be optimized and

that it is bounded in (0, 1).

The logarithm of the likelihood is usually taken in the optimization to avoid numerical underﬂow problem.

UQ[PY]LAB-V1.0-105 - 13 -

UQ[PY]LAB user manual - Python

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



−1

Y. (1.50)

There is, however, no analytical expression for the Gaussian process variance σ

. and it

must be simultaneously optimized with the rest of the hyperparameters θ. The optimization

problem can be formulated as follows:

θ, bσ

= arg min

θ∈D

,σ

∈D



log (det C) + N (log(2π) + N +



Y − F



−1



Y − F





(1.51)

Note that, if the known noise is independent (i.e., a vector σ

), the noise covariance matrix

is simpliﬁed to diag





; whereas, if the known noise variance is homoscedastic (i.e., a

scalar σ

), the noise covariance matrix Σ

is simpliﬁed to σ

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 =



(i)

, y

(i)



, i = 1, . . . , N



into K mutually

exclusive and collectively exhaustive subsets {D

, k = 1, . . . , K} such that

∩ D

= ∅ , ∀(i, j) ∈ {1, . . . , K}

and

[

k=1

= D. (1.52)

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 speciﬁc 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

CV,k

(

(i)

)

∈D



(i)

− µ



(i)

; β, σ

, θ, σ

, D \ D



(

(i)

)

∈D



(i)

− µ

Y ,\D



(i)

; β, σ

, θ, σ



(1.53)

where µ

Y ,\D

denotes the mean of the Kriging predictor (see, for instance, Eq. (1.6) for the

noise-free prediction) on the k-th CV subset



(i)

, y

(i)



∈ D

conditioned on the other K −1

UQ[PY]LAB-V1.0-105 - 14 -

Kriging (Gaussian process modeling)

subsets. In the case of LOO CV, the number of elements in D

is equal to one.

The overall cross-validation error then reads



β, σ

, θ, σ

; Y



k=1

CV,k

. (1.54)

The idea behind CV estimation method is to ﬁnd the set of Kriging parameters β, σ

, θ, and,

if apply, σ

, 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)

θ = arg min

θ∈D

(θ; Y) (1.55)

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 σ

from the predictor

mean (Eq. (1.6)), and the absence of σ

Using

θ, the estimate of β is computed as in Eq. (1.8), while the estimate of σ

is computed

using the following equation (Cressie, 1993; Bachoc, 2013)

bσ

= σ

(

θ) =

k=1

i∈D



(i)

− µ

Y ,\D



(i)

;



Y ,\D



(i)

;



(1.56)

where c

Y ,\D

denotes the normalized variance of the Kriging predictor on the k-th CV subset



(i)

, y

(i)



∈ D

conditional on all points of the other K − 1 subsets. In other words,

Y ,\D

(1.57)

where σ

Y ,\D

denotes the variance of the Kriging predictor (Eq. (1.7)) on the same k-th CV

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τ = arg min

θ∈D

,τ∈(0,1)

(θ, τ; Y) . (1.58)

UQ[PY]LAB-V1.0-105 - 15 -

UQ[PY]LAB user manual - Python

in which µ

Y ,\D



(i)

; θ, τ



given in Eq. (1.25) is used as the Kriging predictor.

Using

θ and bτ, the estimate of β is computed as in Eq. (1.47), while the estimate of σ

total

computed using the following equation

bσ

total

= σ

total



θ, bτ



k=1

i∈D



(i)

− µ

Y ,\D



(i)

;

θ, bτ



Y ,\D



(i)

;

θ, bτ



(1.59)

where c

Y ,\D

above is analogous to the one appeared in Eq. (1.56) but based on Eq. (1.26),

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 σ

is optimized with rest of

the hyperparameters θ. The optimization problem now reads

θ, bσ

= arg min

θ∈D

,σ

∈D



θ, σ

; Y, Σ



(1.60)

where Σ

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 Σ

is simpliﬁed to diag





, where

is the vector of noise variances; while if the noise is homoscedastic (i.e., a scalar σ

), the

noise covariance matrix Σ

is simpliﬁed to σ

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 ﬂat 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 brieﬂy 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 ﬁelds for the past decades. A hybrid approach can also be used, in

UQ[PY]LAB-V1.0-105 - 16 -

Kriging (Gaussian process modeling)

which the ﬁnal 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

Figure 7. The original model is M(x) =



1 + 25x



−1

(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

LOO







i=1



M(x

) − µ

Y ,(−i)

)



Var [Y]







(1.61)

where µ

Y ,(−i)

) denotes the Kriging metamodel that is obtained using all the points of the

experimental design X, except x

1.7.2 Validation error

The validation error is calculated as the relative generalization error on an independent set of

data D

val

n

(i)

val

, y

(i)

val



, i = 1, . . . , N

val

as follows

val

− 1

val







val

i=1



M(x

(i)

val

− M



(i)

val



val

i=1



(i)

val



− ˆµ

val









(1.62)

UQ[PY]LAB-V1.0-105 - 17 -

UQ[PY]LAB user manual - Python

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

(b) Maximum likelihood (ML) 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.

where ˆµ

val

i=1



(i)

val



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.

UQ[PY]LAB-V1.0-105 - 18 -

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 UQ[PY]LAB.

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",

"MetaType": "Kriging"

}

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

(x) =

j=1





(x) + σ



x ; R





, (2.2)

where

θ is obtained by solving an optimization problem

θ = arg min

θ∈D

J(θ) (2.3)

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

UQ[PY]LAB user manual - Python

to obtain a Kriging metamodel are the following:

• An experimental design, X, and the corresponding model responses, Y. If they are not

available, they can be generated by deﬁning a full computational (true) model and a

probabilistic input model

• A trend speciﬁcation, i.e., selection of the type of trend component functions f

(x) , j =

1, . . . , P .

• An appropriate correlation function R(x, x

′

; θ).

• A hyperparameter estimation method that deﬁnes the objective function J(θ) in Eq. (2.3);

this choice also deﬁnes the approach to estimate the Gaussian process variance σ

• 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 conﬁguration example to create a Kriging metamodel in UQ[PY]LAB is

given as was discussed in Section 2.2.

from uqpylab import sessions

import numpy as np

# Start the session

mySession = sessions.cloud(host=UQCloudhost,

token=UQCloudtoken)

# (Optional) Get a convenient handle to the command line interface

uq = mySession.cli

# Reset the session

mySession.reset()

uq.rng(100,'twister') # For reproducible results

# Create experimental design

X = np.arange(0, 15, 2)

Y = X

np.sin(X)

# Define a Kriging metamodel

MetaOpts = {

'Type': 'Metamodel',

'MetaType': 'Kriging',

'ExpDesign':{

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

Note that this step belongs to the general context of metamodeling, i.e., it is not speciﬁc to Kriging (see, for

instance, UQ[PY]LAB User Manual – Polynomial Chaos Expansions)

UQ[PY]LAB-V1.0-105 - 20 -

Kriging (Gaussian process modeling)

}

# Create the Kriging metamodel

myKriging = uq.createModel(MetaOpts)

# Terminate the remote UQCloud session

mySession.quit()

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 ---------------%

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 ﬁgure produced by uq.display is shown in Figure 8.

UQ[PY]LAB-V1.0-105 - 21 -

UQ[PY]LAB user manual - Python

Note: Note that uq.display for Kriging MODEL objects can only be used for model

one- and two-dimensional inputs.

0 5 10

−10

−5

Kriging approximation

95% confidence interval

Observations

𝑋

𝑌

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

2.3.1 Accessing the results

The Kriging MODEL object myKriging created in Section 2.3 uses all the default conﬁgura-

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.

2.3.1.1 Kriging parameters

The values of the basic parameters of the Kriging metamodel, namely θ, σ

, and β are con-

tained in the dictionary myKriging['Kriging']:

{

'beta': 31.66474977533767,

'sigmaSQ': 118181.5463011093,

'theta': 2.905772037286146

}

2.3.1.2 Model evaluations

The experimental design and the corresponding model responses are stored in the

myKriging['ExpDesign'] dictionary:

{

'Sampling': 'User',

'X': [0, 2, 4, 6, 8, 10, 12, 14],

'Y': [0, 1.8185948536513634, -3.027209981231713, -1.6764929891935552,

7.914865972987054, -5.440211108893697, -6.43887501600522,

13.868502979728184],

'NSamples': 8,

UQ[PY]LAB-V1.0-105 - 22 -

Kriging (Gaussian process modeling)

'U': [-1.4288690166235207, -1.0206207261596576, -0.6123724356957946,

-0.20412414523193154, 0.20412414523193154, 0.6123724356957946,

1.0206207261596576, 1.4288690166235207]

}

The ﬁelds 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']:

{

'LOO': 0.5555163005242105

}

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 conﬁguration 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 speciﬁ-

cation of the experimental design. The way to use each option is discussed in Sec-

tion 2.4.1 and list of all available options can be found in Table 5.

• MetaOpts['Trend'] contains the options regarding the term

j=1

(θ)f

(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 hy-

perparameters θ (maximum-likelihood or cross-validation). The choice of estimation

method corresponds to different ways of calculating σ

(θ) as well as J(θ). See Sec-

tion 1.5 for a theoretical introduction of each estimation method.

UQ[PY]LAB-V1.0-105 - 23 -

UQ[PY]LAB user manual - Python

• MetaOpts['Optim'] contains the options related to the method for solving the opti-

mization problem in Eq. (2.3). See Section 1.6 for a brief description of the available

optimization 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.

2.4.1 Speciﬁcation and generation of experimental design

An experimental design for constructing a Kriging metamodel can either be speciﬁed using

existing data or generated from a MODEL and an INPUT objects.

2.4.1.1 Using existing data

If data is stored in the variables X, Y:

MetaOpts['ExpDesign'] = {

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

}

The input and output dimensions of the problem M, N

out

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 an INPUT object, it ﬁrst needs to be created. For the reference problem in

Eq. (2.1) this reads:

InputOpts = {

'Marginals': [{

'Type': 'Uniform',

'Parameters': [0, 15]

}]

}

myInput = uq.createInput(InputOpts)

The input and output dimensions of the problem are automatically inferred from the con-

ﬁguration of the INPUT object myInput and MODEL object myModel . For details about the

conﬁguration options available for INPUT and MODEL objects, please refer to the UQ[PY]LAB

User Manual – the INPUT module and the UQ[PY]LAB User Manual – the MODEL module,

respectively.

Then, a MODEL object is created as follows:

ModelOpts = {

'Type' : 'Model',

'mString' : 'X.

sin(X)'

}

myModel = uq.createModel(ModelOpts)

UQ[PY]LAB-V1.0-105 - 24 -

Kriging (Gaussian process modeling)

Now, the INPUT and MODEL objects are speciﬁed into the Kriging metamodel conﬁguration:

MetaOpts['Input'] = myInput

MetaOpts['FullModel'] = myModel

Finally, the size of the experimental design is speciﬁed, 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 conﬁguration options are the following:

• When considering a constant (ordinary Kriging), linear, or quadratic trend, the ﬁeld

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',

'Degree': q

}

where q is an integer equal to the polynomial degree.

• For a constant trend with ﬁxed value (i.e., simple Kriging), set the following option:

MetaOpts['Trend'] = {

'Type': 'simple',

'Degree': V

}

where V is a real number.

A summary of the available trend options in UQ[PY]LAB along with the corresponding for-

mula 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 Required

'ordinary' β

Not required Not required

'linear' β

i=1

Not required Not required

UQ[PY]LAB-V1.0-105 - 25 -

UQ[PY]LAB user manual - Python

0 2 4 6 8 10 12 14

−10

0 2 4 6 8 10 12 14

Kriging, R: Ordinary Kriging

Kriging, R: 3rd deg. polynomial trend

Kriging, R: Custom trend

Observations

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 conﬁgura-

tions. In the top ﬁgure, black line indicates the true model prediction.

'quadratic' β

i=1

j=1

Not required Not required

'polynomial' see Eq. (1.27) Required Not required

'custom'

i=1

(x) Not required Required

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

ﬁgurations are plotted as in the example script uq_Example_Kriging_03_TrendTypes.

Note: The ﬁeld MetaOpts['Trend']['Type'] determines the trend type of a Krig-

ing metamodel. If the user does not deﬁne a trend type, the default

MetaOpts['Trend']['Type'] = 'ordinary' (i.e., ordinary Kriging) is used.

2.4.2.2 Advanced options

A user can deﬁne arbitrary custom trends via advanced options. First, the trend type is

speciﬁed as follows:

MetaOpts['Trend']['Type'] = 'custom'

UQ[PY]LAB-V1.0-105 - 26 -

Kriging (Gaussian process modeling)

Then the trend can be deﬁned as a string. If, for example, the trend is a sine function sin(x),

then the option is set as follows:

MetaOpts['Trend']['CustomF'] = 'sin'

UQCloud uses the MATLAB engine to interpret the function strings. Therefore, the function

can be any of the built-in MATLAB functions. The list of supported elementary math functions

can be found here.

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 equals to 1, M + 1,

(M+2)(M+1)

, respectively.

• If MetaOpts['Trend']['Type']='polynomial', then depending on the polynomial

degree q, deﬁned in MetaOpts['Trend']['Degree'], P equals to



M+q



• If MetaOpts['Trend']['Type'] is 'custom', then P depends on how

MetaOpts['Trend']['CustomF'] is speciﬁed:

– If MetaOpts['Trend']['CustomF'] is a string or ﬂoat, then P equals to 1.

– If MetaOpts['Trend']['CustomF'] is a list of strings, then P equals the length

of the list. In other words, each element of the list 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 ﬂexibility 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 speciﬁes 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'

}

UQ[PY]LAB-V1.0-105 - 27 -

UQ[PY]LAB user manual - Python

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

'ellipsoidal'.

• Family that speciﬁes 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

and x

′

(see Eq. (1.34)). Each coordinate difference between the rows of x and x

′

is normalized by a positive scalar θ.

– For separable correlation functions, the correlation family is a function of the form

R(x, x

′

; θ) 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 ﬂag (with true or false value) that speciﬁes 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.

Note: In the current version of UQ[PY]LAB there is limited support for advanced cor-

relation functions options.

UQ[PY]LAB-V1.0-105 - 28 -

Kriging (Gaussian process modeling)

0 2 4 6 8 10 12 14

−20

−10

0 2 4 6 8 10 12 14

Kriging, R: Matern 5/2

Kriging, R: Linear

Kriging, R: Exponential

Observations

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 ﬁgure, black line indicates the true model prediction.

UQ[PY]LAB-V1.0-105 - 29 -

UQ[PY]LAB user manual - Python

0 2 4 6 8 10 12 14

−10

0 2 4 6 8 10 12 14

−10

0 2 4 6 8 10 12 14

True model

Kriging, optim.method: BFGS

Kriging, optim. method: HGA

Observations

0 2 4 6 8 10 12 14

True model

Kriging, optim.method: BFGS

Kriging, optim. method: HGA

Observations

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.

• 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 speciﬁed as a vector ν of length N, each of its elements is added to the

corresponding diagonal element of R, that is, R

= 1 + ν

, {i = 1, . . . , N}.

Note: By default, the nugget value in UQ[PY]LAB is set to ν = 10

−10

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 σ

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 sufﬁces to select either 'ML' or 'CV' in the ﬁeld 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,

UQ[PY]LAB uses the N -fold CV (or equivalently, the leave-one-out) approach.

UQ[PY]LAB-V1.0-105 - 30 -

Kriging (Gaussian process modeling)

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 = ⌈N/k⌉ (2.4)

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, UQ[PY]LAB 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 available for cross-validation. Otherwise UQ[PY]LAB will throw an

error.

2.4.5 Optimization methods

Various conﬁguration 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 θ

) is required. By default, θ

= 1 in each dimension. However, a different

initial value can be selected, e.g., θ

= 0.5 (in each dimension):

MetaOpts['Optim']['InitialValue'] = 0.5

Furthermore, different initial value per dimension can also be speciﬁed using an M ×1

vector instead of a scalar.

Note: Note that, by default, the hyperparameters are deﬁned in to the scaled

(auxiliary) space U as described in Section 2.10.

UQ[PY]LAB-V1.0-105 - 31 -

UQ[PY]LAB user manual - Python

• 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 UQ[PY]LAB. To use one of these methods, e.g., GA, the fol-

lowing option is set:

MetaOpts['Optim']['Method'] = 'GA'

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]

are used. To specify different bounds,

for example, [0.01, 1]

, the following options is set:

MetaOpts['Optim']['Bounds'] = [0.01, 1]

Note: Note that, by default, the hyperparameters are deﬁned in to the scaled

(auxiliary) space U as described in Section 2.10.

Additional options which are speciﬁc 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

}

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 ﬁnal 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, the Hybrid Genetic Algorithm (HGA) is

used by default.

If no optimization of the hyperparameters is required, due to, for example, a prescribed value

theta_user is used, the following options are set:

UQ[PY]LAB-V1.0-105 - 32 -

Kriging (Gaussian process modeling)

MetaOpts['Optim']['Method'] = 'none'

MetaOpts['Optim']['InitialValue'] = theta_user

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 UQ[PY]LAB

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), UQ[PY]LAB 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:

import numpy as np

import math

# Create experimental design with noisy responses

uq.rng(100,'twister') # For reproducible results

X = np.linspace(0,14,100)

noise_var = 3.0

Y = X

np.sin(X) + math.sqrt(noise_var)

np.random.randn(X.size)

Then, a Kriging metamodel as in Section 2.3 is speciﬁed, now with the additional option:

# Define a Kriging metamodel

MetaOpts = {

'Type': 'Metamodel',

'MetaType': 'Kriging',

'ExpDesign': {

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

'Regression': { # Estimate the noise

UQ[PY]LAB-V1.0-105 - 33 -

UQ[PY]LAB user manual - Python

'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.

%---------------- 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:

UQ[PY]LAB-V1.0-105 - 34 -

Kriging (Gaussian process modeling)

0 5 10

−15

−10

−5

Kriging approximation

95% confidence interval

Observations

𝑋

𝑌

Figure 12: The output of uq.display of a Kriging MODEL object having a one-dimensional

input, with noisy responses.

MetaOpts['Regression'] = {

'SigmaNSQ': constSigmaNSQ

}

where constSigmaNSQ is the known noise variance. The speciﬁcation 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 array. This is the case in

which the noises are independent, but may differ at each observation point.

• 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 hypothetical experimental design with noisy responses

uq.rng(100,'twister') # For reproducible results

X = np.linspace(0,14,50)

noise_var = 3.0

Y = X

np.sin(X) + math.sqrt(noise_var)

np.random.randn(X.size)

# Define a Kriging metamodel

MetaOpts = {

'Type': 'Metamodel',

'MetaType': 'Kriging',

'ExpDesign': {

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

UQ[PY]LAB-V1.0-105 - 35 -

UQ[PY]LAB user manual - Python

# Impose known homoscedastic noise variance

'Regression': {

'SigmaNSQ': noise_var

}

# Create the Kriging metamodel

myKriging = uq.createModel(MetaOpts)

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 hypothetical experimental design with noisy responses

X = np.linspace(0,14,15)

# Define a vector of noise variance at individual data locations

noise_var = np.array([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 set with the noise

Y = X

np.sin(X) + np.sqrt(noise_var)

np.random.randn(X.size)

# Define a Kriging metamodel

MetaOpts = {

'Type': 'Metamodel',

'MetaType': 'Kriging',

'ExpDesign': {

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

# Impose the known heteroscedastic noise variance

'Regression': {

'SigmaNSQ': noise_var.tolist()

}

# Create the Kriging metamodel

myKriging = 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

example for the different cases of noise variance speciﬁcation is provided in

uq_Example_Kriging_08_Regression.

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 speciﬁed, the experimental design) is having vector-valued (multi-

ple) outputs, UQ[PY]LAB performs an independent Kriging metamodel for each output com-

ponent with shared (common) experimental design and trend basis functions. No additional

conﬁguration is needed to enable this behavior. A Kriging example for model with multiple

outputs can be found in the UQ[PY]LAB example script uq_Example_Kriging_07_MultipleOutputs.

UQ[PY]LAB-V1.0-105 - 36 -

Kriging (Gaussian process modeling)

0 5 10

−15

−10

−5

Kriging approximation

95% confidence interval

Observations

𝑋

𝑌

(a) known homoscedastic noise

0 5 10

−15

−10

−5

Kriging approximation

95% confidence interval

Observations

𝑋

𝑌

(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.1 Accessing the results

Calculating a Kriging metamodel on a model with multiple outputs will result in an output

list of dictionaries. As an example, a model with nine outputs will produce a list with nine

dictionaries for myKriging['Kriging'], each dictionary having the following keys:

beta

sigmaSQ

theta

Each element of the Kriging dictionary is functionally identical to its scalar counterpart in

Section 2.3.1. Similarly, the myKriging['Error'] dictionary becomes a list of dictionaries

with key:

LOO

2.7 Using a validation set

If an independent validation set is provided (see Table 19 in Section 3.1.8), UQ[PY]LAB

automatically computes the validation error according to Eq. (1.62). To provide a validation

set stored in, for example, X_val and Y_val, the following option is set:

MetaOpts['ValidationSet'] = {

'X': X_val,

'Y': Y_val

}

The value of the validation error is stored in myKriging['Error']['Val'] (see Table 22)

UQ[PY]LAB-V1.0-105 - 37 -

UQ[PY]LAB user manual - Python

and will also be displayed when invoking uq.print(myKriging).

2.8 Using Kriging predictor as a model

Regardless of the conﬁguration 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, af-

ter a Kriging MODEL object is created in UQ[PY]LAB, it can be used just like an ordinary

UQ[PY]LAB MODEL object (for details, see the UQ[PY]LAB User Manual – the MODEL mod-

ule).

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:

x = 1

YmuKG = uq.evalModel(myKriging, x)

The variance of the predictor can be also evaluated by using a two-component output:

YmuKG,YvarKG = uq.evalModel(myKriging, x)

As most functions within UQ[PY]LAB, model evaluations are vectorized. Therefore, evalu-

ating multiple points at a time is much faster than repeatedly evaluating one point at a time.

For example:

# We use reshape(-1,1) to ensure that X is a column vector

X = np.arange(0, 15, 0.1).reshape(-1,1)

# Evaluate the mean and variance of the Kriging predictor on X

YmuKG,YvarKG = uq.evalModel(myKriging, X)

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

The Kriging module in UQ[PY]LAB 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 UQ[PY]LAB 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.

X = np.arange(0, 15, 2)

# Calculate the response of the model

Y = X

np.sin(X)

# Define a custom Kriging object

MetaOpts = {

'Type': 'Metamodel',

'MetaType': 'Kriging',

'ExpDesign': {

UQ[PY]LAB-V1.0-105 - 38 -

Kriging (Gaussian process modeling)

'Sampling': 'User',

'X': X.tolist(),

'Y': Y.tolist()

'Kriging': {

'Trend': { #Select ordinary Kriging

'Type': 'ordinary'

# Set the values of beta, sigmaˆ2 and theta

'beta': 69.84,

'sigmaSQ': 2.566e5,

'theta': 9.999,

# Select ellipsoidal, Matern-3/2 correlation function

'Corr': {

'Type': 'ellipsoidal',

'Family': 'matern-3_2'

}

# Create the metamodel

myCustomKriging = uq.createModel(MetaOpts)

# Evaluate the metamodel on some new points

X_new = np.arange(1, 13, 2).reshape(-1,1)

Y_new = uq.evalModel(myCustomKriging, X_new)

If the metamodel has more than one output, it is sufﬁcient to specify the same informa-

tion for each output in MetaOpts['Kriging'], which now needs to be a list of dictio-

naries. Furthermore, to deﬁne a custom regression model, a value to an additional ﬁeld

MetaOpts['Kriging']['sigmaNSQ'] must be speciﬁed.

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, UQ[PY]LAB may use the scaled experimental

design U. Depending on the value of the option MetaOpts['Scaling'] and whether a prob-

abilistic input model has been deﬁned (in MetaOpts['Input']), the transformation X 7→ U

varies. All the possible cases are summarized in Table 3.

Table 3: Available experimental design transformations

Scaling

Input

No INPUT deﬁned INPUT deﬁned

true U = (X − µ

) /σ

U = (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

UQ[PY]LAB-V1.0-105 - 39 -

UQ[PY]LAB user manual - Python

zero mean and a unit variance (i.e., standardized). If no input model has been speciﬁed, then

the mean and variance (µ

, σ

) are empirically estimated from the available data speciﬁed

in MetaOpts['ExpDesign']['X']. If, however, an input model has been speciﬁed via an

INPUT object, say myInput, as follows:

MetaOpts['Input'] = myInput

then (µ

, σ

) are estimated from the moments of the marginal distribution of each compo-

nent in X.

If an additional INPUT object, say my_scaled_input, is set to the MetaOpts['Scaling']

option:

MetaOpts['Scaling'] = my_scaled_input

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

two probability spaces (for more information, refer to the UQ[PY]LAB 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

As mentioned in Section 1.2, a Kriging predictor is a Gaussian process posterior (with respect

to the experimental design). Thus, instead of retrieving the mean and conﬁdence bounds

of the Gaussian process posterior, users may also draw an arbitrary number of sample paths

or realizations from the process. This can be achieved by exploiting the Gaussian property

of any Kriging metamodel. Each sample path of the Gaussian process Y , discretized over

a set of points X



(1)

, . . . , x

)



, follows a multivariate Gaussian distribution, i.e.,

Y ∼ N



, C



where µ

and C

corresponds to the posterior mean and the posterior

covariance matrix of the Kriging predictor evaluated on X

, respectively.

Given the myKriging object that was calculated in Section 2.3, the posterior mean and

covariance matrix of the Kriging predictor over a set of points are obtained as follows:

Ns = 500

Xs = np.linspace(0, 15, Ns).reshape(-1, 1)

mu, _, C = uq.evalModel(myKriging, Xs, nargout=3)

Then, one sample path of the posterior Gaussian process can be generated as follows:

Ys = mu + np.matmul(np.linalg.cholesky(C), \

np.random.multivariate_normal(np.zeros(Ns), np.eye(Ns), 1).T)

A set of 15 such sample paths is plotted in Figure 14 together with the mean of the Kriging

predictor.

UQ[PY]LAB-V1.0-105 - 40 -

Kriging (Gaussian process modeling)

0 5 10 15

−10

Kriging approximation

Observations

Sample Paths

𝑋

𝑌

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).

UQ[PY]LAB will automatically and appropriately account for the set of input parameters

which are declared constant.

To set a parameter to constant, e.g., with const_value as its value, the following options are

speciﬁed (for details, see UQ[PY]LAB User Manual – the INPUT module):

InputOpts = {

'Marginals': [{

'Type': 'Constant',

'Parameters': [const_value]

}]

}

Furthermore, when the standard deviation of an input parameter is set to zero, UQ[PY]LAB

automatically sets the marginal of this parameter to the type Constant. For example, the fol-

lowing 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': 'Constant',

'Parameters': [1, 1]

}]

}

UQ[PY]LAB-V1.0-105 - 41 -

Chapter 3

Reference List

How to read the reference list

Python dictionaries play an important role throughout the UQ[PY]LAB syntax. They offer

a natural way to semantically group conﬁguration options and output quantities. Due to

the complexity of the algorithms implemented, it is not uncommon to employ nested dictio-

naries to ﬁne-tune the inputs and outputs. Throughout this reference guide, a table-based

description of the conﬁguration dictionaries is adopted.

The simplest case is given when a value of a dictionary key is a simple value or a list:

Table X: Input

Name String A description of the ﬁeld 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 key-value pair. At the beginning of each row a symbol is given to inform as to whether

the corresponding key 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 keys 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)

Kriging (Gaussian process modeling)

When the value of one of the keys of a dictionary is a dictionary itself, a link to a table that

describes the structure of that nested dictionary is provided, as in the case of the Options

key in the following example:

Table X: Input

Name String Description

□ Options Table Y Description of the Options

dictionary

Table Y: Input['Options']

Key1 String Description of Key1

□ Key2 Double Description of Key2

In some cases, an option value gives the possibility to deﬁne further options related to that

value. The general syntax would be:

Input = {

'Option1' : 'VALUE1',

'VALUE1' : {

'Val1Opt1' : ... ,

'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 Float Description

Table Z: Input['VALUE2']

□ Val2Opt1 String Description

UQ[PY]LAB-V1.0-105 - 43 -

UQ[PY]LAB user manual - Python

□ Val2Opt2 Float Description

UQ[PY]LAB-V1.0-105 - 44 -

Kriging (Gaussian process modeling)

3.1 Create a Kriging metamodel

Syntax

myKriging = uq.createModel(MetaOpts)

Input

The dictionary variable MetaOpts contains the conﬁguration information for a Kriging

metamodel. 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 identiﬁer for the metamodel

□ ExpDesign Table 5 Option to specify an experimental

design (Section 2.4.1)

□ Input String Option to specify an INPUT object

that describes the inputs of the

metamodel. For example, you specify

the INPUT myInput by supplying

myInput['Name']

□ FullModel String Option to specify the MODEL object

that describes the full computational

model of the metamodel. The object

is used to compute the model

responses. For example, you specify

the MODEL myModel by supplying

myModel['Name']

□ Trend Table 6 Options to specify the Kriging trend

(Sections 1.3 and 2.4.2)

□ Corr Table 8 Options to specify the correlation

function (Section 2.4.3)

□ EstimMethod String

default: 'CV'

Select the method to estimate

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.

UQ[PY]LAB-V1.0-105 - 45 -

UQ[PY]LAB user manual - Python

□ 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 Boolean

default: True

Flag to keep (cached) some

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 Boolean or String

default: true

Conduct Kriging calculation in an

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).

Boolean Flag to scale the experimental

design. If set to True , the

experimental design will be scaled

(speciﬁcally, standardized) before

calculating the Kriging metamodel.

String The name of the INPUT object that

deﬁnes the auxiliary probabilistic

space

□ ValidationSet Table 19 Independent validation data set

(Section 2.7)

3.1.1 Experimental design options

Table 5: MetaOpts['ExpDesign']

Sampling String

default: 'LHS'

Type of sampling

'MC' Monte Carlo sampling

'LHS' Latin Hypercube sampling

'Sobol' Sobol’ sequence sampling

'Halton' Halton sequence sampling

UQ[PY]LAB-V1.0-105 - 46 -

Kriging (Gaussian process modeling)

'User' User-deﬁned experimental design,

deﬁned in

MetaOpts['ExpDesign']['X']

and

MetaOpts['ExpDesign']['Y'].

□ NSamples Integer Number of sample points to draw. It

is required for all types of sampling.

□ X N × M Numpy Array of

Floats

User-deﬁned experimental design

(model inputs) X. It only applies,

and is required, if the type of

sampling is 'User'.

□ Y N × M Numpy Array of

Floats

User-deﬁned 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

default: 'ordinary'

Type of the Kriging trend

(Sections 1.3 and 2.4.2)

'simple' Known trend function (i.e., simple

Kriging), deﬁned in

.Trend.CustomF. deﬁned in

metaopts['Trend']['CustomF'].

No estimation of the coefﬁcients

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, deﬁned in

MetaOpts['Trend']['Degree'].

'custom' User-deﬁned trend (Section 2.4.2.2)

□ Degree Integer

default: 0

Degree of the polynomial trend. It

only applies if the Kriging trend type

is 'polynomial'.

□ CustomF Float , String, or List Field to set up custom trend function

Float Constant trend in simple Kriging

String Function name that returns f(x), an

arbitrary trend function

UQ[PY]LAB-V1.0-105 - 47 -

UQ[PY]LAB user manual - Python

List of strings Respective names of functions

(x) , i = 1, . . . , P , such that each

function evaluated on Xreturns the

i’th column of F

□ TruncOptions Table 7 Polynomial basis truncation options

Table 7: MetaOpts.Trend['TruncOptions']

□ qNorm Float

default: 1

Parameter for the q-norm (or

hyperbolic) truncation scheme. For

details, see the UQ[PY]LAB 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 UQ[PY]LAB User Manual –

Polynomial Chaos Expansions.

□ Custom List of Floats User-deﬁned basis speciﬁed using

multi-indices

3.1.3 Correlation function options

Table 8: MetaOpts['Corr']

□ Family String

default: 'Matern-5_2'

One-dimensional correlation family

(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

ern-3/2 correlation (Eq. (1.32))

'Matern-5_2' Mat

ern-5/2 correlation (Eq. (1.33))

□ Type String

default:

'Ellipsoidal'

Type of the correlation function

'Ellipsoidal' Ellipsoidal correlation (Eq. (1.34))

'Separable' Separable correlation (Eq. (1.35))

□ Isotropic Boolean

default: false

Flag to determine whether the

correlation function is isotropic or

anisotropic (Section 1.4.3)

UQ[PY]LAB-V1.0-105 - 48 -

Kriging (Gaussian process modeling)

□ Nugget Float scalar or List

default: 10

−10

Set of values that are added to the

main diagonal of the correlation

matrix R (Section 1.4.4)

• If scalar, add this quantity to

the diagonal elements of R.

• If List , add each element to

the corresponding diagonal

element of R.

3.1.4 Estimation method options

Table 9: MetaOpts['CV']

□ LeaveKOut Integer (< N )

default: 1

Number of sample points left out in

cross-validation (CV) (Section 2.4.4)

• The number of sample points

in D

(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 Float

default: 1.0

Initial estimate of the correlation

parameters

Scalar Float Initial estimate of the correlation

parameter. If M > 1, the same value

is assigned for all input dimensions.

M × 1 Float Initial estimates of the correlation

parameters for each of the M input

dimensions

□ Bounds 2 × M or 2 × 1 Float

default: [10

−3

, 10]

Bound of admissible values of the

correlation parameters, in the form

[lower_bound; upper_bound].

UQ[PY]LAB-V1.0-105 - 49 -

UQ[PY]LAB user manual - Python

2 ×1 Float Bound of the correlation parameter.

If M > 1, the same value is assigned

for all input dimensions.

2 ×M Float Bounds of the correlation parameters

for each of the M input dimensions

□ Display String

default: 'none'

Option to determine the verbosity of

the optimization process

'none' Nothing will be printed to the

command window

'final' Only the ﬁnal result of the

optimization process will be printed

'iter' State of the optimization process will

be printed after each iteration

□ MaxIter Integer

default: 20

Maximum number of iterations or

generations

□ Tol Float

default: 10

−4

Covergence tolerance of the

optimization

□ Method String

default: 'HGA'

Optimization method (Section 2.4.5)

'none' No optimization. The value

MetaOpts['Optim']['InitialValue']

is used.

'LBFGS' Gradient-based optimization

(L-BFGS) using the fmincon

function of MATLAB

'GA' Genetic Algorithm (GA) optimization

using the ga function of MATLAB

'HGA' Hybrid Genetic Algorithm (HGA)

optimization using the functions ga

and fmincon of MATLAB

'CMAES' Covariance Matrix

Adaptation-Evolution Strategy

(CMA-ES) optimization using the

function uq_cmaes of UQLIB

'HCMAES' Hybrid Covariance Matrix

Adaptation-Evolution Strategy

(HCMA-ES) optimization using the

functions uq_cmaes of UQLIB and

fmincon of MATLAB

⊞ 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

UQ[PY]LAB-V1.0-105 - 50 -

Kriging (Gaussian process modeling)

⊞ CMAES Table 14 Options relevant to the CMA-ES

optimization method

⊞ HCMAES Table 15 Options relevant to the HCMA-ES

optimization method

Note: The convergence tolerance deﬁned under Metaopts.Optim.Tol may slightly dif-

fer depending on the optimization method. For MATLAB optimization functions

(fmincon and ga), check the deﬁnition of the option 'TolFun' of the respective

function.

Table 11: MetaOpts['Optim']['BFGS']

□ nLM Integer (≥ 1)

default: 5

Limited memory size of the L-BFGS

method ( nLM ≥ 1)

Table 12: MetaOpts['Optim']['GA']

□ nPop Integer

default: 30

Population size of each generation

□ nStall Integer

default: 5

Maximum number of stall

generations

Table 13: MetaOpts['Optim']['HGA']

□ nPop Integer

default: 30

Population size of each generation

□ nStall Integer

default: 5

Maximum number of stall

generations

□ nLM Integer (≥ 1)

default: 5

Limited memory size of the L-BFGS

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

default: 30

The population size of each

generation

□ nStall Integer

default: 5

The maximum number of stall

generations

Table 15: MetaOpts['Optim']['HCMAES']

□ nPop Integer

default: 30

Population size of each generation

□ nStall Integer

default: 5

Maximum number of stall

generations

UQ[PY]LAB-V1.0-105 - 51 -

UQ[PY]LAB user manual - Python

□ nLM Integer

default: 5

Limited memory size of the L-BFGS

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',

Boolean, or Float

default: 'none'

Speciﬁcation of the output noise

variance in a regression model

(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 Float Homoscedastic noise variance σ

N × 1 Float Independent heterogeneous

(heteroscedastic) noise variance σ

N × N Float Noise covariance matrix Σ

□ SigmaSQ Table 17 Initial value and bound for the

optimization of the Gaussian process

variance σ

in the case of a

regression model with known noise

variance

Note: If the metamodel has more than one outputs, then MetaOpts['Regression']

is expected to be a list with length equal to the number of outputs and that the

parameters in Table 16 are deﬁned for each MetaOpts['Regression'][idx]

(where idx ranges from 0 to the total number of outputs−1). If only one regres-

sion option is speciﬁed, then the option applies to all output by default.

Table 17: MetaOpts['Optim']['Regression']['SigmaSQ']

□ InitialValue Float

default: 0.5Var [Y]

Initial value of the Gaussian process

variance

□ Bound 2 × 1 Float

default:

[0.1 ×Var [Y] ; 10 × Var [Y]]

Lower and upper bounds for the

Gaussian process variance

where Var [Y] =

− ¯y)

and ¯y =

are the sample variance

and mean of the observed responses, respectively.

UQ[PY]LAB-V1.0-105 - 52 -

Kriging (Gaussian process modeling)

3.1.7 Custom Kriging options

To create a user-deﬁned Kriging model as desribed in Section 2.9, a MetaOpts['Kriging']

options must be set according to Table 18.

Table 18: MetaOpts['Kriging']

Trend Table 6 Kriging trend deﬁnition

beta Float Values of the trend coefﬁcients β in

Eq. (2.2)

sigmaSQ Float Value of σ

in Eq. (2.2)

theta Float Values of θ in Eq. (2.2)

Corr Table 8 Correlation function deﬁnition (See

below Note)

□ SigmaNSQ Scalar, List , or List of

lists with Float entries

Values of the output noise variance

in a regression model (Sections 1.2.2

and 1.2.2)

Scalar Float Homoscedastic noise variance σ

N × 1 Float Independent heteroscedastic noise

variance σ

N × N Float Noise covariance matrix Σ

Note: In the current version of UQ[PY]LAB, only the default correlation function han-

dle uq_eval_Kernel is supported. All other options follow Table 8.

Note: If the metamodel has more than one outputs, then MetaOpts['Kriging'] is

expected to be a list with length equal to the number of outputs and that the pa-

rameters in Table 18 are deﬁned for each MetaOpts['Kriging'][idx] (where

idx ranges from 0 to the total number of outputs−1).

Note: To create a custom Kriging model, the ﬁelds 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, UQ[PY]LAB automatically calculates the vali-

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

Table 19.

Table 19: MetaOpts['ValidationSet']

UQ[PY]LAB-V1.0-105 - 53 -

UQ[PY]LAB user manual - Python

X N × M Float User-speciﬁed validation set

experimental design X

val

Y N × N

out

Float User-speciﬁed validation set response

val

UQ[PY]LAB-V1.0-105 - 54 -

Kriging (Gaussian process modeling)

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 deﬁned in MetaOpts

variable (Section 3.1)

Internal Table 24 Internal ﬁelds

Table 21: myKriging['Kriging']

beta P × 1 Float Values of β(θ) in Eq. (2.2)

sigmaSQ Float Value of σ

(θ) in Eq. (2.2)

theta Scalar or

1 ×M Float

Values of θ in Eq. (2.2)

sigmaNSQ Scalar, List , or

List of lists with

Float entries

Value of noise variance in Gaussian process

regression (Sections 1.2.2 and 1.2.2)

Scalar Float Homoscedastic noise variance σ

N × 1 Float Independent heteroscedastic noise variance σ

N × N Float Noise covariance matrix Σ

Table 22: myKriging['Error']

LOO Float Leave-One-Out (LOO) error, calculated using

Eq. (1.61)

Val Float 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 ﬁelds myKriging['Kriging'] and myKriging['Error'] are

lists 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 = 0, ..., N

out

− 1.

Table 23: myKriging['ExpDesign']

NSamples Float Number of sample points

UQ[PY]LAB-V1.0-105 - 55 -

UQ[PY]LAB user manual - Python

Sampling String Sampling method

X N × M Float Values in the experimental design

U N × M Float Experimental design values in the auxiliary

(scaled) space. For details, see Section 2.10

Y N × N

out

Float Observed model responses (output) Y that

corresponds to the experimental design X

3.2.1 Internal ﬁelds (advanced)

Note: The internal ﬁelds 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 ﬁelds are not documented in this user manual.

Table 24: myKriging['Internal']

Kriging Table 25 Internal ﬁelds with data related to the Kriging

model

Runtime Dictionary Internal ﬁelds with variables that are used

during the calculation of the kriging

metamodel

Error Table 26 Internal ﬁelds related to the metamodeling

error

ExpDesign Table 27 Internal ﬁelds related to the experimental

design

Scaling Boolean or

INPUT object

Scaling of the experimental design (see

Section 2.10)

Regression Table 28 Internal ﬁelds 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 ﬁelds related to the trend

GP Table 30 Internal ﬁelds related to the Gaussian process

Optim Table 31 Internal ﬁelds related to the optimization of

the hyperparameters

sigmaNSQ Float Internal ﬁeld that contains the value of noise

variance in Gaussian process regression (refer

to the entry sigmaNSQ in Table 21)

Cached Dictionary Internal ﬁelds related to cached variables

UQ[PY]LAB-V1.0-105 - 56 -

Kriging (Gaussian process modeling)

Note: The ﬁelds myKriging['Internal']['Kriging'] and

myKriging['Internal']['Error'] are, in general, list objects 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 = 0, ..., N

out

− 1.

Table 26: myKriging.Internal['Error']

LOOmean N × 1 Float 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 Float Standard deviation of the Kriging predictor

from LOO CV

varY 1 × N

out

Float 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 Float Mean of the inputs used in the scaling of the

experimental design (Section 2.10).

• If an experimental design is speciﬁed in

MetaOpts['ExpDesign'], then this

corresponds to the empirical mean of

each input dimension.

• If INPUT object is speciﬁed in

MetaOpts['Input'], then this

corresponds to the ﬁrst moment of the

corresponding marginal of each input

dimension.

stdX 1 × M Float Standard deviation of the inputs used in the

scaling of the experimental design

(Section 2.10).

• If an experimental design is speciﬁed in

MetaOpts['ExpDesign'], then this

corresponds to the empirical standard

deviation of each input dimension.

• If INPUT object is speciﬁed in

MetaOpts['Input'], then this

corresponds to the second moment of

the corresponding marginal of each

input dimension.

UQ[PY]LAB-V1.0-105 - 57 -

UQ[PY]LAB user manual - Python

varY 1 × N

out

Float 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 Boolean Flag that indicates the current Kriging model is

a Gaussian process regression model

EstimNoise Boolean Flag that indicates that the noise variance is

estimated

Tau Dictionary Dictionary that contains the initial value and

bounds used in the optimization of the τ

parameter (Eq. (1.21))

SigmaSQ Dictionary Dictionary that contains the initial value and

bounds used in the optimization of the σ

(Gaussian process variance)

SigmaNSQ Scalar, List , or

List of lists with

Float entries

Estimated or speciﬁed noise variance.

IsHomoscedastic Boolean Flag that indicates the Gaussian process

regression model is homoscedastic, with a

constant noise variance.

Table 29: myKriging.Internal.Kriging['Trend']

F N × P Float Observation matrix, i.e., the basis functions of

the Kriging trend f

’s evaluated on the

experimental design (Eq. (2.2)).

beta p × 1 Float Values of β(

θ) in Eq. (2.2)

Table 30: myKriging.Internal.Kriging['GP']

R N × N Float Correlation matrix of the Gaussian process on

the experimental design points (i.e.,

R = R(X, X))

sigmaSQ Float Value of σ

in Eq. (2.2)

Table 31: myKriging.Internal.Kriging['Optim']

Theta Scalar or

1 ×M Float

Optimum value of θ

Tau Float Ratio of σ

to σ

+ σ

as deﬁned in Eq. (1.21).

Only available in a regression with unknown

homoscedastic noise variance.

UQ[PY]LAB-V1.0-105 - 58 -

Kriging (Gaussian process modeling)

SigmaSQ Float 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 Float Objective function value at the optimum θ

InitialObjFun Float Objective function value at the initial estimate

of θ

nEval Float Number of objective function evaluations

during the optimization process

nIter Float 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.

UQ[PY]LAB-V1.0-105 - 59 -

UQ[PY]LAB user manual - Python

3.3 Kriging predictor

Syntax

Y_mu = uq.evalModel([Model=None, ]X=None[, nargout=1])

Y_mu, Y_sigma2 = uq.evalModel(..., nargout=2)

Y_mu, Y_sigma2, Y_cov = uq.evalModel(..., nargout=3)

Description

Y_mu = uq.evalModel(X) returns the mean of the Kriging predictor (N × N

out

) 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 ×

out

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

[Y_mu,Y_sigma2] = uq.evalModel(..., nargout=2) additionally returns the variance

of the Kriging predictor (N × N

out

[Y_mu,Y_sigma2,Y_cov] = uq.evalModel(..., nargout=3) additionally returns the co-

variance matrices of the Kriging predictor (N × N × N

out

), i.e., one covariance matrix

per output dimension. Note that the diagonal elements of Y_cov (in each output dimen-

sion) are equal to the corresponding elements in Y_sigma2, i.e., Y_sigma2 = diag(Y_cov).

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

model.

UQ[PY]LAB-V1.0-105 - 60 -

Kriging (Gaussian process modeling)

3.4 Printing and visualizing a Kriging metamodel

UQ[PY]LAB offers two commands to conveniently print reports containing contextually rele-

vant information for a given Kriging metamodel object.

3.4.1 Printing the results: uq.print

Syntax

uq.print(myKriging[, outIdx[, 'Option1'[, 'Option2', ...]]])

Description

uq.print( myKriging ) prints a report of the conﬁguration options and results of the

Kriging metamodel object myKriging . If the model has multiple outputs, only the

report on the information about the ﬁrst output dimension are printed.

uq.print( myKriging , outIdx) prints a report on the conﬁguration options and results

of the Kriging metamodel object myKriging for the output dimensions speciﬁed in

the array outIdx.

uq.print( myKriging , outIdx, 'Option1', 'Option2', ...) prints a report only

on selected conﬁguration options or results, the selection of which, are speciﬁed by the

string options 'Option1', 'Option2', etc.See Table 32 for the available options.

Table 32: uq.print options

'beta' Prints out the regression coefﬁcients β

'theta' Prints out the ﬁnal hyperparameters θ (either the optimal value

or the user-speciﬁed value if they are manually speciﬁed)

'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, coefﬁcients)

'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

Examples

uq.print( myKriging ,[1 3]) prints the Kriging metamodeling results for output dimen-

sions 1 and 3.

UQ[PY]LAB-V1.0-105 - 61 -

UQ[PY]LAB user manual - Python

uq.print( myKriging , 3, 'theta', 'R') only prints the hyperparameters θ and cor-

relation matrix R for output dimension 3.

UQ[PY]LAB-V1.0-105 - 62 -

Kriging (Gaussian process modeling)

3.4.2 Visualize the results: uq.display

Syntax

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 ﬁrst output

dimension is plotted.

uq.display( myKriging , outIdx) plots the Kriging model predictions against the input

for the model output dimensions speciﬁed 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 speciﬁed 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

dimension, 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 ﬁgures.

UQ[PY]LAB-V1.0-105 - 63 -

References

Abramovitz, M. and I. A. Stegun (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 misspeciﬁcation. Computational Statistics

and Data Analysis 66, 55–69. 15

Byrd, R. H., M. E. Hribar, and J. Nocedal (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. Ph. D. thesis, Universit

e Blaise Pascal, Clermont-Ferrand, France. 3, 6,

11, 13

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

Addison-Wesley Professional. 16

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

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

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 C. Williams (2006). Gaussian processes for machine learning. Adaptive

computation and machine learning. Cambridge, Massachusetts: MIT Press. 1, 5, 10, 11,

Sacks, J., W. J. Welch, T. J. Mitchell, and H. P. Wynn (1989). Design and analysis of computer

experiments. Statistical Science 4, 409–435. 1, 11

Santner, T., B. Williams, and W. Notz (2003). The design and analysis of computer experiments.

Springer series in Statistics. Springer. 1, 2, 3, 13, 15