UQ[PY]LAB USER MANUAL
POLYNOMIAL CHAOS EXPANSIONS
S. Marelli, N. L
¨
uthen, 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 Quantification in the cloud with UQCloud, Proceedings of the 4th
International Conference on Uncertainty Quantification in Computational Sciences and Engineering (UNCECOMP
2021), Athens, Greece, June 27–30, 2021.
How to cite this manual
S. Marelli, N. L
¨
uthen, B. Sudret, UQ[py]Lab user manual – Polynomial chaos expansions, Report UQ[py]Lab
-V1.0-104, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2024
BIBT
E
X entry
@TechReport{UQdoc_10_104,
author = {Marelli, S. and L\"uthen, N. and Sudret, B.},
title = {{UQ[py]Lab user manual -- Polynomial chaos expansions }},
institution = {Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich,
Switzerland},
year = {2024},
note = {Report UQ[py]Lab - V1.0-104}
}
List of contributors:
Name Contribution
M. Berchier Implementation and documentation of the orthogonal matching pursuit (OMP) regression method
P. Wiederkehr Implementation and documentation of the q-norm-adaptivity for regression methods
A. Hlobilov
´
a Translation from the UQLab manual
Document Data Sheet
Document Ref. UQ[PY]LAB-V1.0-104
Title: UQ[PY]LAB user manual – Polynomial chaos expansions
Authors: S. Marelli, N. L
¨
uthen, B. Sudret
Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich,
Switzerland
Date: 27/05/2024
Doc. Version Date Comments
V0.8 30/09/2022 Initial release
V1.0 27/05/2024 Minor update for UQ[PY]LAB release 1.0
Abstract
Polynomial Chaos Expansions (PCE) are a powerful metamodelling tool that has important
applications in many engineering and applied mathematics fields. These include structural
reliability, sensitivity analysis, Monte Carlo simulation and others. Due to the underlying
complexity of its formulation, however, this technique has seen relatively little use outside of
these fields.
UQ[PY]LAB metamodelling tools provide an efficient, flexible and easy to use PCE module
that allows one to apply state-of-the-art algorithms for non-intrusive, sparse and adaptive
PCE on a variety of applications. This manual for the polynomial chaos expansion metamod-
elling module is divided into three parts:
• A short introduction to the main concepts and techniques behind PCE, with a selection
of references to the relevant literature;
• A detailed example-based guide, with the explanation of most of the available options
and methods;
• A comprehensive reference list detailing all the available functionalities in UQ[PY]LAB.
Keywords: UQCloud, UQ[PY]LAB, metamodelling, Polynomial Chaos Expansions, PCE,
Sparse PCE
Contents
1 Theory 1
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Polynomial chaos expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Building the polynomial basis . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3.1 Classical families of univariate orthonormal polynomials . . . . . . . . 2
1.3.2 Extension to arbitrary distributions . . . . . . . . . . . . . . . . . . . . 3
1.3.3 Basis truncation schemes . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.4 Basis-adaptive PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 A posteriori error estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.1 Normalized empirical error . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2 Leave-one-out cross-validation error . . . . . . . . . . . . . . . . . . . 7
1.4.3 Corrected error estimates . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5 Calculation of the coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.5.1 Projection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.5.2 Least-Squares regression . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5.3 Sparse PCE: Least Angle Regression . . . . . . . . . . . . . . . . . . . . 11
1.5.4 Sparse PCE: Orthogonal Matching Pursuit . . . . . . . . . . . . . . . . 13
1.5.5 Sparse PCE: Subspace Pursuit . . . . . . . . . . . . . . . . . . . . . . . 15
1.5.6 Sparse PCE: Bayesian compressive sensing . . . . . . . . . . . . . . . . 16
1.6 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.6.1 Moments of a PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.6.2 Uncertainty propagation and sensitivity analysis . . . . . . . . . . . . . 19
1.7 Bootstrap-based error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.7.1 Bootstrap PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.7.2 Fast bPCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2 Usage 21
2.1 Reference problem: the Ishigami function . . . . . . . . . . . . . . . . . . . . 21
2.2 Problem set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2.1 Full model and probabilistic input model . . . . . . . . . . . . . . . . . 22
2.3 Set-up of the polynomial chaos expansion . . . . . . . . . . . . . . . . . . . . 22
2.4 Orthogonal polynomial basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.1 Univariate polynomial types . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.2 Truncation schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Experimental design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.6 Calculation of the coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.6.1 Projection: Gaussian quadrature . . . . . . . . . . . . . . . . . . . . . 25
2.6.2 Ordinary Least-Squares (OLS) . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.3 Sparse regression methods (LARS, OMP, SP, BCS) . . . . . . . . . . . . 31
2.7 Basis adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.7.1 Degree-Adaptive PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.7.2 q-norm-Adaptive PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.8 Use of a validation set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.9 Manually specify inputs and computational models . . . . . . . . . . . . . . . 40
2.10 PCE of vector-valued models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.10.1 Accessing the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
2.11 Using a PCE as a predictor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.11.1 Evaluate a PCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.11.2 Local error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.12 Manually specifying PCE parameters (predictor-only mode) . . . . . . . . . . . 43
2.13 Using a PCE with constant input variables . . . . . . . . . . . . . . . . . . . . 44
3 Reference List 45
3.1 Create a PCE metamodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.1.1 Truncation options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.1.2 PCE Coefficients calculation options . . . . . . . . . . . . . . . . . . . 49
3.1.3 Quadrature-specific options . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.4 OLS-specific options . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.5 LARS-specific options . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.6 OMP-specific options . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.1.7 SP-specific options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.1.8 BCS-specific options . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.9 Experimental design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.10 Validation Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.11 Bootstrap options . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2 Accessing the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2.1 Polynomial chaos expansion information . . . . . . . . . . . . . . . . . 53
3.2.2 Experimental design information . . . . . . . . . . . . . . . . . . . . . 54
3.2.3 Error estimates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.2.4 Internal keys (advanced) . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Chapter 1
Theory
1.1 Introduction
In most modern engineering contexts (and in applied sciences in general), uncertainty quan-
tification is becoming an increasingly important field. Deterministic scenario-based predictive
modelling is being gradually substituted by stochastic modelling to account for the inevitable
uncertainty in physical phenomena and measurements. This smooth transition, however,
comes at the cost of dealing with greatly increased amounts of information (e.g. when using
Monte-Carlo simulation), usually resulting in the need to perform expensive computational
model evaluations repeatedly.
Metamodelling (or surrogate modelling) attempts to offset the increased costs of stochastic
modelling by substituting the expensive-to-evaluate computational models (e.g. finite ele-
ment models, FEM) with inexpensive-to-evaluate surrogates. Polynomial chaos expansions
(PCE) are a powerful metamodelling technique that aims at providing a functional approxi-
mation of a computational model through its spectral representation on a suitably built basis
of polynomial functions.
Due to the large scope of the UQ[PY]LAB software framework, its polynomial chaos ex-
pansions module offers extensive facilities for the deployment of a number of non-intrusive
PCE calculation techniques. This part is intended as an overview of the relevant theory and
literature in this field.
1.2 Polynomial chaos expansion
Consider a random vector with independent components X ∈ R
M
described by the joint
probability density function (PDF) f
X
. Consider also a finite variance computational model
as a map Y = M(X), with Y ∈ R such that:
E
Y
2
=
Z
D
X
M
2
(x)f
X
(x)dx < ∞. (1.1)
Then the polynomial chaos expansion of M(X) is defined as:
1
UQ[PY]LAB user manual
Y = M(X) =
X
α∈N
M
y
α
Ψ
α
(X) (1.2)
where the Ψ
α
(X) are multivariate polynomials orthonormal with respect to f
X
, α ∈ N
M
is a multi-index that identifies the components of the multivariate polynomials Ψ
α
and the
y
α
∈ R are the corresponding coefficients (coordinates).
In realistic applications, the sum in Eq. (1.2) needs to be truncated to a finite sum, by intro-
ducing the truncated polynomial chaos expansion:
M(X) ≈ M
P C
(X) =
X
α∈A
y
α
Ψ
α
(X) (1.3)
where A ⊂ N
M
is the set of selected multi-indices of multivariate polynomials. The construc-
tion of such a set is detailed in Section 1.3.3.
1.3 Building the polynomial basis
The polynomial basis Ψ
α
(X) in Eq. (1.3) is traditionally built starting from a set of univariate
orthonormal polynomials ϕ
(i)
k
(x
i
) which satisfy:
D
ϕ
(i)
j
(x
i
), ϕ
(i)
k
(x
i
)
E
def
=
Z
D
X
i
ϕ
(i)
j
(x
i
)ϕ
(i)
k
(x
i
)f
X
i
(x
i
)dx
i
= δ
jk
(1.4)
where i identifies the input variable w.r.t. which they are orthogonal as well as the corre-
sponding polynomial family, j and k the corresponding polynomial degree, f
X
i
(x
i
) is the
i
th
-input marginal distribution and δ
jk
is the Kronecker symbol. Note that this definition of
inner product can be interpreted as the expectation value of the product of the factors.
The multivariate polynomials Ψ
α
(X) are then assembled as the tensor product of their uni-
variate counterparts:
Ψ
α
(x)
def
=
M
Y
i=1
ϕ
(i)
α
i
(x
i
) (1.5)
Due to the orthonormality relations in Eq. (1.4), it follows that also the multivariate polyno-
mials thus constructed are orthonormal:
⟨Ψ
α
(x), Ψ
β
(x)⟩ = δ
αβ
(1.6)
where δ
αβ
is an extension Kronecker symbol to the multi-dimensional case.
1.3.1 Classical families of univariate orthonormal polynomials
The classical families of univariate orthonormal polynomials and the distributions to which
they are orthonormal are given for reference in Table 1 (?). Detailed descriptions of each
of the classical families of polynomials (also called Askey-Scheme orthonormal polynomials)
UQ[PY]LAB-V1.0-104 - 2 -
Polynomial chaos expansions
Table 1: List of classical univariate polynomial families common in polynomial chaos expan-
sion applications.
Type of variable Distribution Orthogonal polynomials Hilbertian basis ψ
k
(x)
Uniform 1
]−1,1[
(x)/2 Legendre P
k
(x) P
k
(x)/
q
1
2k+1
Gaussian
1
√
2π
e
−x
2
/2
Hermite H
e
k
(x) H
e
k
(x)/
√
k!
Gamma x
a
e
−x
1
R
+
(x) Laguerre L
a
k
(x) L
a
k
(x)/
q
Γ(k+a+1)
k!
Beta 1
]−1,1[
(x)
(1−x)
a
(1+x)
b
B(a) B(b)
Jacobi J
a,b
k
(x) J
a,b
k
(x)/J
a,b,k
J
2
a,b,k
=
2
a+b+1
2k+a+b+1
Γ(k+a+1)Γ(k+b+1)
Γ(k+a+b+1)Γ(k+1)
given in the table are abundant in the literature, see e.g. ?. The values of the polynomials
are computed using Equation 1.9.
Note: The computation of the polynomial values via recurrence relation is not always
stable. In practice, Laguerre and Jacobi polynomials of degree over 23 should be
avoided. See ? for details.
1.3.2 Extension to arbitrary distributions
In case the input random variables are not independent, or no standard polynomials are
defined for their marginal distributions, there are two possibilities for the choice of marginal
distributions: performing an isoprobabilistic transform to a space with independent standard
marginals, or computation of a custom set of polynomials which are orthonormal to the
non-standard distribution.
1.3.2.1 Isoprobabilistic transform
It is possible to define an isoprobabilistic transform from the original probabilistic space to
the so-called reduced space, which has the same dimensionality, but where the input random
variables are independent and have standard marginals. Consider an input vector of random
variables Z with joint PDF Z ∼ f
Z
(z). Then, there exists an isoprobabilistic transform T
such that:
X = T (Z), Z = T
−1
(X) (1.7)
where X is a random vector with independent components distributed according to one of
the distributions in Table 1. We can then compute a PCE in terms of X as in Eq. (1.2) and
rewrite it in terms of Z as follows:
Y = M(Z) =
X
α∈N
M
y
α
Ψ
α
(T (Z)) (1.8)
This transform also allows to use any type of orthonormal polynomial with any type of input
marginals, at the cost of an additional isoprobabilistic transform. Note that this type of
transform can be highly non-linear, more so when transforming from a compact to a non-
compact support distribution (e.g., uniform to Gaussian). The added non-linearity can have
UQ[PY]LAB-V1.0-104 - 3 -
UQ[PY]LAB user manual
a significant detrimental effect on the accuracy of the final truncated PCE, because it may
result in a more complex model.
For simplicity and without loss of generality, we will leave out any reference to possible
isoprobabilistic transforms in the following and consider the random vector X as the vector
of independent input variables with marginal distributions defined in Table 1.
1.3.2.2 Polynomials orthonormal with respect to arbitrary distributions
For a given random variable with probability density function f(x), it is possible to construct
an associated set of univariate polynomials π
n
, n = 0, 1, 2, . . ., which are orthogonal to f(x)
(provided certain mild assumptions on f (x) hold, see ?). The corresponding orthonormal
polynomials are given by ˜π
n
=
π
n
√
⟨π
n
,π
n
⟩
, where ⟨g, h⟩ =
R
g(x)h(x)f(x)dx.
In UQ[PY]LAB, the polynomials are computed through the Stieltjes procedure. This pro-
cedure generates a sequence of univariate polynomials, which are of increasing degree and
orthogonal to the given probability distribution f (x), using the following recurrence relation:
p
β
n+1
˜π
n+1
(x) = (x − α
n
)˜π
n
(x) −
p
β
n
˜π
n−1
(x), n = 0, 1, 2, ... (1.9)
where α
n
and β
n
are defined by the so-called Christoffel-Darboux formulae
α
n
=
⟨xπ
n
, π
n
⟩
⟨π
n
, π
n
⟩
, (1.10)
β
n
=
⟨π
n
, π
n
⟩
⟨π
n−1
, π
n−1
⟩
. (1.11)
The numerical integrations for the inner products involved in the Christoffel-Darboux for-
mulae are performed using the MATLAB adaptive integrator (?). For more details on the
procedure refer to ?.
1.3.3 Basis truncation schemes
Given the polynomials in Table 1, it is straightforward to define the total-degree truncation
scheme, which corresponds to all polynomials in the M input variables of total degree less
than or equal to p:
A
M,p
= {α ∈ N
M
: |α| ≤ p} card A
M,p
≡ P =
M + p
p
(1.12)
Note that the total-degree basis grows exponentially with the degree p.
In many applied science problems, not all terms in the basis are equally important. Often, the
important terms in the expansion tend to be the ones where only few variables are involved.
This is known as the sparsity-of-effects principle. Based on this reasoning, two additional
truncation schemes have been proposed which have the effect of excluding terms of high
interaction order.
UQ[PY]LAB-V1.0-104 - 4 -
Polynomial chaos expansions
1.3.3.1 Restriction of maximum interaction
This truncation scheme is based on choosing a subset of the terms defined in Eq. (1.12), such
that the α’s have at most r non-zero elements (low-rank α):
A
M,p,r
= {α ∈ A
M,p
: ||α||
0
≤ r}, (1.13)
where ||α||
0
=
M
P
i=1
1
{α
i
>0}
is the rank of the multi-index α. For a multivariate polynomial
ψ
α
(x), the rank of α corresponds to the number of variables involved in the polynomial, also
called interaction order.
This truncation scheme can be used to significantly reduce the cardinality of the polynomial
basis by limiting the number of interaction terms, which is particularly effective in high
dimension.
1.3.3.2 Hyperbolic truncation
A modification of the standard scheme, the hyperbolic (or q-norm) truncation scheme makes
use of the so-called q-norm to define the truncation (?):
A
M,p,q
= {α ∈ A
M,p
: ||α||
q
≤ p}, (1.14)
where:
||α||
q
=
M
X
i=1
α
q
i
!
1/q
. (1.15)
Note that for q = 1, hyperbolic truncation corresponds exactly to the standard total-degree
truncation scheme in Eq. (1.12). For q < 1, hyperbolic truncation includes all the univari-
ate high-degree terms, but excludes high-degree terms with many interacting variables. An
example of the behaviour of the hyperbolic norm in two dimensions for different values of p
and q is shown in Figure 1.
1.3.4 Basis-adaptive PCE
In real-world applications, it is often not clear which finite basis would yield the best PCE
approximation. On one hand, the basis must contain enough elements to enable an accurate
representation. On the other hand, the available number of experimental design points limits
the size of the basis.
In this case, a popular strategy is basis-adaptive PCE, where a suitable basis is chosen from a
set of candidate bases. Here, starting from a small candidate basis, one gradually generates
new bases by adding new elements (e.g. by increasing the maximum polynomial degree in
the truncation scheme) and calculates the corresponding PCE and (an approximation to) its
generalization error. Finally, the best PCE in terms of generalization error is chosen. For
any adaptive algorithm, a crucial tool is a-posteriori cross-validation error estimation (see
UQ[PY]LAB-V1.0-104 - 5 -
UQ[PY]LAB user manual
0 5
0
2
4
6
p=3, q=1
0 5
0
2
4
6
p=4, q=1
0 5
0
2
4
6
p=5, q=1
0 5
0
2
4
6
p=6, q=1
0 5
0
2
4
6
p=3, q=0.75
0 5
0
2
4
6
p=4, q=0.75
0 5
0
2
4
6
p=5, q=0.75
0 5
0
2
4
6
p=6, q=0.75
0 5
0
2
4
6
p=3, q=0.5
0 5
0
2
4
6
p=4, q=0.5
0 5
0
2
4
6
p=5, q=0.5
0 5
0
2
4
6
p=6, q=0.5
Figure 1: Hyperbolic truncation set for varying values of p (constant in each column) and q
(constant in each row) as defined in Eq. (1.14). For q = 1, hyperbolic truncation reduces to
the standard total-degree truncation scheme (first row). Decreasing the value of q decreases
the number of polynomials of high interaction order included in the expansion.
Section 1.4) because it allows one to estimate the accuracy of the model without having to
run additional expensive model evaluations to generate a proper validation set.
Two commonly used strategies for basis-adaptive PCE, which can also be combined with each
other, are degree adaptivity and q-norm adaptivity. Given a target accuracy ϵ
T
and a maximum
number of iterations NI
max
, the algorithms can be summarized as follows:
1. Generate an initial basis with one or a combination of the truncation schemes in Sec-
tion 1.3.3, with p = p
0
(q = q
0
);
2. Calculate the PCE coefficients and the corresponding generalization error estimate ϵ
(e.g. ϵ
LOO
in Eqs. (1.19) or (1.22));
3. Compare the error with the preset threshold ϵ
T
. If ϵ ≤ ϵ
T
or the number of iterations
NI = NI
max
, stop the algorithm and return the PCE with the lowest generalization
error. Otherwise, set p = p + 1 (increase q) and return to Step 1.
This simple algorithm is effective in letting the maximum degree or the q-norm of the PCE
be driven directly from the available data. For this type of algorithm to properly converge,
an error estimate sensitive to over-fitting should be chosen, e.g. ϵ
LOO
.
UQ[PY]LAB-V1.0-104 - 6 -
Polynomial chaos expansions
1.4 A posteriori error estimation
After the metamodel is constructed (see Section 1.5), its accuracy and predictive quality can
be quantified by estimating the relative generalization error ϵ
gen
. It is defined as:
ϵ
gen
=
E
h
M(X) − M
P C
(X)
2
i
Var [Y ]
(1.16)
In case next to the training set used to construct the metamodel an independent set of inputs
and outputs, also called validation set, [X
Val
, Y
Val
= M(X
Val
)] is available, the validation
error can be calculated as:
ϵ
Val
=
N − 1
N
N
P
i=1
M(x
(i)
Val
) − M
P C
(x
(i)
Val
)
2
N
P
i=1
M(x
(i)
Val
) − ˆµ
Y
Val
2
(1.17)
where ˆµ
Y
Val
=
1
N
N
P
i=1
M(x
(i)
Val
) is the sample mean of the validation set response. This version of
the relative generalization error is useful to compare the performance of different surrogate
models when evaluated on the same validation set.
If no validation set is available, as commonly the case with expensive computational models,
there are two main ways to estimate ϵ
gen
: normalized empirical error and leave-one-out cross-
validation error.
1.4.1 Normalized empirical error
The normalized empirical error ϵ
emp
is an estimator of the generalization error based on the
accuracy with which the metamodel reproduces the experimental design model evaluations
Y. It is given by:
ϵ
emp
=
N
P
i=1
M(x
(i)
) − M
P C
(x
(i)
)
2
N
P
i=1
M(x
(i)
) − ˆµ
Y
2
(1.18)
where ˆµ
Y
is the sample mean of the experimental design response.
The estimator in Eq. (1.18), although inexpensive to calculate, leads to over-fitting: it is
a monotone decreasing function of the polynomial degree p regardless of the size of the
experimental design.
1.4.2 Leave-one-out cross-validation error
The leave-one-out (LOO) cross-validation error ϵ
LOO
is designed to overcome the over-fitting
limitation of ϵ
emp
by using cross-validation, a technique developed in statistical learning the-
ory. It consists in building N metamodels M
P C\i
, each one created on a reduced experimen-
tal design X \x
(i)
= {x
(j)
, j = 1, ..., N, j ̸= i} and comparing its prediction on the excluded
UQ[PY]LAB-V1.0-104 - 7 -
UQ[PY]LAB user manual
point x
(i)
with the real value y
(i)
(see, e.g., ?). The leave-one-out cross-validation error can
be written as:
ϵ
LOO
=
N
P
i=1
M(x
(i)
) − M
P C\i
(x
(i)
)
2
N
P
i=1
M(x
(i)
) − ˆµ
Y
2
. (1.19)
In practice, when the results of a least-square minimization (see Section 1.5.2) are available,
there is no need to explicitly calculate N separate metamodels, but one can use the following
formulation to calculate ϵ
LOO
(see ?, appendix D):
ϵ
LOO
=
N
P
i=1
M(x
(i)
)−M
P C
(x
(i)
)
1−h
i
2
N
P
i=1
M(x
(i)
) − ˆµ
Y
2
, (1.20)
where h
i
is the i
th
component of the vector given by:
h = diag
A(A
T
A)
−1
A
T
, (1.21)
and A is the experimental matrix in Eq. (1.32).
1.4.3 Corrected error estimates
Further empirical corrections exist to ensure that the generalization error estimates from both
ϵ
emp
and ϵ
LOO
are not underestimated. They generally take the form:
ϵ
∗
= ϵ T (P, N ), (1.22)
where ϵ
∗
is the corrected error, ϵ is the original error and T (P, N) is a correction factor that
typically increases with the number of regressors P and tends to 1 when the size of the
experimental design N → ∞.
A correction factor particularly effective in the context of PCE with small experimental de-
signs (?, after ?) is given by:
T (P, N ) =
N
N − P
1 +
tr(C
−1
emp
)
N
!
, (1.23)
with
C
emp
=
1
N
A
T
A. (1.24)
1.5 Calculation of the coefficients
Several methods exist to calculate the coefficients y
α
of the polynomial chaos expansion for
a given basis (Eq (1.3)). In UQ[PY]LAB, only non-intrusive methods are implemented, i.e.
the coefficients are obtained by post-processing of the experimental design, a set consisting of
UQ[PY]LAB-V1.0-104 - 8 -
Polynomial chaos expansions
samples of the input random variables and the corresponding model evaluations.
The two principal strategies to calculate the polynomial chaos coefficients non-intrusively
are projection and regression. Projection methods use the orthogonality of the basis func-
tions to compute the coefficients by numerical integration (quadrature). Regression methods
formulate Eq. (1.2) as a system of linear equations and solve the system by standard lin-
ear regression approaches. One such regression approach is ordinary least-squares regression
(Section 1.5.2), for which the system of linear equations must be overdetermined. Alterna-
tively, sparse regression techniques can be used, which aim to find coefficient vectors with
only few nonzero entries (i.e., sparse solutions), and which are able to find solutions to un-
derdetermined systems of equations.
1.5.1 Projection method
The calculation of the polynomial coefficients y
α
with the projection approach directly fol-
lows from the definition of PCE given in Eq. (1.2) and from the orthonormality of the poly-
nomial basis. Indeed, taking the expectation value of Eq. (1.2) multiplied by Ψ
β
(x) yields:
y
α
= E [Ψ
α
(X) · M(X)] (1.25)
The calculation of the coefficients is therefore reduced to the calculation of the expectation
value in Eq. (1.25). It can be cast as a numerical integration problem which in turn can be
efficiently solved using quadrature methods.
Gaussian quadrature: a standard tool in the numerical evaluation of integrals, Gaussian
quadrature is based on a simple weighted-sum scheme:
y
α
=
Z
Ω
X
M(x)Ψ
α
(x)f
X
(x)dx ≈
N
X
i=1
w
(i)
M(x
(i)
)Ψ
α
(x
(i)
) (1.26)
The set of weights w
(i)
and quadrature points x
(i)
(the experimental design) are derived
from Lagrange polynomial interpolation and guarantees exactness in the evaluation of the
integrals of functions of polynomial complexity (?).
The integration weights w
(i)
and the integration nodes x
(i)
are uniquely determined by the
marginals of the independent components of the input random vector X, and they corre-
spond to the roots of the corresponding polynomial basis functions as reported in Table 1. In
UQ[PY]LAB the Gaussian quadrature nodes and weights are numerically calculated with the
Golub-Welsch algorithm (?, ?).
Standard multivariate Gaussian quadrature is achieved by a tensor-product of univariate in-
tegration rules. Therefore the number of integration nodes (i.e. full-model evaluations)
increases rapidly with the number of input variables. As an example, selecting a max poly-
nomial degree p would require (p + 1) integration points in each dimension, leading to
N = (p + 1)
M
in Eq (1.26). This is the so-called curse of dimensionality.
UQ[PY]LAB-V1.0-104 - 9 -
UQ[PY]LAB user manual
Sparse quadrature: a more recent tool to deal with high-dimensional integration, Smolyak’
sparse quadrature is an alternative approach to the original tensor-product multi-dimensional
quadrature (?). The integration is still performed as given in Eq. (1.26), but the weights are
derived from a combination of lower order standard quadrature terms:
Q
M,l
Smolyak
≡
X
l+1⩽|i|⩽l+M
(−1)
M+l−|i|
·
M − 1
k + M − |i|
!
· Q
i
(1.27)
where:
i = i
1
, i
2
, ..., i
M
, |i| = i
1
+ ... + i
M
∈ N
and
Q
i
= Q
i
1
⊗ .... ⊗ Q
i
M
is a tensor product of the lower order Gaussian quadrature rules identified by the multi-
index i. This method can lead to a substantially reduced number of integration points w.r.t.
classical Gaussian quadrature without sacrificing accuracy in higher dimensions.
1.5.1.1 Error estimation
An a posteriori error estimate of the Gaussian quadrature error in the estimation of the PCE
coefficients in Eq. (1.26) can be calculated by taking the expectation value of the residual
mean-square error E
M(X) − M
P C
(X)
by integrating it with the same quadrature rule
and on the same nodes:
ϵ
res
≈
N
P
i=1
w
(i)
Y
(i)
− y
T
Ψ(x
(i)
)
2
N
P
i=1
M(x
(i)
) − ˆµ
Y
2
(1.28)
where y = {y
α
1
, . . . , y
α
P
}
T
is the vector of polynomial coefficients,
Ψ(x
(i)
) =
Ψ
α
1
(x
(i)
), . . . , Ψ
α
P
(x
(i)
)
T
is a vector containing the values of the polynomial
basis elements at quadrature point x
(i)
and ˆµ
Y
=
1
N
N
P
i=1
M(x
(i)
) is the sample mean of the
set of quadrature points.
1.5.2 Least-Squares regression
A different approach to estimate the coefficients in Eq. (1.3) is to set up a least-squares
minimization problem (?). The infinite series in Eq. (1.2) can be written as a sum of its
truncated version Eq. (1.3) and a residual:
Y = M(X) =
P −1
X
j=0
y
j
Ψ
j
(X) + ε
P
≡ y
T
Ψ(X) + ε
P
(1.29)
where P = card A
M,p
, ϵ
P
is the truncation error, y
α
= {y
0
, . . . , y
P −1
}
T
is a vector contain-
ing the coefficients and Ψ(X) = {Ψ
0
(X), . . . , Ψ
P −1
(X)}
T
is the vector that assembles the
UQ[PY]LAB-V1.0-104 - 10 -
Polynomial chaos expansions
values of all the orthonormal polynomials in X.
The least-square minimization problem can then be set-up as:
ˆ
y = arg min E
y
T
Ψ(X) − M(X)
2
. (1.30)
1.5.2.1 Ordinary Least-Squares
A direct approach to solving Eq. (1.30) is given by Ordinary Least-Squares (OLS). Given a
sample X = {x
(1)
, ..., x
(N)
}
T
of size N of the input random vector X (the experimental de-
sign) and the corresponding model responses Y = {y
(1)
, ..., y
(N)
}
T
, the ordinary least-squares
solution of Eq. (1.30) reads:
ˆ
y = (A
T
A)
−1
A
T
Y, (1.31)
where
A
ij
= Ψ
j
x
(i)
i = 1, . . . , n ; j = 0, . . . , P − 1 (1.32)
is the so-called experimental matrix (or regression matrix) that contains the values of all the
basis polynomials in the experimental design points.
The main advantage of the least-squares minimization method over the projection method
lies in the fact that an arbitrary number of points can be used to calculate the coefficients,
as long as they are a representative sample of the random input vector X and N ≥ P .
A theoretical analysis of the convergence of the least-squares minimization method can be
found in ?.
1.5.3 Sparse PCE: Least Angle Regression
A complementary strategy to favour sparsity in high dimension consists in directly modifying
the least-square minimization problem in Eq. (1.30) by adding a penalty term of the form
λ||y||
1
, i.e. solving:
ˆy = argmin
y∈R
|A|
E
y
T
Ψ(X) − M(X)
2
+ λ||y||
1
, (1.33)
where the regularization term ||ˆy||
1
=
P
α∈A
|y
α
| forces the minimization to favour low-
rank solutions. Several algorithms exist that solve the penalized minimization problem in
Eq. (1.33), including least absolute shrinkage and selection operator (LASSO, ?), forward
stagewise regression (?) and least angle regression, or LAR, (?). In the context of PCE, ?
successfully applied the LAR algorithm to obtain sparse PCE models that are accurate even
with very small experimental designs.
1.5.3.1 The LAR algorithm
The LAR algorithm is a linear regression tool based on iteratively moving regressors from a
candidate set to an active set. The next regressor is chosen based on its correlation with the
current residual. At each iteration, analytical relations are used to identify the best set of
UQ[PY]LAB-V1.0-104 - 11 -
UQ[PY]LAB user manual
regression coefficients for that particular active set, by imposing that every active regressor
is equicorrelated with the current residual.
The optimal number of predictors in the metamodel (i.e. the optimal number of LAR steps)
may be determined using a suitable criterion.
The full LAR algorithm in the context of PCE (?) reads:
Initialization
• y
α
= 0, ∀α ∈ A
M,p,q
;
• Candidate set: Ψ
α
;
• Active set: ∅;
• Residual: r
0
= Y
Iterative algorithm:
1. Find the regressor Ψ
α
j
that is most correlated with the current residual
2. Move all the coefficients of the current active set towards their least-square value until
their regressors are equicorrelated to the residual as some other regressor in the can-
didate set. This regressor will also be the most correlated to the residual in the next
iteration.
3. Calculate and store the error estimate ϵ
j
LOO
for the current iteration
4. Update all the active coefficients and move Ψ
α
j
from the candidate set to the active set
5. Repeat the previous steps until the size of the active set is equal to m = min(P, N − 1)
After the iterations are finished, the candidate set of regressors with the lowest ϵ
LOO
is se-
lected as the best sparse candidate basis. A typical example of the evolution of ϵ
j
LOO
vs. j is
shown in Figure 2.
1.5.3.2 Hybrid LAR
One limitation of the LAR algorithm is that it is defined only for non-constant regressors (due
to the presence of the cross-correlation-based selection at the first step of the algorithm).
Moreover, the hyperparameter of the algorithm is the number of LAR steps. Therefore, the
leave-one-out error can not be calculated as easily as OLS that does not require rebuilding
N models (see Eq. (1.20)). Both limitations can be overcome by introducing the so-called
hybrid-LAR step. At the end of each LAR iterations, the constant regressor is added to the
selected basis, and OLS is performed to calculate the related coefficients as well as the corre-
sponding leave-one-out error. Thus, in this setting, LAR iterations are used to provide a series
of set of basis functions, and OLS is used to build the associated surrogate models. The final
model is chosen as the one having the smallest leave-one-out error.
UQ[PY]LAB-V1.0-104 - 12 -
Polynomial chaos expansions
Figure 2: Typical evolution of ϵ
LOO
vs. LAR iteration. The evolution is in many cases smooth
and consistently convex throughout the iterations. In some cases with very small experimen-
tal designs, however, local minima in the early iterations can be observed.
1.5.3.3 LAR early stop criterion
When dealing with a large number of regressors, each LAR iteration can be time consuming
(mostly due to the calculation of the ϵ
LOO
, which entails a large matrix inversion). An early
stop criterion can be introduced in practical implementations to mitigate the corresponding
costs and speed-up the algorithm. The criterion stems from the observation that in most
real-case scenarios the behaviour of ϵ
j
LOO
is relatively smooth with the iteration number j
(see Figure 2) and convex. An effective and robust early stop criterion for LAR is then to stop
adding regressors after the ϵ
LOO
is above its minimum value for at least 10% of the maximum
number of possible iterations.
1.5.4 Sparse PCE: Orthogonal Matching Pursuit
Orthogonal Matching Pursuit (OMP) is a greedy algorithm proposed by ? as a refinement
to the Matching Pursuit algorithm by ?. OMP works by iteratively retrieving the polynomial
basis elements that are most correlated with the current approximation residual and adding
them to the active set of regressors.
OMP uses a greedy iterative strategy that minimizes the approximation residual at each iter-
ation. Consider the approximation residual R
n
for a polynomial basis with n elements and
the approximation residual R
n+1
for a polynomial basis with n + 1 elements. By projecting
the residual R
n
onto a new polynomial the following equation is obtained:
R
n
=
R
n
, Ψ
α
n+1
Ψ
α
n+1
+ R
n+1
(1.34)
Therefore, the residual R
n+1
is by construction orthogonal to the new polynomial Ψ
α
n+1
.
UQ[PY]LAB-V1.0-104 - 13 -
UQ[PY]LAB user manual
Projecting Eq. (1.34) onto R
n
yields:
R
n
2
=
⟨R
n
, Ψ
α
n+1
⟩
2
+
R
n+1
2
. (1.35)
It follows that minimizing the approximation residual R
n+1
is equivalent to choosing a poly-
nomial such that
⟨R
n
, Ψ
α
n+1
⟩
is maximized. Therefore, each iteration of the OMP algorithm
consists in solving the following problem:
Ψ
α
n+1
= arg max
α∈A
⟨R
n
, Ψ
α
⟩
(1.36)
After the basis element Ψ
α
n+1
has been added to the active set of regressors, all the corre-
sponding polynomial coefficients y
α
are updated via ordinary least squares. This additional
step guarantees that the newly calculated residual is orthogonal to all the regressors in the
current active set.
1.5.4.1 The OMP algorithm
The OMP algorithm is a linear regression tool that minimizes the norm of the approximation
residual at each iteration. The algorithm uses the leave-one-out error estimator in Eq. (1.19)
to adaptively select the best active set. The implementation of the OMP algorithm in the
context of PCE reads:
Initialization:
• y
0
α
= 0, ∀α ∈ A
M,p,q
;
• Candidate set: Ψ
C,0
= Ψ
α
;
• Active set: Ψ
A,0
= ∅;
• Residual: R
0
= Y
Iterative algorithm:
1. Find the polynomial Ψ
α
j
that is most correlated with the current approximation resid-
ual R
j−1
.
2. Add the polynomial Ψ
α
j
to the active set of polynomials, i.e. Ψ
A,j
= Ψ
A,j−1
∪ Ψ
α
j
.
3. Calculate the new polynomial coefficients y
j
α
by projecting the model response Y onto
the active set of polynomials, i.e. calculate an ordinary least squares using the active
set.
4. Calculate the new approximation residual R
j
= Y − Ψ
A,j
y
j
α
.
5. Calculate and store the error estimate ϵ
j
LOO
for the current iteration.
6. Repeat the previous steps until the size of the active set is equal to m = min (P, N).
After the iterative procedure is terminated, the active set of polynomials with the lowest ϵ
LOO
is selected as the best sparse basis.
UQ[PY]LAB-V1.0-104 - 14 -
Polynomial chaos expansions
1.5.4.2 OMP early stop criterion
When the polynomial basis size P is large, each OMP iteration can be computationally ex-
pensive. An early stop criterion can be used to reduce the computational costs. Similar to
the LAR algorithm, the behaviour of ϵ
j
LOO
is relatively smooth and convex (see Figure 3).
The proposed early stop criterion for OMP stops adding regressors after the ϵ
LOO
is above its
minimum value for at least 10% of the maximum number of possible iterations.
Figure 3: Typical evolution of ϵ
LOO
vs. OMP iterations. The evolution is in many cases
relatively smooth and consistently convex throughout the iterations. In some cases with very
small experimental designs, however, local minima in the early iterations can be observed.
1.5.5 Sparse PCE: Subspace Pursuit
Subspace pursuit (SP) is a sparse regression algorithm developed by ? and introduced for
PCE by ?. Its single hyperparameter is the number K of nonzero elements in the coeffi-
cient vector. Starting from an initial set of K regressors that are highly correlated with
the model evaluations, it computes the corresponding coefficients by ordinary least-squares
(OLS), identifies K more regressors that are highly correlated with the residual, computes
the OLS solution on the combined set of 2K regressors, and removes again the K regressors
with the smallest-in-magnitude coefficients. This is repeated until convergence.
Theoretical guarantees for this algorithm were derived by ?. Numerically, it has been found
that SP performs particularly well for low-dimensional problems with moderate to large ex-
perimental design sizes (?).
UQ[PY]LAB-V1.0-104 - 15 -
UQ[PY]LAB user manual
1.5.5.1 The SP algorithm
Denote by A
†
= (A
T
A)
−1
A
T
the pseudoinverse of a matrix A. Then, the OLS solution y to
an overdetermined system Ay = Y can be written as y = A
†
Y.
Initialization:
• Given K, the desired number of nonzero coefficients (2K ≤ min{N, P })
• Given the set of candidate regressors {ψ
α
: α ∈ A} and the associated regression
matrix Ψ
• Initial active set: A
0
= {K indices corresponding to the largest-in-magnitude entries in Ψ
T
Y}
• Denote by Ψ
A
0
the submatrix corresponding to A
0
• Residual vector: r
0
= Y − Ψ
A
0
(Ψ
A
0
)
†
Y
Iterative algorithm: In step j = 1, 2, . . .,
1. Define the temporary set
S = A
j−1
∪ {K indices corresponding to the largest-in-magnitude entries in Ψ
T
r
j−1
}
containing 2K regressors.
2. Compute the corresponding coefficients by OLS: y = (Ψ
S
)
†
Y.
3. Update the active set
A
j
= {K indices corresponding to the largest-in-magnitude entries in y}.
4. Update residual: r
j
= Y − Ψ
A
j
(Ψ
A
j
)
†
Y.
5. If ∥r
j
∥
2
≥ ∥r
j−1
∥
2
, stop the iteration and return A
j−1
and the associated coefficient
vector and leave-one-out error. Otherwise, set j = j + 1 and continue with step 1.
The hyperparameter K of SP can be determined, e.g., by leave-one-out cross-validation
(LOO). In this case, the above algorithm is run for a range of reasonable values for K. The
final solution is the one with the lowest LOO error. Note that the LOO error can be computed
inexpensively as in Eq. (1.20), since SP uses OLS on each subset of regressors.
1.5.6 Sparse PCE: Bayesian compressive sensing
A wholly different approach to the sparse regression problem is Bayesian compressive sens-
ing (BCS), where the regression problem is rewritten in a Bayesian framework. Likelihood,
priors, and auxiliary parameters are chosen to encourage sparsity in the regression coeffi-
cients. The coefficients are computed to be the maximum-a-posteriori (MAP) estimate given
the model evaluations.
One of the first publications on BCS was the relevance vector machine (RVM) by ?. ? pro-
posed the use of BCS in the context of PCE. Here, we use the BCS formulation of ? (Fast
Laplace), which is a generalization of RVM attaining a provably sparser solutions. An illus-
tration of the BCS framework is shown in Fig. 4.
UQ[PY]LAB-V1.0-104 - 16 -
Polynomial chaos expansions
Each of the quantities and auxiliary variables is considered to be random with a certain
parametrized distribution, except for σ
2
, which is chosen beforehand (see below), and ν = 0.
More precisely, the assumptions on parameters and distributions are as follows. The model
evaluations Y are considered as i.i.d. realizations of a random variable parametrized by coef-
ficients y and a noise variance parameter σ
2
: p(Y|y, σ
2
) = N(Y|Ψy, σ
2
1
N
) (likelihood).
The noise variance is a given (fixed) hyperparameter of the algorithm, while the coeffi-
cients y are again random variables, each normally distributed with its own variance γ
i
:
p(y
i
|γ
i
) = N(y
i
|0, γ
i
). The γ
i
’s are i.i.d. following an exponential distribution with shared
parameter λ: p(γ
i
|λ) = Exp(γ
i
|
λ
2
). Finally, λ is the last layer of the hierarchical framework,
drawn from a Gamma-distribution p(λ|ν) = Γ(λ|
ν
2
,
ν
2
) with ν → 0 (resulting in an uninfor-
mative, improper prior p(λ) ∝
1
|λ|
).
The goal of the algorithm is to compute the posterior distribution p(y, γ, λ|Y) for a set of
given model evaluations Y. From this, the MAP estimate for y can be determined. However,
the prior distributions of BCS are chosen to encourage sparsity, but are not conjugate, so
that an analytical solution of the problem is not feasible. Instead, a fast approximate algo-
rithm is employed, which iteratively updates the different auxiliary quantities as well as the
coefficient vector estimate y.
It has been found that BCS performs especially well for high-dimensional problems and in
cases where only a small set of data is available (?).
γ
y
σ
2
Y
p(y|γ) =
Q
i
N(y
i
|0, γ
i
)
p(Y|y, σ
2
) = N(Y|Ψy, σ
2
1)
λ
p(γ
i
|λ) =
λ
2
exp
−
λ
2
γ
i
p(λ|ν) = Γ
λ|
ν
2
,
ν
2
ν
Figure 4: Illustration of the BCS setup by ?. All quantities are random variables with distribu-
tions parametrized by hyperparameters. The hierarchical structure and the choice of priors
and hyperpriors results in a sparsity-encouraging posterior distribution for the coefficients y.
1.5.6.1 The BCS algorithm
The idea of the BCS algorithm as proposed by ? (called Fast Laplace) is to maintain a vector
γ of coefficient variances. If γ
i
= 0, the associated coefficient y
i
must be zero as well, and
the corresponding term in the basis is inactive. Else, if γ
i
> 0, the term is active. In each
iteration, one regressor is chosen, and to be either added to the basis (γ
i
is set to a value > 0),
deleted from the basis (γ
i
= 0), or reassessed by means of a re-estimation of its variance γ
i
.
UQ[PY]LAB-V1.0-104 - 17 -
UQ[PY]LAB user manual
The objective function L is the logarithm of p(γ, λ, Y), which has an analytical expression
and can be differentiated to obtain the update equations. We refer to ? and (?, Appendix A)
for further details.
Initialization:
• Given model evaluations Y and the corresponding matrix of regressor evaluations Ψ
• Given the noise variance σ
2
• Given stopping threshold η
• Initialize with the constant regressor as the only regressor in the model
Iterative algorithm:
1. For each of the regressors, compute the hypothetical updated value of γ
i
if this regres-
sor alone was updated, and the associated hypothetical change in L(γ). Choose the
regressor i that accounts for the largest increase in L.
2. If none of the regressors results in an increase in L, or if the increase in L divided by
the current value of L has been smaller than η twice in a row, stop.
3. Add/remove/re-estimate: depending on the old and the new value of γ
i
, the change
corresponds to either removal, addition, or re-estimation of this regressor. All other
quantities have to be updated accordingly. Continue with step 1.
The hyperparameter σ
2
can be chosen by k-fold cross-validation as follows. A number of
candidate values for σ
2
is generated. The available data Y is divided into k equally sized
chunks. Each chunk is in turn treated as validation set, while the remaining data is used
to train the model. Averaging over the k validation errors gives the corresponding cross-
validation error. This is done for each value of σ
2
. The σ
2
with the smallest cross-validation
error is chosen. Finally, BCS is run one more time with the chosen value of σ
2
and the full
experimental design to get the final solution.
Throughout this document, algorithms like degree-adaptivity (Section 1.3.4) are described
using the LOO error ϵ
LOO
, but hold equivalently for BCS by replacing the LOO error with the
k-fold CV error ϵ
CV
.
1.6 Post-processing
The polynomial chaos expansion can be post-processed to obtain further information about
the quantities of interest.
1.6.1 Moments of a PCE
Due to the orthonormality of the polynomial basis, the first two moments of a PCE are en-
coded in its coefficients. In particular, the mean value of a PCE reads:
µ
P C
= E
M
P C
(X)
= y
0
(1.37)
UQ[PY]LAB-V1.0-104 - 18 -
Polynomial chaos expansions
where y
0
is the coefficient of the constant basis term Ψ
0
= 1.
Similarly, the variance of a PCE reads:
(σ
P C
)
2
= E
(M
P C
(X) − µ
P C
)
2
=
X
α∈A
α̸=0
y
2
α
(1.38)
where the summation is over the coefficients of the non-constant basis elements only.
1.6.2 Uncertainty propagation and sensitivity analysis
When the polynomial coefficients are known, it is straightforward to evaluate the metamodel
on new samples of the input random vector X. In fact, it is sufficient to directly apply
Eq. (1.2) by first evaluating the multivariate polynomials on the new sample and summing
them weighted by their coefficients. The computational costs to perform this operation are
limited to the evaluation of the univariate polynomials on the new sample and a small num-
ber of matrix multiplications, hence making this operation very efficient. This property can
be used effectively for calculating the PDF of the model response accurately by using large
Monte-Carlo samples of the inputs and, e.g., Kernel smoothing techniques.
Another important property of PCE is that the coefficients encode important information
about the ANOVA decomposition of the surrogate model, which can be exploited to effectively
calculate global sensitivity indices at very limited costs. The reader is referred to the UQLAB
User Manual – Sensitivity Analysis module for further information on the relation between
global sensitivity analysis and polynomial chaos expansions.
1.7 Bootstrap-based error estimates
The PCE module of UQ[PY]LAB offers the possibility to estimate the accuracy of a local
prediction through bootstrap resampling.
1.7.1 Bootstrap PCE
? have shown how the bootstrap resampling technique (?) can be applied to PCE to provide
a local error estimator, calling this method bPCE. Resampling with substitution is used on
the experimental design X , thus generating a set of bootstrap replications. Each replication
contains the same number of sample points as the original experimental design X . Each of
the B replications {X
(1)
, . . . , X
(B)
} is used to calculate a corresponding PCE. This yields
a set of B different PCEs with coefficients y
(b)
α
, b = 1, . . . , B and, consequently, a set of B
responses at each prediction point x. Such sets of responses can be used to calculate the
local variance (or other statistics of interest such as quantiles) of the PCE predictor, due to
the finite size of the experimental design. An example of bPCE on a 1-dimensional function is
illustrated in Figure 5, where the replications are represented as thin orange lines, and 95%
confidence bounds are derived from empirical quantiles on the bootstrap samples.
UQ[PY]LAB-V1.0-104 - 19 -
UQ[PY]LAB user manual
Figure 5: Function y = x sin(x) approximated by a PCE with confidence bounds from boot-
strapping.
1.7.2 Fast bPCE
In practice, setting up a PCE for each bootstrap sample can be time consuming, especially
when using sparse expansion techniques in high dimension or when an already large exper-
imental design is available. Since for most applications local error estimates do not need to
be very accurate, fast bPCE is used instead in UQ[PY]LAB. In this approach, a sparse PCE
basis is established only once using a sparse regression algorithm (e.g., LAR or OMP) on the
original experimental design X . The coefficients y
(b)
α
, b = 1, . . . , B are then calculated by
OLS regression for each bootstrap sample.
UQ[PY]LAB-V1.0-104 - 20 -
Chapter 2
Usage
In this section a reference problem will be set up to showcase how each of the techniques
described in Part 1 can be deployed in UQ[PY]LAB.
2.1 Reference problem: the Ishigami function
Polynomial chaos expansion aims at approximating a computational model with a polyno-
mial surrogate. To this end, we will make use of a well-known benchmark for polynomial
chaos expansions: the Ishigami function (?, www.sfu.ca/
˜
ssurjano/ishigami.html).
It is an analytical 3-dimensional function characterized by non-monotonicity and high non-
linearity, given by the following equation:
f(x) = sin(x
1
) + a sin
2
(x
2
) + bx
4
3
sin(x
1
) (2.1)
where the parameters are set to a = 7 and b = 0.07 in this example (see e.g. ?).
The input random vector consists of three i.i.d. uniform random variables X
i
∼ U(−π, π). An
example UQ[PY]LAB script that showcases several of the currently available PCE expansion
techniques on the Ishigami function can be found in the example file:
1-Calculation strategies
2.2 Problem set-up
Recalling Eq. (1.3), truncated PCE reads:
M(X) ≈
X
α∈A
y
α
Ψ
α
(X)
The main ingredients that need to be set-up in a PCE analysis are:
• A model to surrogate Y = M(X);
• A probabilistic input model (random input vector X);
• A truncated polynomial basis defined by A;
21
UQ[PY]LAB user manual
• The polynomial coefficients {y
α
, α ∈ A}.
2.2.1 Full model and probabilistic input model
The model in Eq. (2.1) is implemented as a Python file in ishigami.py.
To surrogate it using UQ[PY]LAB, we need to first configure a basic MODEL object:
MetaOpts = {'Type': 'Model',
'ModelFun': 'ishigami.model'}
myModel = uq.createModel(MetaOpts)
For more details about the configuration options available for a model, refer to the UQ[PY]LAB
User Manual – the MODEL module.
The three independent input variables can be defined as:
IOpts = {'Marginals': uq.Marginals(3, 'Uniform', Parameters=[-np.pi, np.pi])}
myInput = uq.createInput(IOpts)
For more details about the configuration options available for an INPUT object, refer to the
UQ[PY]LAB User Manual – the INPUT module.
2.3 Set-up of the polynomial chaos expansion
The PCE module creates a MODEL object that can be used as any other model. Its configura-
tion options, however, are generally more complex than for a basic full model definition like
the one in Section 2.2.1.
The basic options common to any PCE metamodelling MODEL read:
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
The input dimension of the problem M is automatically retrieved by the configuration of the
INPUT module given above. The additional configuration options needed to properly create a
PCE object in UQ[PY]LAB are given in the following sections.
2.4 Orthogonal polynomial basis
2.4.1 Univariate polynomial types
For most practical applications, it is not necessary in UQ[PY]LAB to manually specify the uni-
variate polynomial types that form the multivariate polynomial basis Ψ
α
(X) via Eq. (1.5).
The default behaviour of UQ[PY]LAB is to choose univariate polynomials depending on the
distributions of the input variables, according to Table 2. Internally, a linear isoprobabilistic
transform (see Section 1.3.2.1) to a distribution from the same family with standard param-
eters is automatically performed (e.g., U(a, b) is transformed to U(0, 1)). For non-classical
UQ[PY]LAB-V1.0-104 - 22 -
Polynomial chaos expansions
input distributions, the univariate orthonormal polynomials are computed numerically by
applying the Stieltjes procedure (see Section 1.3.2.2).
Table 2: Default univariate polynomial types used in UQ[PY]LAB w.r.t. input distributions
Input PDF f
X
i
(X
i
) Univariate polynomial family
U(a, b) Legendre
N(µ, σ) Hermite
Γ (λ, κ) Laguerre(λ, κ)
B(r, s, a, b) Jacobi(r, s)
log N(µ, σ) Hermite
Other Arbitrary
Note that the Jacobi and Laguerre polynomials are defined parametrically with the Beta and
Gamma distribution parameters respectively. When a distribution does not belong to the
above, the recurrence terms will be numerically computed if the integral of the distribu-
tion itself can be estimated accurately with numerical integration. Otherwise the Hermite
polynomials will be used for distributions that do not have bounded support and Legendre
polynomials when distributions have bounded support.
It is possible, however, to manually force the univariate polynomial families to the desired
value by specifying the PolyTypes option in the input. As an example, to force the use
of Hermite polynomials in the first dimension, Legendre polynomials in the second, and
numerically compute the orthonormal polynomials for the third direction one has to specify:
MetaOpts['PolyTypes'] = ['Hermite', 'Legendre', 'Arbitrary']
In case of constant input variables, the corresponding PolyTypes entry would be 'Zero'.
In case Laguerre or Jacobi polynomials are selected with PolyTypes, it is necessary for the
PolyTypesParams option to be defined. The definition of the parameters of the polynomial
families are consistent with that of their respective distributions and the redundant param-
eters, such as the bounds of the beta distribution for Jacobi polynomials, are ignored. For
example one can specify:
MetaOpts['PolyTypes'] = ['Hermite', 'Jacobi', 'Laguerre']
MetaOpts['PolyTypesParams'] = [[], [2, 3, 0, 1], [3, 4]]
The dimension of the MetaOpts['PolyTypes'] list must agree with that of the input model.
2.4.2 Truncation schemes
The default truncation strategy in UQ[PY]LAB is the standard total-degree truncation scheme
in Eq. (1.12), with maximum degree p = 3. To specify a desired maximum polynomial degree
it is sufficient to add the Degree key to the MetaOpts configuration variable. To specify a
maximum polynomial degree of e.g. 10 one can add:
MetaOpts['Degree'] = 10
UQ[PY]LAB-V1.0-104 - 23 -
UQ[PY]LAB user manual
2.4.2.1 Basis truncation
Additionally, one can configure any of the truncation strategies described in Section 1.3.3.1
and 1.3.3.2 with the optional TruncOptions key .
A q-norm truncation with q = 0.75 and maximum rank r = 2 can be specified as follows:
MetaOpts['TruncOptions'] = {
'qNorm': 0.75,
'MaxInteraction': 2
}
The two truncation schemes are not mutually exclusive and they can be specified either one-
at-a-time or both together.
2.4.2.2 User-specified basis
It is also possible to directly specify the set of multi-indices A that will be used to gener-
ate the multivariate polynomial basis. This can be accomplished by manually specifying
the P × M matrix of polynomial degrees in the MetaOpts['TruncOptions']['Custom']
key. As an example, one can specify a basis with M = 3, p = 2 and basis elements
Ψ
0,0,0
(x), Ψ
0,2,0
(x), Ψ
1,0,0
(x) and Ψ
0,1,1
(x), as follows:
MetaOpts['TruncOptions'] = {
'Custom': [[0, 0, 0], [0, 2, 0], [1, 0, 0], [0, 1, 1]]
}
In this case, it is not necessary to specify the degree of the basis. It is computed automatically
from the user-specified basis.
2.5 Experimental design
For projection-based PCE, the experimental design is determined by the choice of quadrature
scheme and degree (see Section 2.6.1).
For regression methods, there is more freedom in the choice of experimental design. At least
the number of experimental design points has to be specified:
MetaOpts['ExpDesign'] = {
'NSamples': 500
}
By default, samples are generated by Latin Hypercube sampling (LHS).
Several other options are available for the creation of the experimental design. A summary
of the most common is given in the following.
• Specify a sampling strategy: it is possible to specify another sampling strategy by
adding a Sampling option. The following specifies sampling from a Sobol’ pseudoran-
dom sequence:
MetaOpts['ExpDesign']['Sampling'] = 'Sobol'
UQ[PY]LAB-V1.0-104 - 24 -
Polynomial chaos expansions
• Manually specify an experimental design: it is common to create PCE from already
existing data. There are two ways to import data in a UQ[PY]LAB PCE MODEL:
– Specify the values of X and Y directly. Assuming the existing data is stored in the
local variables X_ED and Y_ED, where each row of these variables corresponds to
one point from the experimental design, and X_ED has as many columns as there
are input random variables, they can be imported in UQ[PY]LAB as follows:
MetaOpts['ExpDesign'] = {
'X': X_ED.tolist(),
'Y': Y_ED.tolist()
}
Note: When an experimental design is specified manually, there is no need to
create a MODEL object as in Section 2.2. However, an INPUT module with
an input random vector compatible with the provided experimental design
must be defined. This is an intrinsic property of PCE: f
X
is needed to
calculate the PCE coefficients.
A comprehensive list of the options available for the calculation of the experimental design
of a PCE can be found in Table 11.
2.6 Calculation of the coefficients
The remaining ingredient needed to complete the PCE is the set of polynomial coefficients
y
α
. In this section, the techniques introduced in Section 1.5 are deployed in UQ[PY]LAB.
2.6.1 Projection: Gaussian quadrature
Calculating the PCE coefficients with Gaussian quadrature does not require any special con-
figuration. Due to the very high non-linearity of the Ishigami function, a relatively high
polynomial degree of p = 14 is needed to achieve satisfactory accuracy.
A projection-based PCE can be created with the following lines of code (note that for quadrature-
based calculation of the coefficients no truncation scheme is necessary, as all the coefficients
up to the specified degree are calculated simultaneously anyway):
# Reporting the previous configuration options as a reminder
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
# Specification of 14th degree, Gaussian quadrature-based projection
MetaOpts['Degree'] = 14
MetaOpts['Method'] = 'quadrature'
# Creation of the metamodel:
myPCE_Quadrature = uq.createModel(MetaOpts)
UQ[PY]LAB-V1.0-104 - 25 -
UQ[PY]LAB user manual
Once the model is created, a report with basic information about the PCE can be printed out
as follows:
uq.print(myPCE_Quadrature)
which produces the following output:
%------------ Polynomial chaos output ------------%
Number of input variables: 3
Maximal degree: 14
q-norm: 1.00
Size of full basis: 680
Size of sparse basis: 661
Full model evaluations: 3375
Quadrature error: 3.8654334e-13
Mean value: 3.5000
Standard deviation: 3.7208
Coef. of variation: 106.309%
%--------------------------------------------------%
Similarly, a visual representation of the spectrum of the resulting non-zero coefficients can
be visualized graphically as follows:
uq.display(myPCE_Quadrature)
which produces the image in Figure 6.
0 200 400 600
−20
−15
−10
−5
0
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 661
𝛼
log
10
(|𝑦
𝛼
|)
Figure 6: Graphical representation of the logarithmic spectrum of the PCE coefficients com-
puted by quadrature. Most of the coefficients of the 661 basis elements are close to 0.
2.6.1.1 Accessing the results
Coefficients and basis
UQ[PY]LAB-V1.0-104 - 26 -
Polynomial chaos expansions
All the information needed to evaluate Eq. (1.3) are available in the dictionary
myPCE_Quadrature['PCE'] having the following keys:
{
'Basis': {'PolyTypes': ['Legendre', 'Legendre', 'Legendre'], ...},
'Coefficients': [3.499999999999992, -1.1006820455072841e-15, ...],
'Moments': {'Mean': 3.499999999999992, 'Var': 13.844587949669927}
}
The myPCE_Quadrature['PCE']['Coefficients'] list contains the coefficients for all of
the 680 Ψ
α
(X) basis elements. The corresponding basis elements are given by the
myPCE_Quadrature['PCE']['Basis'] dictionary having the following keys:
{
'PolyTypes': ['Legendre', 'Legendre', 'Legendre'],
'PolyTypesParams': [],
'PolyTypesAB': [['Empty'], ['Empty'], ['Empty']],
'Indices': [[0, 0, 0], [0, 0, 1], [0, 1, 0], [1, 0, 0], ...],
'MaxCompDeg': [14, 14, 14],
'MaxInteractions': 3,
'Degree': 14,
'qNorm': 1
}
The myPCE_Quadrature['PCE']['Basis']['PolyTypes'] list contains the M univariate
polynomial families ϕ
(i)
in Eq. (1.5). The myPCE_Quadrature['PCE']['Basis']['Indices']
list of lists contains the α multi-indices in Eq. (1.3) in row-vector format (in other words, the
index set in A
M,p,q
set in Eq. (1.14)). To each nested list of
myPCE_Quadrature['PCE']['Basis']['Indices'] corresponds a coefficient in the list
myPCE_Quadrature['PCE']['Coefficients'].
The additional keys contain respectively:
• MaxCompDeg: the maximum univariate polynomial degree for each input variable for
the basis elements with non-zero coefficients.
• MaxInteractions: the maximum rank of the basis elements with non-zero coeffi-
cients.
• Degree: the maximum degree of the basis elements with non-zero coefficients.
Finally, the myPCE_Quadrature['PCE']['Moments'] dictionary contains mean and vari-
ance of the model as calculated from the PCE (Eqs. (1.37),(1.38)).
For more details about the available information in the PCE output, refer to Section 3.2.1.
UQ[PY]LAB-V1.0-104 - 27 -
UQ[PY]LAB user manual
Model evaluations
The quadrature points and their corresponding model evaluations are stored in the dictionary
myPCE_Quadrature['ExpDesign']:
{
'Sampling': 'Quadrature',
'NSamples': 3375,
'X': [[-3.103870036414838, -3.103870036414838, ...], ...],
'U': [[-0.9879925180204853, -0.9879925180204853, ...], ...],
'W': [3.635655546811611e-06, 8.318690786800488e-06, ...],
'Y': [-0.3777935957963988, -0.3112642082659244, ...]
}
The value associated with the Sampling key contains information about the generation
of the sample of model evaluations. In the case of quadrature-based projection method, it
can only have the 'Quadrature' value. The NSamples key contains the total number of
full-model evaluations that were run during the calculation.
The remaining keys contain, respectively:
• X: the quadrature nodes where the model is evaluated;
• U: the same points as X, but rescaled and transformed onto the domain of definition of
the orthogonal polynomials;
• Y: the full model evaluation at each of the quadrature notes;
• W: the quadrature weight of each point as in Eq. (1.26).
A posteriori error estimates
The dictionary myPCE_Quadrature['Error'] contains the normalized quadrature error
estimate from equation Eq. (1.28).
2.6.1.2 Advanced options
There are several advanced options for the calculation of PCE coefficients with the projec-
tion method, namely selecting the Smolyak’ sparse quadrature method in Section 1.5.1 and
specifying the quadrature level.
• Smolyak’ scheme: the Smolyak’ quadrature scheme can be enabled by adding the
following option:
MetaOpts['Quadrature'] = {'Type': 'Smolyak'}
Note that up to dimension M = 4 Smolyak’ scheme requires more nodes than full
quadrature for the same level of accuracy.
• Quadrature level: the quadrature level (by default set to l = p + 1, where p is the
maximum polynomial degree) can be set to the desired value (e.g. l = 15) as:
UQ[PY]LAB-V1.0-104 - 28 -
Polynomial chaos expansions
MetaOpts['Quadrature']['Level'] = 15
For a comprehensive list of the options available for the quadrature method, see Table 5.
2.6.2 Ordinary Least-Squares (OLS)
The calculation of PCE coefficients with Ordinary Least-Squares on a sample of N = 1, 000
model evaluations can be enabled with the following configuration:
# Reporting the previous configuration options as a reminder
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
# Specification of 14th degree, OLS-based PCE
MetaOpts['Degree'] = 14
MetaOpts['Method'] = 'OLS'
# Specification of the experimental design
MetaOpts['ExpDesign'] = {'NSamples': 1000}
# Creation of the metamodel:
myPCE_OLS = uq.createModel(MetaOpts);
Note that UQ[PY]LAB will create the experimental design and evaluate the model response
on it. Once the PCE is calculated, a report with basic information about the PCE results can
be printed on screen by:
uq.print(myPCE_OLS)
which produces the following report:
%------------ Polynomial chaos output ------------%
Number of input variables: 3
Maximal degree: 14
q-norm: 1.00
Size of full basis: 680
Size of sparse basis: 680
Full model evaluations: 1000
Leave-one-out error: 2.3251511e-08
Modified leave-one-out error: 2.5484956e-04
Mean value: 3.5000
Standard deviation: 3.7208
Coef. of variation: 106.310%
%--------------------------------------------------%
A visual representation of the spectrum of the resulting PCE coefficients can be created as
follows:
UQ[PY]LAB-V1.0-104 - 29 -
UQ[PY]LAB user manual
uq.display(myPCE_OLS)
which produces the image in Figure 7.
0 200 400 600
−8
−6
−4
−2
0
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 680
𝛼
log
10
(|𝑦
𝛼
|)
Figure 7: Graphical representation of the logarithmic spectrum of the PCE coefficients com-
puted by OLS. Again, most of the coefficients of the 680 basis elements are close to 0.
2.6.2.1 Accessing the results
Coefficients and basis
The coefficients and basis can be accessed from the dictionary myPCE_OLS['PCE']:
{
'Basis': {'PolyTypes': ['Legendre', 'Legendre', 'Legendre'], ...},
'Coefficients': [[3.499985658657572, 4.3041446260903194e-05, ...],
'Moments': {'Mean': 3.499998374934148, 'Var': 13.844721550794361},
}
Model evaluations
The model evaluations used to calculate the PCE coefficients with OLS can be accessed from
the myPCE_OLS['ExpDesign']dictionary :
{
'NSamples': 1000,
'Sampling': 'LHS',
'X': [[-1.5814659610815798, 1.0260808091053413, ...], ...],
'U': [[-0.5033962500754168, 0.3266116655616931, ...], ...],
'Y': [3.0741638768935946, 0.07656583917580881, ...]
}
The Sampling key has the 'LHS' value. It represents the sampling strategy used to create
the experimental design, see Section 2.5 for advanced options for the creation of the experi-
mental design.
A posteriori error estimates
The a posteriori error estimates are stored in the myPCE_OLS['Error']dictionary :
UQ[PY]LAB-V1.0-104 - 30 -
Polynomial chaos expansions
{
'LOO': 2.325151052510105e-08,
'ModifiedLOO': 0.0002548495560192674,
'normEmpError': 1.1070587842866738e-12,
}
where normEmpError and LOO corresponds to the empirical error estimate in Eq. (1.18)
and to the modified leave-one-out error in Eq. (1.20), respectively.
Note that the normEmpError is much smaller than LOO, as it does not account for over-
fitting.
For a comprehensive overview of the outputs available for the OLS method, see Table 23.
2.6.2.2 Advanced options
There are no OLS specific options per se. However, when combined with degree adaptive PCE
described in Section 1.3.4 and 2.7.1, two parameters can be set to tune the convergence of
the algorithm:
• Target accuracy: by default set to 0, it corresponds to ϵ
T
in Section 1.3.4. It can be
manually set to any value (e.g. 10
−8
) as follows:
MetaOpts['OLS'] = {'TargetAccuracy': 1e-8}
• Disable the modified LOO estimator: by default, ϵ
LOO
is calculated with the modified
estimator in Eq. (1.22) and the correction factor in (1.23). It is possible, however, to
enable the classical LOO estimator in Eq. (1.19) as follows:
MetaOpts['OLS']['ModifiedLOO'] = 0
Additional configuration options are available for the creation of the experimental design
from which the least-square regression is performed. They are listed in Section 2.5. A de-
tailed list of the available configuration options for OLS can be found in Table 6.
2.6.3 Sparse regression methods (LARS, OMP, SP, BCS)
Configuring a sparse PCE with a sparse regression solver is very similar to setting up a PCE
with OLS. All four sparse regression solvers available in UQ[PY]LAB (LARS, OMP, SP, BCS)
are used in basically the same way.
The code needed to create a basic LARS-based PCE with 1,000 sample points in the experi-
mental design is as follows:
# Reporting the previous configuration options as a reminder
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
# Specification of 14th degree LARS-based PCE
MetaOpts['Degree'] = 14
UQ[PY]LAB-V1.0-104 - 31 -
UQ[PY]LAB user manual
MetaOpts['Method'] = 'LARS'
# Specification of the experimental design
MetaOpts['ExpDesign'] = {'NSamples': 1000}
# Creation of the metamodel:
myPCE_sparse = uq.createModel(MetaOpts);
Instead of LARS, we could use OMP, SP, or BCS by assigning 'OMP', 'SP', or 'BCS' to the
key Method.
A report with basic information about the PCE results can be printed on screen by:
uq.print(myPCE_sparse)
which produces the following report:
%------------ Polynomial chaos output ------------%
Number of input variables: 3
Maximal degree: 14
q-norm: 1.00
Size of full basis: 680
Size of sparse basis: 39
Full model evaluations: 1000
Leave-one-out error: 8.6063410e-12
Modified leave-one-out error: 9.3348271e-12
Mean value: 3.5000
Standard deviation: 3.7208
Coef. of variation: 106.309%
%--------------------------------------------------%
A visual representation of the spectrum of the non-zero coefficients can visualized graphically
as follows:
uq.display(myPCE_sparse)
which produces the image in Figure 8a. Notice how the LARS solution only produces 39
non-zero coefficients, which is much sparser than the OLS solution with its 680 non-zero
coefficients (Figure 7). Similar observations can be made for the other three sparse solvers
OMP (Figure 8b), SP (Figure 8c), and BCS (Figure 8d).
2.6.3.1 Accessing the results
Coefficients and basis
The coefficients and basis can be accessed from the structure myPCE_sparse['PCE']:
{
'Basis': {'PolyTypes': ['Legendre', 'Legendre', 'Legendre'], ...},
'Coefficients': [3.4999999974846885, 0, 0, 1.6254179760353538, ...],
'Moments': {'Mean': 3.4999999974846885, 'Var': 13.844586608982498}
}
UQ[PY]LAB-V1.0-104 - 32 -
Polynomial chaos expansions
0 200 400 600
−6
−5
−4
−3
−2
−1
0
1
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 39
𝛼
log
10
(|𝑦
𝛼
|)
(a) LARS
0 200 400 600
−6
−5
−4
−3
−2
−1
0
1
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 68
𝛼
log
10
(|𝑦
𝛼
|)
(b) OMP
0 200 400 600
−6
−5
−4
−3
−2
−1
0
1
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 102
𝛼
log
10
(|𝑦
𝛼
|)
(c) SP
0 100 200 300
−3
−2
−1
0
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 20
𝛼
log
10
(|𝑦
𝛼
|)
(d) BCS
Figure 8: Graphical representation of the logarithmic spectra of the PCE coefficients com-
puted by the sparse regression methods LARS, OMP, SP, and BCS. The sparsity of these
solutions compared to the OLS solution in Figure 7 is clear.
Model evaluations
The experimental design structure containing the model evaluations used to calculate the
PCE coefficients is identical to that in OLS in Section 2.6.1.1.
{
'NSamples': 1000,
'Sampling': 'LHS',
'ED_Input': 'ED_Input',
'X': [[0.19377715485433855, 1.270031433237528,, ...], ...],
'U': [[0.06168118410670331, 0.404263560963674, ...], ...],
'Y': [7.304797547964152, 16.550041072666275, ...]
}
The Sampling field has the 'LHS' value. It represents the sampling strategy used to create
the experimental design. See Section 2.5 for advanced options for the creation of the experi-
mental design.
A posteriori error estimates
The a posteriori error estimates are stored in the myPCE_sparse['Error'] dictionary in
UQ[PY]LAB-V1.0-104 - 33 -
UQ[PY]LAB user manual
the same format as in Section 2.6.2.1:
{
'LOO': 8.606340987907494e-12,
'ModifiedLOO': 9.334827148115448e-12,
'normEmpError': 7.731314785646453e-12
}
Note that the LOO error for sparse PCE (≈ 10
−11
) is significantly smaller than for OLS (≈
10
−8
), even though the experimental design has the same size N = 1, 000.
Note that BCS relies on the k-fold cross-validation error instead of LOO. Therefore, in case of
BCS the error reported in myPCE_sparse.Error['LOO'](or
myPCE_sparse['Error.ModifiedLoo']) is not the LOO error, but the (modified) k-fold CV
error.
For a comprehensive overview of further outputs available for the LARS method, see Table 24;
for the OMP method, see Table 25.
2.6.3.2 Shared advanced options
Albeit the default settings are optimal for most real case scenarios, the sparse solvers allow
for the customization of several parameters that can be used to fine-tune the coefficients
estimation.
For all three sparse regression method relying on the LOO error (i.e., LARS, OMP, and SP),
by default the modified LOO error estimate of Eq. (1.22) is used. For BCS, which relies on
the k-fold CV error, the modification factor of Eq. (1.23) can be applied to the CV error, but
by default it is not. This can be changed as follows:
• Enable or disable the modified LOO (resp. CV) estimator: by default, ϵ
LOO
is calcu-
lated with the modified estimator in Eq. (1.22), using the correction factor in (1.23). It
is possible, however, to enable the classical LOO (resp. CV) estimator in Eq. (1.19) as
follows:
MetaOpts['LARS'] = {'ModifiedLoo': 0}
For OMP or SP, use the key OMP or SP instead of LARS. For BCS, use the field BCS,
and note that while the option is called ModifiedLoo, it actually applies to the k-fold
cross-validation error. Note that the classical estimator tends to be less sensitive to
over-fitting, hence generally producing denser and less accurate PCE models.
Additional configuration options are available for the creation of the experimental design for
which the sparse regression is performed. They are listed in Section 2.5.
2.6.3.3 Advanced options for LARS
Advanced options for LARS can be specified with the MetaOpts['LARS'] structure as fol-
lows:
UQ[PY]LAB-V1.0-104 - 34 -
Polynomial chaos expansions
• Disable the early stop criterion: with some models, LARS can stop prematurely and
yield inaccurate results. To disable the early stop criterion in Section 1.5.3.3, add:
MetaOpts['LARS'] = {'LarsEarlyStop': False}
Note: Disabling this option can significantly increase the computational time nec-
essary to calculate the coefficients. However, for small experimental de-
signs (i.e. N = 50 points), it is disabled by default to improve accuracy.
• Disable applying OLS for calculating LOO error: by default, ϵ
LOO
used for model
selection during LARS iteration is calculated as in Eq. (1.20) (or Eq. (1.22)). As men-
tioned in Section 1.4.2, the use of Eq. (1.20) requires the PCE model M
P C
(x) built
with OLS. As a result, at the end of each LARS iteration step, a least-square minimiza-
tion is applied to the selected basis functions. This OLS operation can be disabled by
setting
MetaOpts['LARS'] = {'HybridLoo': 0}
Nevertheless, Eq. (1.20) is still used for model selection. In this case, M
P C
(x) in
Eq. (1.20) corresponds to the PCE model whose coefficients are calculated by LARS
update at each iteration.
• Store the LARS iterations in memory: by default UQ[PY]LAB will not cache all
the LARS iterations after the algorithm ends, because it may require significant mem-
ory resources. This behaviour can be changed as follows: If this option is active,
a large array containing all of the coefficients for each LARS iteration is saved in:
myPCE_sparse['Internal']['PCE']['LARS']['coeff_array'].
A detailed list of the available configuration options for LARS can be found in Table 7.
2.6.3.4 Advanced options for OMP
Advanced options for OMP can be specified with the MetaOpts.OMP structure as follows:
• Disable the early stop criterion: with some models, OMP can stop prematurely and
yield inaccurate results. To disable the early stop criterion in Section 1.5.4.2, add:
MetaOpts['OMP'] = {'OmpEarlyStop': False}
Note: Disabling this option can significantly increase the computational time nec-
essary to calculate the coefficients. However, for small experimental de-
signs (i.e. N = 50 points), it is disabled by default to improve accuracy.
• Store the OMP iterations in memory: by default UQ[PY]LAB will not cache all the
OMP iterations after the algorithm ends, because it may require significant memory
resources. This behaviour can be changed as follows:
UQ[PY]LAB-V1.0-104 - 35 -
UQ[PY]LAB user manual
MetaOpts['OMP'] = {'KeepIterations': 1}
If this option is active, a large array containing all of the coefficients for each OMP
iteration is saved in: myPCE_OMP['Internal']['PCE']['OMP']['coeff_array'].
A detailed list of the available configuration options for OMP can be found in Table 8.
2.6.3.5 Advanced options for SP
SP has one hyperparameter, the number K of nonzero coefficients. It can be specified explic-
itly by setting
MetaOpts['SP'] = {'NNZ': K}
If K is not specified, it is determined by cross-validation from a range of 10 evenly spaced
values between 1 and min{
N
2
,
P
2
} (K = 1 excluded). By default, this is done by LOO as
described in 1.5.5.1. The user can also request k-fold cross-validation (default are k = 5
folds) by
MetaOpts['SP'] = {'CVMethod': 'kfold'}
The number of folds can be set in the optional key MetaOpts['SP']['NumFolds'].
A detailed list of the available configuration options for SP can be found in Table 9.
2.6.3.6 Advanced options for BCS
The hyperparameter σ
2
of BCS is determined by k-fold cross-validation (default: k = 10)
from the following range of values:
sigma2_range = N
*
np.var(Y)
*
np.power(10,np.linspace(-16,-1,10)). The num-
ber of folds can be changed in the optional key MetaOpts['BCS']['NumFolds'].
2.7 Basis adaptivity
All regression methods (OLS, LARS, OMP, SP, and BCS) provide a cross-validation error
estimator – the LOO error ϵ
LOO
in case of OLS, LARS, OMP, and SP, and the k-fold cross-
validation error ϵ
CV
in case of BCS – which can be used to develop basis-adaptive PCE as
described in Section 1.3.4.
In the following, everything described in terms of ϵ
LOO
holds equivalently for ϵ
CV
in case of
BCS.
2.7.1 Degree-Adaptive PCE
Degree-adaptive PCE is automatically enabled in UQ[PY]LAB if the MetaOpts['Degree']
option is an array of values instead of a single one. The following code creates a degree-
adaptive sparse PCE (LARS) with p ∈ [1, 30] from an experimental design with N = 256 for
the Ishigami function:
UQ[PY]LAB-V1.0-104 - 36 -
Polynomial chaos expansions
# Reporting the previous configuration options as a reminder:
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
# Specification of degree-adaptive LARS
MetaOpts['Degree'] = np.arange(1,31).tolist() # range of degrees to be tested
MetaOpts['Method'] = 'LARS'
# Specification of the experimental design
MetaOpts['ExpDesign'] = {
'Sampling': 'Sobol',
'NSamples': 256
}
# Creation of the metamodel:
myPCE_LARSAdaptive = uq.createModel(MetaOpts)
Despite the smaller experimental design, the degree-adaptive PCE converges to a maximal
PCE degree p = 18 with the lowest ϵ
LOO
amongst the examples presented in this section: The
a posteriori error estimates are stored in the myPCE_LARSAdaptive['Error'] dictionary
{
'LOO': 1.267083028494045e-17,
'ModifiedLOO': 2.863587078712277e-17,
'normEmpError': 3.981500516884674e-18
}
The resulting coefficients spectrum is shown in Figure 9.
0 500 1000
−12
−10
−8
−6
−4
−2
0
Mean
p=1
p=2
p=3
p>3
NNZ Coeffs: 73
𝛼
log
10
(|𝑦
𝛼
|)
Figure 9: Graphical representation of the logarithmic spectrum of the PCE coefficients for
degree-adaptive PCE. The analysis converged to p = 18 with N = 256.
It is therefore recommended to always specify in the PCE options a range of polynomial
degrees when using least-square methods, so as to allow adaptive PCE to adaptively choose
the best polynomial degree given the experimental design specifications.
UQ[PY]LAB-V1.0-104 - 37 -
UQ[PY]LAB user manual
2.7.1.1 Accessing the results
The outputs of a degree-adaptive PCE analysis are unchanged from their non-adaptive coun-
terparts, because only the iteration with the best ϵ
LOO
is stored. The best degree is stored in
myPCE['PCE']['Basis']['Degree'].
2.7.1.2 Advanced options
The default behaviour of the degree-adaptive scheme is to automatically stop increasing the
maximum degree if the ϵ
LOO
has not decreased for at least two subsequent iterations of the
algorithm. Experience shows that once over-fitting is detected with ϵ
LOO
on an experimental
design, further increasing the size of the polynomial basis results in worse PCE models.
In some rare cases, however, the algorithm can stop prematurely due to a local minimum in
the ϵ
LOO
vs. p curve. This can be prevented by setting the MetaOpts['DegreeEarlyStop']
flag to False :
MetaOpts['DegreeEarlyStop'] = False
When disabled, all the degrees specified in the MetaOpts['Degree'] list will be calculated,
and the best candidate will be chosen only at the end.
Note: Disabling this option can significantly increase the computational costs of the
PCE coefficients calculation, as the size of the polynomial basis (and hence the
number of coefficients that need to be calculated) increases very rapidly with the
maximum polynomial degree.
2.7.2 q-norm-Adaptive PCE
Due to the sparsity-of-effects principle (see Section 1.5.3) it is advantageous to use a trun-
cation scheme when searching for an optimal basis. To this end, UQ[PY]LAB offers a q-
norm-adaptive PCE set-up (for the q-norm truncation scheme see Section 1.3.3.2). The q-
norm adaptivity also makes use of the ϵ
LOO
error estimator and is automatically enabled in
UQ[PY]LAB if the option is an list of values instead of a single one. The q-norm adaptiv-
ity can be combined with the degree adaptivity. The following code creates a q-norm- and
degree-adaptive sparse PCE (LARS) with q-norm ∈ [0.5, 0.6, . . . , 1] and p ∈ [1, 30] from an
experimental design with N = 256 for the Ishigami function:
# Reporting the previous configuration options as a reminder:
MetaOpts= {
'Type': 'Metamodel',
'MetaType': 'PCE',
}
# Specification of degree- and q-norm-adaptive LARS
MetaOpts['TruncOptions'] = {
'qNorm': np.arange(0.5,1.01,0.1).tolist() # q-norms to be tested
}
MetaOpts['Degree'] = np.arange(1,31).tolist()
MetaOpts['Method'] = 'LARS'
UQ[PY]LAB-V1.0-104 - 38 -
Polynomial chaos expansions
# Specification of the experimental design
MetaOpts['ExpDesign'] = {
'Sampling': 'Sobol',
'NSamples': 256
}
# Creation of the metamodel:
myPCE_LARSAdaptive = uq.createModel(MetaOpts)
Note: When setting up a low-dimensional q-norm-adaptive PCE, small q-norm incre-
ments will not affect the basis.
2.7.2.1 Accessing the results
As in the case of degree-adaptive PCE, the outputs of a q-norm-adaptive PCE analysis are
unchanged from their non-adaptive counterparts, because only the iteration with the best
ϵ
LOO
is stored. The best q-norm is stored in myPCE['PCE']['Basis']['qNorm'].
2.7.2.2 Advanced options
The q-norm-adaptive scheme starts for every polynomial degree with the smallest specified
q-norm. The stopping criterion is a little more complex than the one for the degree-adaptive
scheme. For small polynomial degrees a small increase in the q-norm might not allow for
additional basis functions and will consequently not affect the set up PCE or its ϵ
LOO
. For
this reason, UQ[PY]LAB only regards iterations as interesting if they affect the ϵ
LOO
or the
basis size of the PCE. The scheme automatically stops increasing the q-norm once ϵ
LOO
did
not decrease in two subsequent interesting iterations. It also stops if the ϵ
LOO
stays the same
but the basis size increases twice.
The early stopping mechanism can be disabled by setting the MetaOpts['qNormEarlyStop']
flag to False :
MetaOpts['qNormEarlyStop'] = False
When disabled, for each degree, all q-norms specified in the
MetaOpts['TruncOptions']['qNorm'] list will be calculated, and the best candidate will
be chosen only at the end.
Note: As in the case of the degree-adaptive scheme, disabling this option can dramati-
cally increase the computational costs of the PCE coefficients calculation. Espe-
cially for high degrees and q-norm arrays with many entries, the computational
cost increases rapidly. However, when the number of coefficients is kept small (e.g.
less than 50 points) the early stop mechanism should be turned off to ensure the best
fit can be found. In any case, when the experimental design is small, there are
not many calculations to be done to construct the metamodel.
UQ[PY]LAB-V1.0-104 - 39 -
UQ[PY]LAB user manual
2.8 Use of a validation set
If a validation set is provided (see Table 12 in Section 3.1.10.), UQ[PY]LAB automatically
computes the validation error given in Eq. (1.17). To provide a validation set, the following
command shall be used:
MetaOpts['ValidationSet'] = {
'X': XVal,
'Y': Yval
}
The value of the validation error is stored in myPCE['Error']['Val'] next to the other
error measures (see Table 19 in Section 3.2.3) and will also be displayed when typing
uq.print(myPCE).
2.9 Manually specify inputs and computational models
The UQ[PY]LAB framework allows one to create many INPUT and MODEL objects in the same
session (see, e.g. UQ[PY]LAB User Manual – the MODEL module and UQ[PY]LAB User Manual
– the INPUT module). The default behaviour of the PCE module is to use as probabilistic input
(resp. computational model) the last created INPUT (resp. MODEL) object. This behaviour
can be altered by manually specifying the desired objects in the configuration as follows:
• Specify an INPUT object: an INPUT object myInput can be specified with:
MetaOpts['Input'] = myInput
• Specify a MODEL object: a MODEL object myModel can be specified with:
MetaOpts['FullModel'] = myModel
2.10 PCE of vector-valued models
The examples presented so far in this chapter dealt with scalar-valued models. In case the
model (or the experimental design, if manually specified) has multi-component outputs (de-
noted by N
out
), UQ[PY]LAB performs an independent PCE for each output component on the
shared experimental design. No additional configuration is needed to enable this behaviour.
A PCE with multi-component outputs can be found in the UQ[PY]LAB example in:
4-Multiple outputs
2.10.1 Accessing the results
Running a PCE calculation on a multi-component output model will result in a list of dic-
tionaries. As an example, a model with 9 outputs will produce the following output list of
9 dictionaries, each with the following keys:
UQ[PY]LAB-V1.0-104 - 40 -
Polynomial chaos expansions
{
'Basis',
'Coefficients',
'Moments'
}
Each of elements of the PCE dictionary is functionally identical to its scalar counterpart in
Section 2.6.1.1, 2.6.2.1 and 2.6.3.1.
Similarly, the myPCE['Error'] dictionary becomes a list of 9 dictionaries, each with the
following keys:
{
'LOO',
'ModifiedLOO',
'normEmpError'
}
2.11 Using a PCE as a predictor
Regardless on the strategy used to calculate the PCE coefficients or truncating the basis, the
PCE model in Eq. (1.3) can be used to predict new points outside of the experimental design.
Indeed, after a PCE MODEL object is created in UQ[PY]LAB it can be used just like an ordinary
model (for details, see the UQ[PY]LAB User Manual – the MODEL module).
2.11.1 Evaluate a PCE
Consider the Ishigami example in Section 2.1. After calculating the coefficients with any
of the methods described in Section 2.6, one can evaluate the PCE metamodel on point
x = {0.3, 1.0, 2.2} as follows:
X = [0.3, 1.0, 2.2]
# Evaluate the metamodel on the same input vector
YPC = uq.evalModel(myPCE_Quadrature, X)
print(YPC)
producing the following output:
[[5.9443194]]
which can be compared to the true model:
As most functions within UQ[PY]LAB, model evaluations are vectorized, i.e. evaluating mul-
tiple points at a time is much faster than repeatedly evaluating one point at a time. To eval-
uate the response of the PCE metamodel on an input sample of size N = 10
5
in UQ[PY]LAB,
one can write (for details on how to use the input module to sample distributions, see the
UQ[PY]LAB User Manual – the INPUT module):
X = uq.getSample(N=1e5)
Y = uq.evalModel(myModel, X)
YPC = uq.evalModel(myPCE_Quadrature, X)
UQ[PY]LAB-V1.0-104 - 41 -
UQ[PY]LAB user manual
The histogram and scatter plots of the Y and YPC vectors are given in Figure 10. Due to the
high accuracy of the model, the original function and the metamodel are virtually indistin-
guishable.
−10 −5 0 5 10 15
0
1000
2000
3000
4000
5000
Y
Y
PC
Y
Counts
−10 −5 0 5 10 15
−10
−5
0
5
10
15
Y
Y
PC
Figure 10: Histogram and scatter plots of true vs. metamodelled responses of the Ishigami
function to a sample of the input of size n = 10
5
.
2.11.2 Local error estimates
Getting local error estimates from a PCE consists of two parts: setting up a bootstrap PCE
(bPCE; see Section 1.7.1) and using the bPCE to compute the local error estimates of the PCE
point-wise predictions.
2.11.2.1 Set up bPCE
To set up a bPCE, the number of bootstrap replications B needs to be specified before creating
the PCE metamodel as follows:
MetaOpts['Bootstrap'] = {'Replications': 100}
The PCE produced by uq.createModelwill be similar to the one set up without the bootstrap
specification. The additional information is stored in myPCE['Internal']['PCE']['Bootstrap']
dictionary and used to compute the error estimates when evaluating the metamodel.
2.11.2.2 Apply bPCE to get error estimates
After its creation, the bPCE can be applied on a sample using uq.evalModelto get predic-
tions, just like a non-bootstrap PCE (see Section 2.11.1). The extended capabilities can be
employed by calling more output arguments:
YPC, YPC_var, YPC_replications = uq.evalModel(myPCE, X, nargout=3)
Here, YPC_var provides the sample variance of the predictions at each point. YPC_replications
is a N × B × N
out
list of lists that contains the bPCE predictions for each of the B bootstrap
replications. Using the bPCE predictions and the Python function np.quantile empirical
quantiles can be employed to get error margins:
UQ[PY]LAB-V1.0-104 - 42 -
Polynomial chaos expansions
# 95% confidence interval for the prediction on the fifth point:
loBound, upBound = np.quantile(YYPC_BootSample(5,:),[0.025 0.975]);
2.12 Manually specifying PCE parameters (predictor-only mode)
It is also possible to use the PCE module in UQ[PY]LAB to build custom PCE-based models
that can be used as MODEL objects as in Section 2.11. This allows, e.g., to import a meta-
model calculated with other software within the UQ[PY]LAB framework, or even or to create
one ad-hoc. In the following, we exemplify how to create a custom PCE with the following
characteristics:
• standard normal input variables;
• up to second degree Hermite polynomials polynomial basis;
• only three non-zero coefficients: y
[0,0]
= 5, y
[0,1]
= 1, y
[1,1]
= 3.
from uqpylab import sessions
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()
## Set the random seed for reproducibility
uq.rng(100, 'twister');
# Create an Input object
inputOpts = {'Marginals': uq.StdNormalMarginals(2)}
myInput = uq.createInput(inputOpts)
# Create a custom PCE
MetaOpts = {
'Input': myInput['Name'],
'Type': 'Metamodel',
'MetaType': 'PCE',
'Method': 'Custom'
}
# Basis: polynomial families
MetaOpts['PCE'] = {'Basis': {'PolyTypes': ['Hermite', 'Hermite']}}
# Basis: polynomial alpha indices
MetaOpts['PCE']['Basis']['Indices'] = [[0, 0], [0, 1], [1, 1]]
# PCE coefficients (same order as MetaOpts['PCE']['Basis']['Indices'])
MetaOpts['PCE']['Coefficients'] = [[5], [1], [3]]
# Create the metamodel
myPCE = uq.createModel(MetaOpts)
# Evaluate the model on a sample of the input
X = uq.getSample(N=1000)
Y = uq.evalModel(X=X)
UQ[PY]LAB-V1.0-104 - 43 -
UQ[PY]LAB user manual
Note that UQ[PY]LAB takes care automatically of any isoprobabilistic transformation be-
tween the probabilistic input model and the space onto which the specified polynomial fam-
ilies are orthogonal.
When the desired metamodel has more than one output, it is sufficient to specify the same
information for each of the outputs by adding an index i to the MetaOpts['PCE'][i])
dictionary.
2.13 Using a PCE with constant input variables
In some analyses, one may need to assign a constant value to one or to a set of input variables.
When this is the case, the PCE metamodel is built by internally removing the constant variable
from the inputs. This process is transparent to the users as they shall still evaluate the model
using the full set of parameters (including those which were set constant). UQ[PY]LAB will
automatically and appropriately account for the set of input variables which were declared
constant.
To set a parameter to constant, the following command can be used (See UQ[PY]LAB User
Manual – the INPUT module):
inputOpts['Marginals'] = {
'Type': 'Constant',
'Parameters': value
}
Furthermore, when the standard deviation of a input is set to zero, UQ[PY]LAB automatically
sets this variable’s marginal to the type Constant. For example, the following uniformly
distributed variable whose upper and lower bounds are identical is automatically set to a
constant with value 1:
inputOpts['Marginals'] = {
'Type': 'Uniform',
'Parameters': [1, 1]
}
UQ[PY]LAB-V1.0-104 - 44 -
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 configuration options and output quantities. Due to
the complexity of the algorithms implemented, it is not uncommon to employ nested dictio-
naries to fine-tune the inputs and outputs. Throughout this reference guide, a table-based
description of the configuration 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 field is put here
which corresponds to the following syntax:
Input = {
'Name' : 'My Input'
}
The columns, from left to right, correspond to the name, the data type and a brief description
of each 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)
45
UQ[PY]LAB user manual
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 define 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
□ Val2Opt2 Float Description
UQ[PY]LAB-V1.0-104 - 46 -
Polynomial chaos expansions
3.1 Create a PCE metamodel
Syntax
myPCE = uq.createModel(MetaOpts)
Input
The dictionary variable MetaOpts contains the following keys:
Table 3: MetaOpts
Type 'Metamodel' Select the metamodelling tool
MetaType 'PCE' Select polynomial chaos expansion
□ Input INPUT object Probabilistic input model (See
Section 2.9)
□ Name String Unique identifier for the metamodel
□ Display String
default: 'standard'
Level of information displayed by the
methods.
'quiet' Minimum display level, displays
nothing or very few information.
'standard' Default display level, shows the most
important information.
'verbose' Maximum display level, shows all the
information on runtime, like updates
on iterations, etc.
□ Degree Integer scalar Maximum polynomial degree
List of integers Set of polynomial degrees for
degree-adaptive polynomial chaos
(Section 2.7.1)
□ DegreeEarlyStop Boolean
default: True
Toggle polynomial degree early stop
criterion on / off (Section 2.7.1).
□ PolyTypes 1 × M List of strings List of polynomial families to be used
to build the PCE basis. The default
choice is given in Table 2. If one of
the polynomial families is Jacobi or
Laguerre the corresponding
parameters should be set with
PolyTypesParams.
□ PolyTypesParams 1 × M List of floats Set of parameters to be used to build
the PCE basis. It is only used when
PolyTypes contains Jacobi or
Laguerre polynomials. See
Section 2.4 for usage example.
UQ[PY]LAB-V1.0-104 - 47 -
UQ[PY]LAB user manual
□ TruncOptions Table 4 Basis truncation (Section 1.3.3)
□ qNormEarlyStop Boolean
default: True
Toggle hyperbolic truncation norm
early stop criterion on / off
(Section 2.7.2).
□ Method String
default: 'LARS'
Coefficients calculation method
(Section 2.6)
'Quadrature' Quadrature (Section 1.5.1)
'OLS' Ordinary Least-Squares Regression
(Section 1.5.2.1)
'LARS' Least-Angle Regression
(Section 1.5.3)
'OMP' Orthogonal Matching Pursuit
(Section 1.5.4)
'SP' Subspace Pursuit (Section 1.5.5)
'BCS' Bayesian Compressive Sensing
(Section 1.5.6)
'Custom' User-defined PCE coefficients and
basis (no calculations)
⊞ Quadrature Table 5 Quadrature options (Section 2.6.1.2)
⊞ OLS Table 6 OLS-specific options
(Section 2.6.2.2)
⊞ LARS Table 7 LARS-specific options (Section 2.6.3)
⊞ OMP Table 8 OMP-specific options (Section 2.6.3)
⊞ PCE Table 15 Custom-PCE parameters
(Section 2.12). Use the same format
as the default output of the
calculation
□ FullModel MODEL object UQ[PY]LABmodel used to create an
experimental design (Section 2.9)
□ ExpDesign Table 11 Experimental design-specific options
(Section 2.5)
□ ValidationSet Table 12 Validation set components
(Section 2.8)
□ Bootstrap Table 13 Bootstrapping options
(Section 2.11.2)
3.1.1 Truncation options
The truncation strategies described in Section 1.3.3 can be specified with the TruncOptions
key as described in Section 2.4.2. The full list of available options is reported in Table 4.
Table 4: MetaOpts['TruncOptions']
UQ[PY]LAB-V1.0-104 - 48 -
Polynomial chaos expansions
□ qNorm Float
default: 1
Hyperbolic truncation scheme
(Section 1.3.3). Corresponds to
0 < q ≤ 1 in Eq. (1.14)
List of Floats Set of hyperbolic norms to be tested
for q-norm-adaptive polynomial
chaos (Section 2.7.1).
□ MaxInteraction Integer
default: M
Maximum rank truncation: limit
basis terms to MaxInteraction
variables (Section 1.3.3)
□ Custom P × M List of Integers Manual specification of the A index
set in Eq. (1.3)
3.1.2 PCE Coefficients calculation options
Method-specific options for the calculation of the PCE coefficients are reported in Tables 5 to 9.
3.1.3 Quadrature-specific options
Table 5: MetaOpts['Quadrature']
□ Level Integer
default: p + 1
Quadrature level
□ Type String
default: dependent on M
Quadrature type (full or sparse)
'Full' if M < 4 Full tensor-product quadrature
'Smolyak' if M ≥ 4 Smolyak’ sparse quadrature
□ Rule String
default: 'Gaussian'
Quadrature rule
'Gaussian' Gaussian quadrature
3.1.4 OLS-specific options
Table 6: MetaOpts['OLS']
□ TargetAccuracy Float
default: 0
Early stop leave-one-out error
threshold for degree-adaptive PCE
3.1.5 LARS-specific options
Table 7: MetaOpts['LARS']
□ LarsEarlyStop Boolean
default: dependent on N
Enable early stop during the LARS
adaptive basis selection
(Section 1.5.3.3).
False if N < 50
UQ[PY]LAB-V1.0-104 - 49 -
UQ[PY]LAB user manual
True if N ≥ 50
□ TargetAccuracy Float
default: 0
Early stop leave-one-out error
threshold.
□ KeepIterations Boolean
default: False
Store additional information about
LARS iterations.
Warning: memory consuming.
□ HybridLoo Boolean
default: True
Enable/Disable applying OLS at each
LARS iteration for calculating the
LOO error
□ ModifiedLoo Boolean
default: 1
Enable/Disable using the modified
Leave-One-Out error for model
selection (see Section 1.4.3)
3.1.6 OMP-specific options
Table 8: MetaOpts['OMP']
□ OmpEarlyStop Boolean
default: dependent on N
Enable early stop during the OMP
adaptive basis selection
(Section 1.5.4.2).
False if N < 50
True if N ≥ 50
□ TargetAccuracy Float
default: 0
Early stop leave-one-out error
threshold.
□ KeepIterations Boolean
default: False
Store additional information about
OMP iterations.
Warning: memory consuming.
□ ModifiedLoo Boolean
default: 1
Enable/Disable using the modified
Leave-One-Out error for model
selection (see Section 1.4.3)
3.1.7 SP-specific options
Table 9: MetaOpts['SP']
□ NNZ Integer
default: []
User-specified value for the
hyperparameter K (desired number
of nonzero coefficients)
□ CVMethod String
default: 'loo'
Cross-validation method for
determining the hyperparameter K.
Only used if NNZ is not specified or
empty.
'loo' Leave-one-out cross-validation
'kfold' k-fold cross-validation
UQ[PY]LAB-V1.0-104 - 50 -
Polynomial chaos expansions
□ NumFolds Integer
default: 5
Number of folds for k-fold
cross-validation. Only used if
CVMethod = 'kfold'.
□ ModifiedLoo Boolean
default: 1
Enable/Disable using the modified
Leave-One-Out error for model
selection (see Section 1.4.3)
3.1.8 BCS-specific options
Table 10: MetaOpts['BCS']
□ NumFolds Integer
default: 10
Number of folds for k-fold
cross-validation.
□ ModifiedLoo Boolean
default: 0
Enable/Disable using the modified
cross-validation error for model
selection (see Section 1.4.3)
3.1.9 Experimental design
If a model is specified, UQ[PY]LAB can automatically create an experimental design for PCE.
The available options are listed in Table 11.
Table 11: MetaOpts['ExpDesign']
⊕ Sampling String
default: 'LHS'
Sampling type
'MC' Monte Carlo sampling
'LHS' Latin Hypercube sampling
'Sobol' Sobol sequence sampling
'Halton' Halton sequence sampling
□ NSamples Integer The number of samples to draw. It is
required when Sampling is
specified.
⊕ X N ×M Float User defined experimental design X.
If specified, Sampling is ignored.
⊕ Y N ×N
Out
Float User defined model response Y. If
specified, Sampling is ignored.
3.1.10 Validation Set
If a validation set is provided, UQ[PY]LAB automatically calculates the validation error of
the created PCE. The required information is listed in Table 12.
Table 12: MetaOpts['ValidationSet']
X N ×M Float User-specified validation set X
V al
UQ[PY]LAB-V1.0-104 - 51 -
Polynomial chaos expansions
3.2 Accessing the results
Syntax
myPCE = uq.createModel(MetaOpts)
Output
Regardless on the configuration options given at creation time in the MetaOpts dictionary,
all PCE metamodels share the same output dictionary, given in Table 14.
Table 14: myPCE
Name String Unique name of the PCE metamodel
Options Table 3 Copy of the MetaOpts dictionary used to
create the metamodel
PCE Table 15 Information about all the elements of Eq. (1.3)
ExpDesign Table 18 Experimental design used for calculating the
coefficients
Error Table 19 Error measures of the metamodelling
calculation results (Section 1.4)
Internal Table 20 Internal state of the MODEL object (useful for
debug/diagnostics)
3.2.1 Polynomial chaos expansion information
All the information needed to evaluate and post-process a PCE are contained in the myPCE['PCE']
dictionary. They include a basis and a set of coefficients (see Eq. (1.3)). Their format is given
in Tables 15.
Note that in case the model considered has a N
out
-dimensional output, each output variable
Y
i
is treated separately and stored in myPCE['PCE'][i].
Table 15: myPCE['PCE'][i]
Coefficients P × 1 Float Truncated PCE coefficients.
Moments Table 16 Post-processed moments of the PCE
(Section 1.6.1).
Basis Table 17 Information about the truncated polynomial
basis.
Table 16: myPCE['PCE'][i]['Moments']
Mean Float Mean of the PCE (Eq. (1.37))
Var Float Variance of the PCE (Eq. (1.38))
UQ[PY]LAB-V1.0-104 - 53 -
UQ[PY]LAB user manual
Table 17: myPCE['PCE'][i]['Basis']
Degree Float Maximum total polynomial degree of the basis
Indices P × M Float
(Sparse)
Truncated set of indices A in Eq. (1.3)
(Section 1.3)
PolyTypes 1 × M List of
strings
Polynomial family for each input variable. In
the current version of UQ[PY]LAB, each
element can be one of 'Legendre',
'Hermite', 'Laguerre' or 'Jacobi'
MaxCompDeg 1 × M Float Maximum degree in each input variable of
polynomials with non-zero coefficients
MaxInteractions Float Maximum rank of the polynomials with
non-zero coefficients
qNorm Float The q-norm of the basis (1.0 unless q-norm
adaptivity was used, see Section 2.7.2)
3.2.2 Experimental design information
The experimental design and the corresponding model responses onto which the PCE coeffi-
cients are calculated are stored in the myPCE['ExpDesign'] dictionary. They are accessible
as follows:
Table 18: myPCE['ExpDesign']
NSamples Float The number of samples
Sampling String The sampling method
ED_input INPUT object The input module that represents the reduced
polynomial input (X in Section 1.3.2.1)
X N × M Float The experimental design values
U N × M Float The experimental design values in the reduced
space
Y N × N
out
Float The output Y that corresponds to the input X
W N × 1 Float The Gaussian quadrature weights
corresponding to each quadrature node (only
available when the coefficients are calculated
with the 'Quadrature' method)
3.2.3 Error estimates
The error estimates described in Sections 1.4.1 – 1.4.3 and 2.8 are available in the myPCE['Error']
output key, as described in Table 19.
UQ[PY]LAB-V1.0-104 - 54 -
Polynomial chaos expansions
Table 19: myPCE['Error']
LOO Float Leave-One-Out error (see Section 1.4.2).
ModifiedLOO Float Modified Leave-One-Out error (see
Section 1.4.3).
normEmpError Float Normalized Empirical Error (see Section 1.4.1)
Val Float Validation error (see Eq. (1.17) and
Section 2.8). Only available if a validation set
is provided (see Table 12).
3.2.4 Internal keys (advanced)
Additional information that can be useful to the advanced user is stored in the myPCE['Internal']
key. Both runtime information and complex data structures used internally by the UQ[PY]LAB
PCE module are stored in this dictionary. The general dictionary of the myPCE['Internal']
key is reported in Table 20. Note that not all the keys are always available, as they depend
on the original configuration options.
Table 20: myPCE['Internal']
Input INPUT object The probabilistic input model used to build the
PCE
FullModel MODEL object Full computational model used to calculate the
model response (if available)
Error Table 21 Additional information about the PCE error
estimation given in myPCE['Error']
PCE Table 22 Additional information on the PCE calculation.
Runtime Table 27 Temporary variables and configuration flags
used during the calculation of the PCE
coefficients
Table 21: myPCE['Internal']['Error']
LOO_lars Float LOO error as calculated by LARS
LOO_omp Float LOO error as calculated by OMP
Table 22: myPCE['PCE']
Degree Float Final PCE degree
DegreeEarlyStop Boolean Polynomial degree early stop criterion
Method String Algorithm used to calculate the coefficients
OLS Table 23 OLS-specific information
UQ[PY]LAB-V1.0-104 - 55 -
Polynomial chaos expansions
Table 23: myPCE['Internal']['PCE']['OLS']
TargetAccuracy Float Degree-adaptive early stop threshold
LOO List LOO error for each of the tested degrees in
degree-adaptive mode.
normEmpError List Normalized empirical error for each of the
tested degrees in degree-adaptive mode.
Table 24: myPCE['Internal']['PCE']['LARS']
TargetAccuracy Float Degree-adaptive early stop threshold
LarsEarlyStop Boolean Early stop in LARS iterations flag
HybridLoo Boolean Enable/Disable applying OLS at each LARS
iteration for calculating the LOO error
ModifiedLoo Boolean Enable/Disable the “modified LOO” error
estimation in Eq. (1.20)
LOO List LOO error for each of the tested degrees in
degree-adaptive mode.
normEmpError List Normalized empirical error for each of the
tested degrees in degree-adaptive mode
KeepIterations Boolean Enable storage of LARS iterations (warning:
memory intensive)
coeff_array List of lists with
Float entries
Matrix of coefficients as calculated by each
iteration of LARS (requires
KeepIterations = 1)
a_scores List with Float
entries
Array of scores for each iteration of LARS
(score = 1-LOO). The final basis selected by
LARS is the one with the maximum a_score
loo_scores List with Float
entries
Array of LOO erro values for each iteration of
LARS
lars_idx List with
Integer entries
Array of indices representing the regressor
chosen at each LARS iteration
Table 25: myPCE['Internal']['PCE']['OMP']
TargetAccuracy Float Degree-adaptive early stop threshold
OmpEarlyStop Boolean Early stop in OMP iterations flag
ModifiedLoo Boolean Enable/Disable the “modified LOO” error
estimation in Eq. (1.20)
LOO List with Float
entries
LOO error for each of the tested degrees in
degree-adaptive mode.
normEmpError List with Float
entries
Normalized empirical error for each of the
tested degrees in degree-adaptive mode
UQ[PY]LAB-V1.0-104 - 57 -
UQ[PY]LAB user manual
KeepIterations Boolean Enable storage of OMP iterations (warning:
memory intensive)
coeff_array List of lists with
Float entries
Matrix of coefficients as calculated by each
iteration of OMP (requires
KeepIterations = 1)
a_scores List with Float
entries
Array of scores for each iteration of OMP
(score = 1-LOO). The final basis selected by
OMP is the one with the maximum a_score
loo_scores List with Float
entries
Array of LOO error values for each iteration of
OMP
omp_idx List with
Integer entries
Array of indices representing the regressor
chosen at each OMP iteration
Table 26: myPCE['Internal']['PCE']['Basis']
Truncation Dictionary Dictionary with the truncation options used to
generate the basis. See Table 4
Table 27: myPCE['Internal']['Runtime']
isInitialized Boolean A flag that determines whether the current
metamodel has been initialized
M Float The INPUT dimension
MnonConst Integer The number of non-constants in the input
nonConstIdx List with
Integer entries
The indices of the constant variables
isCalculated Boolean A flag that determines whether all the
necessary quantities of the metamodel have
been calculated
Nout Integer The Output dimension
current_output Integer The current output (This is used during the
calculation of the metamodel)
degree_index Integer Index of the PCE being considered in the
current status of the calculation
UQ[PY]LAB-V1.0-104 - 58 -
