Least Squares Polynomial Chaos Expansion: A Review of Sampling Strategies

TL;DR

This paper reviews various sampling strategies for least squares polynomial chaos expansion, introduces a hybrid method, and compares their performance across multiple numerical examples to guide practitioners.

Contribution

It presents a comprehensive review of sampling methods for PCE, introduces a novel hybrid sampling technique, and provides empirical comparisons to inform method selection.

Findings

Alphabetic-coherence-optimal sampling outperforms others in high-order and low oversampling scenarios.

Hybrid sampling method combines advantages of ODE and coherence-optimal techniques.

Performance varies significantly depending on problem complexity and sampling strategy.

Abstract

As non-institutive polynomial chaos expansion (PCE) techniques have gained growing popularity among researchers, we here provide a comprehensive review of major sampling strategies for the least squares based PCE. Traditional sampling methods, such as Monte Carlo, Latin hypercube, quasi-Monte Carlo, optimal design of experiments (ODE), Gaussian quadratures, as well as more recent techniques, such as coherence-optimal and randomized quadratures are discussed. We also propose a hybrid sampling method, dubbed alphabetic-coherence-optimal, that employs the so-called alphabetic optimality criteria used in the context of ODE in conjunction with coherence-optimal samples. A comparison between the empirical performance of the selected sampling methods applied to three numerical examples, including high-order PCE's, high-dimensional problems, and low oversampling ratios, is presented to provide a road map for practitioners seeking the most suitable sampling technique for a problem at hand. We observed that the alphabetic-coherence-optimal technique outperforms other sampling methods, specially when high-order ODE are employed and/or the oversampling ratio is low.

Tables

Table 1. Table 1: Correspondence of Wiener-Askey polynomial chaos and probability distribution of the random variables [ 127 ] .

${\Xi }_k$	Polynomial type	Support
Gaussian	Hermite	(-$\infty$,+$\infty$)
Gamma	Laguerre	(0,+$\infty$)
Beta	Jacobi	[a,b]
Uniform	Legendre	[a,b]

Table 2. Table 2: LIB model parameters [ 104 , 103 , 102 ] .

Parameter	Value
$R_{sp0}$	0.0272 $\Omega$
$R_{sp1}$	$1.087\times {10}^{-16}$ $\Omega$
$R_{sp2}$	34.64
$R_s$	0.0067 $\Omega$
$R_p$	10000 $\Omega$
$C_{b0}$	19.8 F
$C_{b1}$	1745 F
$C_{b2}$	-1.5 F
$C_{b3}$	-200.2 F
$C_s$	115.28 F
$C_{sp}$	316.69 F
$q_{max}$	31100 C
$C_{max}$	30807 C

Table 3. Table 3: Equations describing the electrical equivalent circuit model for LIB’s [ 104 , 103 , 102 ] .

	Governing equation
Voltages	$V_b=q_b/C_b,V_{sp}=q_{sp}/C_{sp},V_s=q_s/C_s,V_p=V_b-V_{sp}-V_s$
Currents	$I_p=V_p/R_p,I_b=I_p+I,I_{sp}=I_b-V_{sp}/R_{sp},I_s=I_b-V_s/R_s$
Charges	$\dot{q_b}=-I_b,{\dot{q}}_{sp}=I_{sp},\dot{q_s}=I_s$

Equations

$$u(\mathbf{\Xi })=\sum _{j=1}^{\infty }{c}_j{\psi }_j(\mathbf{\Xi }),$$

$$u(\mathbf{\Xi })=\sum _{j=1}^P{c}_j{\psi }_j(\mathbf{\Xi })+ϵ(\mathbf{\Xi }),$$

$$\mathbf{u}=\mathbf{\Psi }\mathbf{c}+\mathbf{ϵ},$$

$$\min_{\color[rgb]{0,0,0}{\bm{c}}}\|\bm{W}\bm{u}-\bm{W}\bm{\Psi}\bm{c}\|_{2},$$

$$(\mathbf{W}\mathbf{\Psi }{)}^T(\mathbf{W}\mathbf{\Psi })\hat{\mathbf{c}}=(\mathbf{W}\mathbf{\Psi }{)}^T\mathbf{W}\mathbf{u}.$$

$$\mathbf{M}=\frac{1}{N}{\mathbf{\Psi }}^T\mathbf{\Psi },$$

$$\mathbf{M}=\frac{1}{N}(\mathbf{W}\mathbf{\Psi }{)}^T(\mathbf{W}\mathbf{\Psi }),$$

$$\mathcal{R}\left(\mathbf{x},t,\mathbf{\Xi };u\right)=0,(\mathbf{x},t,\mathbf{\Xi })\in \mathcal{D}\times [0,T]\times \Omega .$$

$${\psi }_{\mathbf{j}}(\mathbf{\Xi })=\prod _{k=1}^d{\psi }_{j_k}({\Xi }_k),$$

$$P=\frac{(p+d)!}{p!d!}.$$

$$c_j=\mathbb{E}[u{\psi }_j].$$

$$(1-\delta )\le {\lambda }_{\min }(\mathbf{M})\le {\lambda }_{\max }(\mathbf{M})\le (1+\delta );$$

$$\mathcal{K}(\mathbf{M})=\frac{{\lambda }_{\max }(\mathbf{M})}{{\lambda }_{\min }(\mathbf{M})}\le \frac{1+\delta }{1-\delta };$$

$$\Vert \hat{u}{\Vert }_{L_2(\Omega ,f)}=\Vert \hat{\mathbf{c}}\Vert \le \frac{(1+\delta {)}^{1/2}}{1-\delta }\frac{1}{\sqrt{N}}\Vert \mathbf{Wu}\Vert ,$$

$$\frac{\Vert {\hat{\mathbf{c}}}_{\delta }-\hat{\mathbf{c}}\Vert }{\Vert \hat{\mathbf{c}}\Vert }\le \mathcal{K}(\mathbf{M})\frac{\Vert \mathbf{W}\delta \mathbf{u}\Vert }{\Vert \mathbf{W}\mathbf{u}\Vert }.$$

$$\mu(\bm{Y}):=\sup_{\color[rgb]{0,0,0}{\bm{\xi}\in\Omega}}\sum_{j=1}^{P}|w(\bm{\xi})\psi_{j}(\bm{\xi})|^{2}$$

$$\mathbb{P}\left(\Vert \mathbf{M}-\mathbf{I}\Vert \le \delta \right)\ge 1-2P\exp (-c_{\delta }N{\mu }^{-1}),$$

$$\mathbb{P}(\mathcal{E})\ge 1-2P\exp (-0.1N{\mu }^{-1}),$$

$$\mathbb{E}\left(\Vert u-\hat{u}{\Vert }_{L_2(\Omega ,f)}^2;\mathcal{E}\right)\le \mathbb{E}({ϵ}^2)\left(1+\frac{4\mu }{N}\right),$$

$$\mathbb{E}(X;\mathcal{E})={\int }_{\mathcal{E}}X(\mathbf{\xi })f(\mathbf{\xi })d\mathbf{\xi }=\mathbb{E}(X|\mathcal{E})\mathbb{P}(\mathcal{E})$$

$$B(\mathbf{\xi }):=\sqrt{\sum _{j=1}^P|{\psi }_j(\mathbf{\xi }){|}^2}.$$

$$f_{\mathbf{Y}}(\mathbf{\xi }):=c^{2}f(\mathbf{\xi })B^2(\mathbf{\xi })$$

$$f_Y({\xi }_k):=\frac{1}{\pi \sqrt{1-{\xi }_k^2}},$$

$$w(\mathbf{\xi }):=\prod _{k=1}^d(1-{\xi }_k^2{)}^{1/4}.$$

$$\frac{|{\mathbf{\Psi }}_{\text{new}}^T{\mathbf{\Psi }}_{\text{new}}|}{|{\mathbf{\Psi }}_{\text{old}}^T{\mathbf{\Psi }}_{\text{old}}|}=1+\Delta ({\mathbf{\psi }}_i,{\mathbf{\psi }}_j),$$

$$\Delta ({\mathbf{\psi }}_i,{\mathbf{\psi }}_j)=d({\mathbf{\psi }}_j)-\left[d({\mathbf{\psi }}_i)d({\mathbf{\psi }}_j)-(d({\mathbf{\psi }}_i,{\mathbf{\psi }}_j){)}^2\right]-d({\mathbf{\psi }}_i);$$

$$d({\mathbf{\psi }}_i)={\mathbf{\psi }}_i^T({\mathbf{\Psi }}_{\text{old}}^T{\mathbf{\Psi }}_{\text{old}}{)}^{-1}{\mathbf{\psi }}_i;$$

$$d({\mathbf{\psi }}_i,{\mathbf{\psi }}_j)={\mathbf{\psi }}_i^T({\mathbf{\Psi }}_{\text{old}}^T{\mathbf{\Psi }}_{\text{old}}{)}^{-1}{\mathbf{\psi }}_j.$$

$$D_N^{\ast }:=\underset{J\in {\mathcal{J}}^{\ast }}{\sup }\left|\frac{1}{N}\cdot \#\{i:{\mathbf{\xi }}^{(i)}\in J,1\le i\le N\}-{\lambda }_d(J)\right|,$$

$$J=\prod _{k=1}^d[0,v_k{)}^d,0\le v_k\le 1,$$

$$\text{relative error}=\frac{\Vert {\mathbf{u}}_v-{\mathbf{\Psi }}_v\hat{\mathbf{c}}{\Vert }_2}{\Vert {\mathbf{u}}_v{\Vert }_2},$$

$$u({\mathbf{\xi }}^{(i)})=\sum _{j=1}^P{c}_j{\psi }_j({\mathbf{\xi }}^{(i)}),i=1,\cdots ,N,$$

Full text

††journal: Elsevier

1 Introduction

Reliable and rigorous simulations of real world engineering problems involve characterization of the often uncertain system inputs, e.g., material properties, initial, or boundary conditions, and quantification of their impact on the output quantities of interest (QoI’s). This is the subject of uncertainty quantification (UQ), an emerging field in computational engineering and science, which aims at providing tools for assessing the credibility of model predictions and facilitating decision making under uncertainty.

A major class of UQ techniques are probabilistic where uncertain parameters are represented by random variables or processes. Among the probabilistic UQ approaches, stochastic spectral methods based on polynomial chaos (PC) expansions [122, 15, 40, 125] have received special attention due to their advantages over traditional UQ techniques; see, e.g., [40, 117, 86, 27, 42]. Let $\left(\Omega,\mathcal{T},\mathcal{P}\right)$ be a complete probability space where $\Omega$ is the sample set (or design space in the context of design of experiments) and $\mathcal{P}$ is a probability measure on the $\sigma-$field $\mathcal{T}$. Also assume that the system input uncertainty is described by the random vector $\bm{\Xi}=\left(\Xi_{1},\cdots,\Xi_{d}\right):\Omega\rightarrow\mathbb{R}^{d}$, $d\in\mathbb{N}$, consisting of $d$ independent identically distributed (i.i.d.) random variables $\Xi_{k}$ with realizations $\bm{\xi}=(\xi_{1},\cdots,\xi_{d})$. The marginal probability density function (pdf) of $\Xi_{k}$ and the joint pdf of $\bm{\Xi}$ are denoted by $f(\xi_{k})$ and $f(\bm{\xi})=\prod_{k=1}^{d}f(\xi_{k})$, respectively. The PC representation of $u(\bm{\Xi})$, a scalar QoI with finite variance, is given by

[TABLE]

where $\{\psi_{j}(\bm{\Xi})\}$ is a set of multivariate polynomials orthonormal with respect to $f(\bm{\xi})$ and $c_{j}$ are the deterministic PC coefficients to be determined. For computational purpose, the infinite series in (1) may be truncated by retaining the first $P$ terms as

[TABLE]

where $\epsilon(\bm{\Xi})$ represents the random truncation error. A detailed description of the construction of the PC basis $\{\psi_{j}(\bm{\Xi})\}$ and the truncated expansion in (2) will be discussed in Section 2.2.

The main task in PC-based methods is to compute the PC coefficients either intrusively or non-intrusively [40, 72]. In an intrusive approach, the governing equations are projected onto the subspace spanned by the PC basis via the Galerkin formulation [40], often requiring some modifications to the existing deterministic solvers that may not be desirable. On the other hand, non-intrusive methods are based on sampling and treat existing deterministic solvers as a black box. Sampling-based techniques such as least squares approximation (LSA) [51, 44], pseudo-spectral collocation [72, 21], Monte Carlo sampling [72], and compressive sampling [27, 46] have been employed in the literature to identify the PC coefficients.

In this work, we consider the LSA approach to determine the PC coefficients. From the statistical point of view, LSA is similar to the regression of the exact solution $u(\bm{\Xi})$ in the PC basis [51]. Let us denote the $i$th realization of $\bm{\Xi}$ by $\bm{\xi}^{(i)}$, a.k.a a design point, generated via random sampling, for instance, according to $f(\bm{\xi})$, and the corresponding realization of ${u}(\bm{\Xi})$ by ${u}(\bm{\xi}^{(i)})$. Given the pair $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ and $\{{u}(\bm{\xi}^{(i)})\}_{i=1}^{N}$, with $N$ being the number of independent samples, the discrete representation of (2) can be written as

[TABLE]

where $\bm{u}=(u({\bm{\xi}}^{(1)}),\cdots,u({\bm{\xi}}^{(N)}))^{T}\in\mathbb{R}^{N}$ contains the realizations of the QoI $u$, $\bm{\Psi}(i,j)=\psi_{j}({\bm{\xi}}^{(i)})\in\mathbb{R}^{N\times P}$ is the measurement matrix containing samples of the PC basis, $\bm{c}=(c_{1},\cdots,c_{P})^{T}\in\mathbb{R}^{P}$ is the vector of unknown PC coefficients, a.k.a. the estimators, and the vector $\bm{\epsilon}\in\mathbb{R}^{N}$ contains unknown truncation errors between the solution realizations $u$ and the truncated PC approximations of $u$. A diagonal positive-definite matrix $\bm{W}$ is also introduced such that $\bm{W}(i,i)=w(\bm{\xi}^{(i)})$ is a function of the sample points, which in turn depends on the sampling strategy, and will be discussed later in Section 4. The unknown coefficients $\bm{c}$ may be approximated by solving the least squares problem

[TABLE]

where $\|\cdot\|_{2}$ is the standard Euclidean norm. The solution to (4) is computed from the normal equation

[TABLE]

The computational cost of constructing the PCE via (4) is controlled by the number of samples $N$ and may be high when evaluating the QoI is expensive. In general, a stable solution $\bm{c}$ to (4) requires $N\geq P$ (often $N\gg P$) samples of $u$ [44] and, as we shall see in Section 2.2, $P$ grows quickly when either the dimensionality of the problem, $d$, or the total order of the PCE, $p$, is increased. Consequently, accurate approximation of $\bm{c}$ with as small $N$ as possible becomes a critical task. This has motivated researchers to explore various strategies for constructing the sample set $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ beyond the standard MC sampling according to $f(\bm{\xi})$. In general, one may categorize the sampling techniques for the LSA problem into two groups: random and deterministic. Unlike random sampling, deterministic sampling techniques attempt to select the sample points $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ in a deterministic fashion.

Among the early work on sampling techniques for PCE via LSA, we mention [17] in which Latin hypercube (LH) sampling is used instead of the standard MC sampling. Recently, an importance sampling technique referred to as asymptotic was studied in [44, 38, 87] based on asymptotic (in order) analysis of $d$-dimensional Legendre and Hermite polynomials. This approach leads to sampling from a $d$-dimensional Chebyshev distribution for the case of Legendre polynomials and sampling uniformly from a $d$-dimensional ball (with a radius determined by the total order of the PC expansion) for Hermite polynomials. Additionally, for arbitrary dimension and order, a different importance sampling approach dubbed coherence-optimal was proposed in [44] which ensures an accurate computation of $\bm{c}$ with a number of solution realizations that depends linearly on $P$ (up to a logarithmic factor). A related sampling method has also been presented in [19]. Random sampling from tensor product of Gaussian quadrature points has been also proposed in [133], where the authors prove an asymptotic (in order) stability of this method with a number of points which scales linearly (up to a logarithmic factor) with $P$.

A quasi-Monte Carlo (QMC) method using the so-called low-discrepancy points was employed in [79] to deterministically select the sample points for the least squares PCE such that a stable and accurate solution is achieved with a number of samples $N$ proportional to the square of the size of PC basis, $P$, up to logarithmic factors. Another deterministic sampling approach based on Weil’s theorem was proposed in [134] which leads to a stable approximation in the Chebyshev basis with a number of samples that scales quadratically in $P$.

Optimal sampling of the design points $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ for general regression models has received extensive attention in statistics community [36, 37, 9, 95, 3] and is referred to as the optimal design of experiments (ODE). The alphabetic optimality criteria for ODE, a.k.a. classical optimality, either focus on minimizing the error in estimating the unknown coefficients, or parameters, $\bm{c}$, e.g., D-, A-, and E-optimal designs, or the error in the prediction of $u$, e.g., I-optimal design. As explained in Section 4.5, these optimality criterion are mostly based on some functional of the information matrix $\bm{M}$ given by

[TABLE]

when the standard LSA is considered, and

[TABLE]

for the case of weighted LSA problem (4). For example, D-optimal criterion takes the determinant of $\bm{M}^{-1}$, denoted by $|\bm{M}^{-1}|$, and constructs $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ such that $|\bm{M}^{-1}|$ is minimized. Among the limited work utilizing ODE for PCE, we mention [131] in which D-optimality criterion was used to generate the sample set. It was shown that a Fedorov-genetic algorithm used to generate the D-optimal sample set could compute the PC coefficients accurately with a smaller cost in comparison to other common sampling techniques. In [12], $D$-optimal criterion has been employed to efficiently estimate the PCE-based Sobol’ indices via the LSA approach. Very recently, a quasi-optimal sampling approach for the weighted LSA problem has been proposed in [109] which takes a function that involves $|\bm{M}|$ as the optimality measure.

Given the growing interest in designing optimal sampling strategies for PCE’s, the aims of this manuscript are three fold. Firstly, we provide a comprehensive review of both standard and advanced (random and deterministic) techniques for choosing samples $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$. Secondly, we propose the so-called alphabetic-coherence-optimal sampling technique, a hybrid scheme which combines the coherence-optimal approach of [44] with classical alphabetic optimality criteria. Thirdly, we provide a comparison study of the performance of various sampling strategies on three example problems. Our ultimate goal is to provide practitioners with a collection of available options and a set of comparative studies as how these various sampling techniques perform.

The remainder of this paper is organized as follows. Section 2 summarizes the problem formulation, where we start with an abstract problem and provide details on the construction of PC basis. Following that, in Section 3, some definitions and preliminaries used in describing the sampling techniques are presented. Standard MC, asymptotic, coherence-optimal and random quadrature sampling techniques, as well as space-filling methods, such as quasi-MC and LH sampling, are reviewed in Section 4. Moreover, Section 4 provides a comprehensive review on the application of various alphabetic optimality criteria in the context of least squares PCE’s. Finally, three numerical examples concerning the recovery of randomly generated functions, a nonlinear stochastic Duffing oscillator, and the estimation of the remaining useful life of Lithium-ion batteries are considered in Section 5 to study the empirical performance of a number of sampling methods.

2 Problem statement and polynomial chaos expansion

2.1 Problem statement

In the probabilistic framework, model uncertainties are represented using random variables and processes, leading to non-deterministic governing equations,

[TABLE]

Here, $\mathcal{R}$ denotes the set of linear or non-linear governing equations, $\bm{x}$ the spatial variable, $t$ time, $\bm{\Xi}$ the set of $d$ independent random variables characterizing model uncertainties, as introduced in Section 1, and $u$ the solution. One of the goals in uncertainty propagation is to construct an approximation of the mapping $\bm{\Xi}\rightarrow u(\bm{x},t,\bm{\Xi})$, here, using a set of samples $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ and the corresponding solution realizations $\{u(\bm{x},t,\bm{\xi}^{(i)})\}_{i=1}^{N}$, or samples of a QoI depending on $u$. For the interest of a simpler presentation, we assume such a construction is desired for a fixed spatial location $\bm{x}$ and time instance $t$, and hereafter suppress the dependence of $u$ on $\bm{x}$ and $t$.

2.2 Polynomial chaos expansion

PCE provides a framework for approximating the solution $u$ in (8) of a stochastic system, by projecting it onto a basis of polynomials of the random inputs. It was first proposed by Wiener [122] and introduced to the engineering field by Ghanem and Spanos [40] for problems with Gaussian random inputs, and later extended to non-Gaussian random inputs by using polynomials of the Askey scheme (generalized PCE) [127]. The PC representation of the finite variance QoI, $u$, is given in (1) with the associated truncated expansion in (2). To specify the truncation, we first describe the construction of the multi-dimensional basis functions $\psi_{j}(\bm{\Xi})$ adopted here. Let $\{\psi_{j_{k}}(\Xi_{k}):\ j_{k}\in\mathbb{N}_{0}\}$, with $\mathbb{N}_{0}:=\mathbb{N}\cup\{0\}$, denote a complete set of univariate polynomials of order (degree) $j_{k}$ orthonormal with respect to $f(\xi_{k})$; See Table 1. The multi-dimensional basis functions are (tensor) products of the univariate functions, i.e.,

[TABLE]

where the multi-index $\bm{j}=(j_{1},\cdots,j_{d})$ specifies the order of $\psi_{\bm{j}}(\bm{\Xi})$ in each of the dimensions $\Xi_{k}$. For convenience, in (1) we have considered an integer indexing of each $\psi_{\bm{j}}(\bm{\Xi})$ so that there is a one-to-one correspondence between the elements of $\{\psi_{\bm{j}}(\bm{\Xi})\}$ and $\{\psi_{j}(\bm{\Xi})\}$. For computation, in (2), we consider a truncation of (1) where we limit the total order of $\{\psi_{\bm{j}}(\bm{\Xi})\}$ to $p$; that is, we consider only multi-indices ${\bm{j}}$ for which $\sum_{k=1}^{d}j_{k}\leq p$. With this truncation, the number $P$ of basis functions retained in (2) is given by the factorial relation

[TABLE]

The truncated PC expansion in (2) converges in the mean-square sense as $p\ \ (\text{hence}\ \ P)\rightarrow\infty$, when $u$ has finite variance and the coefficients $c_{\bm{i}}$ are computed from the projection equation

[TABLE]

In practice, it may not be feasible or efficient to compute $c_{j}$ directly from (11) and alternative approaches, such as LSA in (4), may be preferred.

3 Background on stability and convergence of least squares polynomial chaos expansion

In order to keep this paper self-contained, we next present a brief summary of results on stability and convergence of the solution to least squares problem (4), highlighting the role of some key quantities.

3.1 Stability of least squares approximation

One of sampling design goals is to improve the stability of (4), in that small perturbations in the samples of $u$ do not lead to large variations in the computed $\bm{c}$ and, consequently, the PC approximation of $u$. Following [18, 44], a key stability quantity is how close $\bm{M}$ is to the identity, given by $\|\bm{M}-\bm{I}\|$, where $\|\cdot\|$ denotes the spectral norm of a matrix. Due to the orthonormality of the PC basis, $\bm{M}$ is on average identity; however, it deviates from the identity for finite samples sizes $N$. A small $\|\bm{M}-\bm{I}\|$ ensures $\bm{M}$ is invertible and has small condition number, and that (5) leads to a stable approximation of $\bm{c}$. These statements are quantified in the following corollary, reporting results from [97, 18, 44].

Corollary 1**.**

Let $\hat{u}(\bm{\Xi})=\mathop{\sum}\limits_{j=1}^{P}\hat{c}_{j}\psi_{j}(\bm{\Xi})$. Let $\hat{\bm{c}}$ and $\hat{\bm{c}}_{\delta}$ be the solution to (5) corresponding to data $\bm{u}$ and $\bm{u}+\delta\bm{u}$, respectively, where $\delta\bm{u}$ is a perturbation in $\bm{u}$. Let $\|\bm{M}-\bm{I}\|\leq\delta$ for some $0\leq\delta<1$. Then,

[TABLE]

where $\mathcal{K}(\bm{M})$ is the condition number of $\bm{M}$. Additionally,

[TABLE]

Proof.

See A for the derivation of (12)-(15). ∎

In particular, (14) and (15) indicate that small changes in $\bm{u}$ do not lead to large changes in $\hat{\bm{c}}$ – hence the stability of the least squares problem – as long as $\delta$ is small. As we shall explain in Section 4, the results in Corollary 1 form the basis of designing sampling strategies for PCE via LSA.

3.2 Coherence parameter definition

Let $w(\bm{\xi})$ denote a weight function depending on the probability density function $f_{\bm{Y}}(\bm{\xi})$ used to generate $\bm{\xi}$, possibly different from $f(\bm{\xi})$, which we shall specify in Sections 4.2 and 4.3. We note that, for the sake of a simpler notation, we refer to all realized input random vectors by $\bm{\xi}$, regardless of the sampling distribution used to generate these realizations. In [44], the coherence parameter $\mu$ given by

[TABLE]

is used to bound (probabilistically) the deviation of $\bm{M}$ from the identity for stability purposes as discussed in Section 3.1. In particular, [18] shows that for $0\leq\delta<1$,

[TABLE]

where $c_{\delta}=\delta+(1-\delta)\log(1-\delta)$ and $\mathbb{P}(\cdot)$ denotes the probability of an event. Notice that (17) suggests bounding $\|\bm{M}-\bm{I}\|$ with high probability requires a sampling rate that linearly depends on $\mu$. More precisely, the following theorem, adopted from [44, Section 2.3], demonstrates the dependence of sample size $N$ on the coherence parameter $\mu$ – or indirectly the sampling distribution – for the convergence of PCE with bounded polynomials, e.g., of Legendre type.

Theorem 1** ([44]).**

Let $\hat{u}(\bm{\Xi})=\mathop{\sum}\limits_{j=1}^{P}\hat{c}_{j}\psi_{j}(\bm{\Xi})$, where $\hat{\bm{c}}=(\hat{c}_{1},\cdots,\hat{c}_{P})^{T}$ is the solution to (5). For $\mathcal{E}$, a sampling event that occurs with probability

[TABLE]

it follows that

[TABLE]

where

[TABLE]

denotes the expectation restricted to the event $\mathcal{E}$.

When $\mu$ in (16) is infinite, i.e., for unbounded polynomials such as Hermite polynomials, a coherence parameter defined on an appropriate subset of $\Omega$ is instead considered. The interested reader is referred to [44, Section 2.2] for further details.

4 Sampling techniques

In this section, we start reviewing the standard MC technique followed by the coherence-optimal and asymptotic sampling methods presented in Sections 4.2 and 4.3, respectively. Similar to MC sampling, the last two sampling strategies are random. However, unlike standard MC sampling, they seek sampling distributions, possibly different from $f(\bm{\xi})$, that lead to some optimality. We then turn our attention to methods that directly generate discrete sample points $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$, starting by stochastic collocation via deterministic and randomized quadrature points in Section 4.4. In Section 4.5, we provide a detailed description of the alphabetic optimal designs and present a sequential technique to construct these designs. Finally, space-filling designs, such as QMC and LH sampling techniques as two instances of widely-used, random sampling strategies in designing computer experiments are presented in Section 4.6.

4.1 Standard Monte Carlo sampling

The natural method of constructing the sample set $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$, known as the standard MC sampling, is to generate $\bm{\xi}$ according to the orthogonality measure $f(\bm{\xi})$ of the PC basis, for which the weights are $w(\bm{\xi})=1$. For instance, standard MC sampling of the $d$-dimensional Legendre polynomials corresponds to sampling from the uniform distribution on $[-1,1]^{d}$. For other polynomials, the corresponding sampling distributions are listed in Table 1.

As shown in [44, Theorem 3.1], a standard MC sampling of the $d$-dimensional Legendre polynomials of total order $p$ gives a coherence parameter $\mu\leq\exp(2p)$, while $\mu$ associated with standard MC sampling of the $d$-dimensional Hermite polynomials is bounded by $\mu\leq C_{p}\eta_{p}^{p}$, with $C_{p}$ and $\eta_{p}$ constants depending on $p$. Notice that these bounds are independent of the dimension $d$, thus suggesting the use of standard MC sampling when $d$ is large and the QoI admits a low order PCE. For high-order PCE’s, $\mu$ is large and therefore $N$ needs to be large, according to Theorem 1. This is why, for those cases, alternative sampling distributions are needed to reduce $N$. Convergence analyses of the standard MC method for univariate and multivariate Legendre polynomials are also presented in [80] and [81], respectively.

4.2 Coherence-optimal sampling

For a fixed sample size $N$, the coherence parameter $\mu$ plays a key role in the stability and convergence of least squares PC approximation, as indicated in (17) – along with Corollary 1 – and Theorem 1. Heuristically, a smaller $\mu$ results in a more stable and accurate approximation. This has motivated the design of a random sampling strategy, dubbed coherence-optimal in [44], which seeks to find a sampling measure – instead of discrete sample points $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ – to minimize $\mu$. Specifically, given a fixed PC basis $\{\psi_{j}(\bm{\xi})\}$, define

[TABLE]

Notice that $B^{2}(\bm{\xi})$ represent a uniformly least upper bound on the sum of the squares of the PC basis functions. As discussed in [45, 19], sampling inputs according to the alternative measure

[TABLE]

with the weight function $w(\bm{\xi})=c^{-1}B(\bm{\xi})^{-1}$ leads to the minimum possible $\mu=P$ compared to any other sampling measure. In (20), $c^{2}=\left(\int_{\Omega}f(\bm{\xi})B^{2}(\bm{\xi})d\bm{\xi}\right)^{-1}=P^{-1}$ is a normalization constant. An analytic expression for $B(\bm{\xi})$ is not generally available; however, for a given set of PC basis, $B(\bm{\xi})$ can be evaluated at any arbitrary values of $\bm{\xi}$. This has motivated using a Markov Chain Monte Carlo approach to sample from $f_{\bm{Y}}(\bm{\xi})$ in (20). A detailed description of this technique can be found in [44, Section 4.3.1] and [45, Algorithm 1].

Following [44], the coherence-optimal sampling approach ensures a stable computation of $\bm{c}$ with a number of solution realizations that depends linearly (up to a logarithmic factor) on the number of PC coefficients $P$, i.e., $N\sim\mathcal{O}(P\log P)$. Moreover, through the application of this approach to various numerical examples, it was empirically observed that the coherence-optimal sampling leads to either similar or considerably more accurate results in comparison to sampling from $f(\bm{\xi})$ [44].

4.3 Asymptotic sampling

For asymptotically large $p$, the coherence-optimal density $f_{\bm{Y}}(\bm{\xi})$ associated with Legendre or Hermite polynomials, respectively, converges to or can be approximated with known distributions, [45, 44]. Specifically, for the case of Legendre polynomials, the asymptotic density coincides with that of i.i.d. Chebyshev random variables over $[-1,1]$, i.e.,

[TABLE]

with the corresponding weights given by

[TABLE]

The coherence parameter of this asymptotic sampling is bounded by $\mu\leq 3^{d}$, thus suggesting the advantage of Chebyshev sampling for high order PCE in low enough dimensions, [44, Theorem 3.2].

An asymptotic analysis of the Hermite polynomials which utilizes Hermite Functions suggests sampling uniformly from a $d$-dimensional ball of radius $\sqrt{2}\sqrt{2p+1}$, with weights given by $w(\bm{\xi}):=\exp(-\|\bm{\xi}\|^{2}_{2}/4)$, [44, Section 3.2]. The corresponding coherence parameter is given by $\mu=\mathcal{O}({(2p)^{d/2}}/{\Gamma(d/2+1)})$, with $\Gamma$ being the Gamma function, which features a weaker dependence on $p$ in comparison to $d$. We refer the interested reader to [44, Section 3.2] for a detailed description of the uniform sampling from a $d$-dimensional ball.

4.4 Subsampled quadrature points

A widely used sampling-based approach to compute the PC coefficients is the stochastic collocation (SC) method [74, 126, 5, 20, 21, 61]. The main idea behind the SC technique is to sample the output QoI at particular points, e.g., Gaussian quadrature points, in the stochastic space and then approximate the solution via interpolation or its (integral) statistics via quadrature integration. SC may also be used to compute the PC coefficients by approximating the $d$-dimensional integrals (projection equations) in (11) [21] using tensor product or sparse grids constructed from, for instance, one-dimensional Gaussian quadrature points [21]. This is known as pseudo-spectral approximation.

Although it has been shown that SC based on full tensor product grids is an effective technique for low-dimensional stochastic problems [126, 6], it suffers from the so-called curse-of-dimensionality for large input dimensions $d$: the number of collocation points $N$ in full tensor product grids grows exponentially fast as a function of $d$, i.e., $N=N_{1}^{d}$, where $N_{1}$ is the number of quadrature points in one dimension. For high-dimensional problems, sparse combination of tensor product rules, as for instance introduced by Smolyak [112], have been proposed to alleviate the issue of curse-of-dimensionality [126, 90].

In the context of least squares PCE’s, Zhou et al. [133] proposed to generate $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ from $N$ randomly (uniformly) selected points of a tensor product Gaussian quadrature grid. It was shown in [133] that this sampling approach, referred to as randomized quadratures, is asymptotically (in order) stable with a number of samples which scales linearly (up to a logarithmic factor) with $P$. Very recently, a technique based on deterministic selection of samples from full tensor grids is proposed in [107] as an alternative to the randomized quadrature approach.

4.5 Optimal design of experiments

Planning an experimental procedure in order to determine the relation between a process inputs, a.k.a. factors, affecting the process outputs prior to performing the actual experiment is studied under the subject of design of experiments (DOE). Conducting experiments, either physical or computational, can be very expensive and time-consuming; hence, it is of interest to extract as much information as possible from a given amount of experimental effort, which is addressed in the context of ODE. One of the earliest work on ODE was presented in 1918 by Smith [111], but the major contributions took place a few decades later when Kiefer and co-workers explained the theory behind ODE and proposed a framework for practical implementation of the $D$-optimal designs [11].

In this work we are interested in approximating the map between model inputs – that are random – and output QoIs in a basis of orthogonal polynomials. In particular, this is achieved via the LSA using selected realizations of inputs and outputs. Although in the context of least squares regression – as opposed to least squares approximation (see Section 4.5.2) – a selection of input realizations for building polynomial models with certain optimality goal has been expensively studied over the last few decades; see, e.g., [36, 37, 9, 95, 3]. The following section provides a brief review and interpretation of a major class of model-based ODE, known as alphabetic optimal design, in relation to least squares PC approximation.

4.5.1 Alphabetic optimality criteria

As discussed in Section 3.1, the information matrix $\bm{M}$ plays a key role in the stability of the least squares solution $\hat{\bm{c}}$, described by its deviation from the identity, condition number, or largest singular value. This has motivated the development of several criteria in ODE to choose sample points $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ in a manner that some functional of $\bm{M}$, $\phi(\bm{M})$, known as alphabetic or classical optimality criterion, is minimized [84]. Among the most popular alphabetic optimality criteria, we mention:

$D$-optimality: An optimal design obtained by minimizing the determinant of the inverse of information matrix, i.e., $\phi_{D}=|\bm{M}^{-1}|^{1/P}$.

• 2.

$A$-optimality: Takes the trace of the inverse of the information matrix as the overall measure of the average variance of the estimators. That is, $\phi_{A}=\text{Tr}(\bm{M}^{-1})$ is minimized, where $\text{Tr}(\cdot)$ denotes the trace operator.

• 3.

$E$-optimality: Minimizes the largest eigenvalue $\lambda_{\max}$ of the inverse of the information matrix, i.e., $\phi_{E}=\lambda_{\max}(\bm{M}^{-1})$ is minimized.

• 4.

$K$-optimality [16, 128]: Minimizes the condition number of the information matrix, i.e, $\phi_{K}=\mathcal{K}(\bm{M})$.

• 5.

$I$-optimality: Minimizes $\phi_{I}=\text{Tr}(\tilde{\bm{M}}\bm{M}^{-1})$, where $\tilde{\bm{M}}=\int_{\Omega}\bm{\psi}^{T}(\bm{\xi})\bm{\psi}(\bm{\xi})f(\bm{\xi})d{\bm{\xi}}$ (see for example [48, Page 5] or [85, Page 473]) and $\bm{\psi}(\bm{\xi})$ is a row vector containing the PC basis functions. Due to the orthonormality of the PC basis, $\tilde{\bm{M}}$ is the $P\times P$ identity matrix and, therefore, $I$-optimality is the same as $A$-optimality.

Remark 1**.**

Instead of minimizing $\phi_{D}$, $\phi_{A}$ and $\phi_{E}$, one may alternatively maximize $|\bm{M}|^{1/P}$, $\text{Tr}(\bm{M})$, $\lambda_{\min}(\bm{M})$ to obtain $D$-, $A$- and $E$-optimal designs, respectively [1, 37].

In the ODE literature, D-, A-, and E-optimal designs are referred to as estimation-oriented optimal designs [60], meaning that these designs are focused on precise estimation of the coefficients ${\bm{c}}$. On the other hand, prediction-oriented optimal designs, such as $I$-optimal, a.k.a. Q-optimal [36], or V-optimal [1], aim at reducing the variance of the approximation $\hat{u}$. $I$-optimal designs have received increasing attention in the context of response surface designs, where the ultimate goal is to make predictions [60]. There are more optimality criteria that are labeled by other letters in the literature and have been introduced over the time. A discussion on those optimal measures is beyond the scope of the present work and we refer the interested reader to [66, 37, 95, 3, 60] for more details.

4.5.2 Interpretation of alphabetic optimality criteria in statistical regression setting

We highlight a distinction between the least squares approximation of this work and those in the statistical regression literature, despite both using the same least squares construction as in (5). In the latter setting, the vector $\bm{\epsilon}$ models the uncertainty or noise in the data/measurements and is often modeled by a random vector with known statistics, e.g., i.i.d. zero mean Gaussian vector, and $\bm{\xi}$ is assumed to be deterministic. In the former, however, $\bm{\xi}$ represents the uncertainty in the underlying model and, given the samples $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$, $\bm{\epsilon}$ is a deterministic error that is unknown. The interpretation of the model (3) in the statistical regression setting leads to particular interpretations of the alphabetic optimality criteria. For example, $D$- and $A$-optimal designs minimize the geometric and arithmetic mean of the variation in the estimated coefficients $\hat{\bm{c}}$, respectively, [30]. $I$-optimal design minimizes the average variance of $\hat{u}$ over possible realizations of $\bm{\xi}$. These variations are a consequence of the uncertainty in the data introduced by $\bm{\epsilon}$.

4.5.3 Constructing optimal designs

In Section 1, we categorized ODE under the deterministic sampling techniques. This is in fact true only for limited cases of low-order polynomial models with a small number of inputs for which analytic expressions for the optimal designs are available, see, e.g., [55, 25, 65]. In general, however, the ODE optimization problems are solved numerically, often, involving search through a large, but finite, number of randomly generated candidate designs. The selected optimal design, therefore, depends on the choice of the candidate ensembles. This suggests that one may consider ODE as (quasi-)random sampling.

Heuristic exchange algorithms are among the earliest search methods proposed to construct optimal designs [113, 73]. Exchange algorithms were originally developed for $D$-optimal designs since, due to the availability of an update equation for the determinant of the information matrix (see B), they are computationally more feasible in comparison to other criteria [37]. Moreover, it has been shown that $D$-optimal designs perform well relative to other criteria [3].

The original Fedorov exchange algorithm (FEA) proposed in [36] starts by generating a discrete candidate set which includes $N_{c}>N$ realizations of $\bm{\xi}$, generated for instance according to $f(\bm{\xi})$, and the corresponding $N_{c}\times P$ candidate measurement matrix $\bm{\Psi}_{c}$. Next, a non-singular initial design $\bm{\Psi}$ of size $N\geq P$ is selected randomly from $\bm{\Psi}_{c}$ and FEA iteratively modifies the current design by exchanging a row $\bm{\psi}_{i}$ in $\bm{\Psi}$ with a row $\bm{\psi}_{j}$ in $\bm{\Psi}_{c}$ which corresponds to the maximum reduction in the $D$-optimal measure $\phi_{D}=|\bm{M}^{-1}|^{1/P}$ of $\bm{\Psi}$. Let $\bm{\Psi}_{\text{old}}$ and $\bm{\Psi}_{\text{new}}$ denote the measurement matrix $\bm{\Psi}$ before and after the exchange, respectively. Using (40) the update formula for the determinant is

[TABLE]

where the Fedorov’s $\Delta$-function is given by

[TABLE]

Consequently, when the $D$-optimal design is constructed by maximizing $|\bm{M}|^{1/P}$, the couple $(\bm{\psi}_{i},\bm{\psi}_{j})$ which correspond to the largest value of the Fedorov’s $\Delta$-function are exchanged to provide the largest improvement in the optimality criterion. The algorithm stops when the largest value of the $\Delta$-function is smaller than a user specified threshold, e.g., $10^{-6}$ [36]. Convergence proof of the FEA can be found in [36].

At each iteration of the original FEA only one exchange is performed. A modified version of FEA with multiple exchanges in each iteration was proposed in [22]. Other popular exchange algorithms include DETMAX algorithm [82], Wynn’s algorithm [124], $k$-exchange algorithm [57], $kl$-exchange algorithm [4], and coordinate exchange algorithm of [78] which does not require a candidate set. Performance comparison between different exchange algorithms can be found in [22, 88, 57, 93]. A discussion on details of these methods is beyond the scope of the present work and we refer the interested reader to the above references for further details.

For either high-dimensional problems or when high order PCE’s are required, the number of unknown coefficients $P$ may be very large. Consequently, the candidate measurement matrix should have a large number of rows to avoid sub-optimal solutions. This results in a large number of exchanges to be performed and eventually high computational costs. Moreover, for criteria other than $D$-optimality, there is no update formula similar to (23) for the simultaneous inclusion and exclusion of the design points that can be used to expedite the exchange procedure.

In order to relax the computational complexity of constructing the alphabetic optimal designs, sequential (greedy) augmentation techniques have been developed [28, 114, 109, 12]. In particular, starting from a (non-singular) design with $n<N$ rows, Dykstra proposed to choose the $(n+1)-th$ row from the candidate design $\bm{\Psi}_{c}$ such that the $D$-optimality criterion is maximally improved [28]. As no exchanges are performed, this algorithm is considerably faster than FEA. Moreover, since only rows are being added to the design, one may employ the simple update formula for the trace of $\bm{M}^{-1}$ (see C) to accelerate the construction of $A$-optimal designs via sequential augmentation. In fact, this technique has been employed in [114] to construct $A$-optimal Bayesian designs. There, it was suggested that the sequentially generated optimal designs may be improved further by FEA, but a noticeable improvement was not observed for the problem at hand [114]. A similar sequential algorithm has been used in [109] to construct designs that maximize the optimality criterion $\mathcal{S}(\bm{\Psi}):=\left(\sqrt{|\bm{\Psi}^{T}\bm{\Psi}|}/\prod_{j=1}^{P}\|\bm{\Psi}(:,j)\|_{2}\right)^{1/P}$, where $\bm{\Psi}(:,j)$ is the $j$th column of $\bm{\Psi}$. In that work, when $n<P$ rows are selected, only the first $n$ columns of $\bm{\Psi}$ are considered; otherwise, $\bm{\Psi}$ (hence $\bm{M}$) will be singular. Stated differently the rows and columns of $P$ are increased simultaneously by one until $P$ rows are selected. When $n\geq P$, only new rows are added until $N$ rows are selected. Very recently, an accurate estimation of the PCE-based Sobol’ indices via the LSA technique has been reported in [12], where the $D$-optimal sample set is generated using a similar greedy approach.

In the present work, we employ the sequential augmentation scheme to construct $A$-, $E$-, and $D$-optimal designs as outlined in Algorithm 1. The following points are worth highlighting regarding this algorithm:

As new rows are added to the design, the computation of $D$- and $A$-optimality criteria may be accelerated by employing the update formulas given in (40) and (41), respectively. There is no update formula for the $E$-optimal designs, hence, its construction is considerably more expensive.

• 2.

Repeating samples are not allowed as they lead to singular designs. Hence, it is crucial to search among the remaining rows of $\bm{\Psi}_{c}$ when adding new rows to $\bm{\Psi}$.

• 3.

In general, the exact solution to alphabetic optimal design is NP-hard [73] and thus infeasible for large $P$. A design obtained by the greedy Algorithm 1 may therefore be sub-optimal. Moreover, the alphabetic optimality criteria $\phi$ are non-convex and typically exhibit multiple local minima/maxima. In order to increase the chance of finding the global optima, it is advised to run them multiple times with different (random) initial designs.

In addition to the greedy algorithms discussed here, other techniques such as simulated annealing (SA) [43, 77], genetic algorithms (GA) [83, 10], semi-definite programming (SDP) [99, 2] and pattern search [50, 48] have been employed in the literature to construct alphabetic optimal designs. Moreover, software packages such as JMP [41], POBE [75], AlgDesign [121], PopED [39, 91] and Gosset [49] have been developed to construct ODE using the algorithms discussed here. We note that the majority of the software packages for the construction of ODE are developed for scenarios where low-order (mostly quadratic) polynomials are employed, hence, limiting the applicability of such packages for high-order PCE’s.

4.6 Space-filling designs

Space-filling designs form a class of DOEs that have received extensive attention in the context of computer experiments [106, 67, 94, 62]. The main idea behind space-filling designs, as suggested by the name, is to evenly spread input samples over the sample space $\Omega$. This is motivated by the heuristic that when the realization of a QoI at an input parameter is deterministic (this may not hold for physical experiments), evaluating the QoI at closely-spaced points may not lead to new information about the QoI, particularly in the scenarios where only a small number of realizations can be afforded. In the following, we will review two important space-filling designs, namely quasi-Monte Carlo (QMC) and Latin hypercube (LH) sampling techniques.

4.6.1 Quasi-Monte Carlo sampling

To alleviate the slow convergence rate of standard MC technique, especially for high-dimensional integrations in finance applications [120], researchers have proposed a deterministic version of the MC sampling method based on the low-discrepancy sequences [26] known as quasi-Monte Carlo (QMC) technique [89, 14, 8, 79]. It has been shown that for integrating functions of $d$-dimensional uniform inputs, a convergence rate of approximately $\mathcal{O}(N^{-1}(\log N)^{d})$, which is asymptotically better than that of MC, can be expected from the QMC technique [14].

Unlike pseudo-random sequences used in standard MC, low-discrepancy sequences, a.k.a. quasi-random, such as Hammersley set, Niederreiter, Halton, Faure, and Sobol sequences [76], are constructed such that the maximal degree of uniformity is provided for the sample points drawn from a uniform distribution. A quantitative assessment of the degree of uniformity of a set of sample points is achieved by calculating its discrepancy. Among the many definitions proposed for the discrepancy of a sequence, star discrepancy $D_{N}^{*}$, defined as follows, is one of the most widely used ones due to its simplicity.

Star discrepancy is used in the context of low-discrepancy sequences for quasi-Monte Carlo sampling techniques. Given the sample set $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ in $d$-dimensional unit cube $[0,1)^{d}$, the star discrepancy $D_{N}^{*}$ of the set $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ is defined as

[TABLE]

where $\mathcal{J}^{*}$ is the class of all subintervals $J$ of $[0,1)^{d}$ given by

[TABLE]

$\#\mathcal{A}$ denotes the number of elements of the set $\mathcal{A}$, and $\lambda_{d}(J)$ is the Lebesgue measure of $J$ [29]. A well uniformly distributed set of points has a small star discrepancy [89, 14].

In [79], a stable and accurate solution to the LSA problem in (4) is achieved with low-discrepancy sample points and a number of samples $N$ which is proportional to the square of $P$ (up to a logarithmic factor). To the best of our knowledge, the application of QMC sampling technique for least squares PCE’s has not been extended to problems with non-uniform random inputs. For a detailed description of QMC sampling approach for multivariate Legendre polynomials we refer the interested reader to [79].

Another class of space-filling designs, known as uniform design, which facilitate the same goal of constructing samples with maximum degree of uniformity was proposed in [32]. Uniform designs are constructed based on alternative and computationally less expensive measures of the uniformity rather than the star discrepancy $D_{N}^{*}$ used in QMC [31]. Several methods such as good lattice [35], optimization searching [123], collapsing [33], combinatorial construction [132], and Latin square [34] have been proposed to construct uniform designs.

4.6.2 Latin hypercube sampling

Latin hypercube (LH) sampling, first introduced in [70], aims at selecting the sample set $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ such that these points cover all portions of the design space $\Omega$. LH sampling attempts to reduce the variance of standard MC estimators by selecting $N$ sample points from $N$ equiprobable partitions (hypercubes) of the sample space, which is known as stratification of the probability distribution. In [115, 92], it is shown that the convergence rate of the LH sampling is never worse than that of standard MC when $N\rightarrow\infty$. However, under some conditions discussed in [70], e.g., when certain monotonicity conditions do not hold, LH may not be advantageous over the standard MC in reducing the variance of an estimator. Algorithm 2 summarizes the main steps to generate $\{\bm{\xi}^{(i)}\}_{i=1}^{N}$ via the LH sampling technique.

Several modified versions of LH sampling such as orthogonal [129, 130, 13] and column-wise orthogonal [116], sliced and optimal sliced [96], and orthogonal-maximin [64] have been proposed to improve the efficiency and accuracy of the LH sampling approach. Algorithm 2 is used when the random inputs are independent from each other. For the correlated random inputs, some extended/modified versions of LH sampling have also been proposed in the literature [53, 52, 101].

In addition to the space-filling designs we reviewed here, many other have been introduced over the past few decades. Among them we mention maximin designs [56], a.k.a. sphere packing designs, in which the minimum distance between pairs of designs points is maximized, maximum entropy designs [108], which takes the entropy as a measure of the uniformity of the samples, Gaussian process integrated mean squared-error (IMSE) designs minimizing the variance of prediction [54, 98], and maximum projection designs [63].

Remark 2**.**

Although alphabetic optimal designs are developed for real world physical experiments, but these designs have also been successfully employed for computer experiments. As an example, we mention the work in [59], which compares $D$- and $I$-optimal designs with space-filling designs such as Latin hypercube, uniform, sphere packing, and maximum entropy, in terms of reducing the prediction variance for high-order polynomials. It was shown that $D$- and $I$-optimal designs result in prediction variances smaller than those obtained by space-filling designs. A similar comparison was also made in [58].

5 Numerical examples

In this section, we study the empirical performance of six different sampling methods by considering three numerical examples. In the first example, the recovery of manufactured PCE’s is investigated. Next, a nonlinear Duffing oscillator under free vibration which requires high order PC expansions is considered. Finally, we examine the application of a subset of the reviewed sampling techniques to the prediction of the remaining useful life of a Lithium-ion battery.

We consider standard MC (Section 4.1), coherence-optimal (Section 4.2) and Latin hypercube (Section 4.6.2) techniques denoted by Standard, coh-opt and LH sampling, respectively. In addition, we introduce hybrid sampling methods, denoted by alphabetic-coherence-optimal, that combine coherence-optimal technique with alphabetic optimality criteria. To construct these sample sets, the $N_{c}$ realizations of $\bm{\Xi}$ in the third line of Algorithm 1 are generated via the coherence-optimal technique of Section 4.2. Consequently, the corresponding $N_{c}\times P$ candidate measurement matrix $\bm{\Psi}_{c}$ to be used in the sequential construction of $D$-, $A$-, and $E$-optimal designs will be already optimal in a sense that it yields to a small coherence parameter $\mu$, hence the names D-coh-opt, A-coh-opt, and E-coh-opt, respectively. We note that the larger $N_{c}$ in Algorithm 1, the higher the quality of the constructed design and the associated computational cost.

To verify the accuracy of the least squares PC model constructed with $N$ independent samples, we compute the relative error of the resulting PC expansion in predicting independent realizations of the QoI as

[TABLE]

where $\bm{u}_{v}$ is the vector of $N_{v}$ additional realizations of QoI (not used in computing $\hat{\bm{c}}$) and $\bm{\Psi}_{v}$ is the measurement matrix corresponding to $\bm{u}_{v}$.

Remark 3**.**

Because of the random nature of the discussed sampling techniques, the results presented in the following are obtained by averaging 60 independent replications using $N$ samples to compute $\bm{c}$ and $N_{v}$ realizations to evaluate the relative error in (27).

5.1 Manufactured functions

In the first example, we consider the reconstruction of manufactured PC expansions, where the coefficients $\bm{c}$ are prescribed to generate $\bm{u}$. Specifically, we consider the model

[TABLE]

where each $c_{j}$ is generated independently from a standard normal distribution. Each realization $u(\bm{\xi}^{(i)})$ is then generated by evaluating the right-hand-side of (28) at a set of $N$ samples of $\bm{\xi}$. Additionally, we assume that there is independent, normally distributed additive noise in the evaluation of ${u}(\bm{\xi}^{(i)})$ with a standard deviation of $0.03\cdot|{u}(\bm{\xi}^{(i)})|$. The PC expansion is then recovered successfully if the least squares solution yields a relative error $\leq 0.02$.

Figs. 1 and 2 represent the probability of a successful recovery of the manufactured solution as a function of the number of samples $N$ to the number of unknowns $P$ for the case of Legendre and Hermite polynomials, receptively. Three different scenarios, i.e., a high-dimensional case with $(p,d)=(2,15)$, a moderate case of $(p,d)=(4,4)$, and a high-order case with $(p,d)=(15,2)$ are considered in these figures. In this example, we consider $N_{c}=4N$ for each value of the oversampling ratio $N/P$. The following observations are worth highlighting regarding Figs. 1 and 2:

For the high-order case $(p,d)=(15,2)$, both standard and LH sampling techniques demonstrate a poor performance in recovering the solution such that for the case of Hermite polynomials, they both fail to recover the solution even for large sample sizes, e.g., $N=10P$. These two methods perform significantly better for the high-dimensional, low-order case $(p,d)=(2,15)$, but are yet worse than the other sampling techniques for the moderate case $(p,d)=(4,4)$.

• 2.

Coherence-optimal approach results in a more accurate solution recovery in comparison to the standard and LH sampling methods, especially, for the high-order case $(p,d)=(15,2)$.

• 3.

In all cases, alphabetic-coherence-optimal techniques outperform the other sampling methods. They are significantly better than the standard and LH techniques for the high-order case $(p,d)=(15,2)$. Their advantage becomes more significant when the oversampling ratio $N/P$ is small.

• 4.

Among the alphabetic-coherence-optimal techniques, D-coh-opt is slightly better than A-coh-opt, while E-coh-opt is the worst of three.

• 5.

The main difference between the alphabetic-coherence-optimal technique and the alphabetic optimal design is that in the former approach the candidate ensemble $\bm{\Psi}_{c}$ is generated based on the coherence-optimal technique while in the latter $\bm{\Psi}_{c}$ is constructed via random sampling from the orthogonality measure $f(\bm{\xi})$. To highlight the importance of constructing $\bm{\Psi}_{c}$ via coherence-optimal approach, we also report the results of the classical $D$-optimal design, denoted by D-opt in Figs. 1 and 2, in which the candidate measurement matrix of Algorithm 1 is constructed by drawing samples from $f(\bm{\xi})$. Although D-opt sampling performs fairly good for the the high-dimensional case in comparison to the alphabetic-coherence-optimal techniques, its accuracy is drastically reduced for the higher order expansions considered.

5.2 Nonlinear Duffing oscillator

The second problem of interest in this study is to quantify the uncertainty in the displacement solution $u(\bm{\Xi},t)$ of a nonlinear single-degree-of-freedom Duffing oscillator [71] under free vibration described by,

[TABLE]

where the uncertain parameters $\{\omega_{i}\}_{i=1}^{3}$, influencing the motion of the oscillator, are given by

[TABLE]

with $\{\Xi_{i}\}_{i=1}^{3}$ being i.i.d. uniform random variables $U(-1,1)$.

While this is a relatively low-dimensional problem ($d=3$), due to the presence of nonlinearities, high order PC expansions are required in order to maintain the approximation’s accuracy at large time instances $t$. This high-order approximation requirement is illustrated in Fig. 3, where we report the mean and standard deviation of the relative error of PCE of $u(\bm{\Xi},4)$ obtained using the standard MC samples. In this experiment, we assumed $N=20P$ to generate sample size independent solutions. We considered an additional $N_{v}=10000$ validation samples to compute the relative error in (27).

We considered $p=9$ for this problem to examine the performance of various sampling strategies. The pair $(p,d)=(9,3)$ leads to $P=220$ PC coefficients to be approximated.

Fig. 4 demonstrates the mean and standard deviation of the relative error in predicting $N_{v}=10000$ independent realizations of the displacement obtained by different sampling methods. Three oversampling ratios of 1.1, two and three which correspond to $N=242,440$ and $660$ samples, respectively, have been considered. Following [87], we assume $N_{c}=\lfloor 1.5P\log(P)\rfloor$, where $P$ is the number of unknown PC coefficients and is given in (10).

Again, the long-time integration issue we discussed earlier is observed here as well; that is, independent of the sampling technique, the relative error slightly increases over time for a fixed $p$. One possible approach to address this issue is to employ Multi-Element generalized Polynomial Chaos (ME-gPC) in which the total order of expansion is adjusted over the time. Discussion on ME-gPC is beyond the scope of this paper and we refer the interested reader to [119].

As it is shown in Fig. 4, the alphabetic-coherence-optimal techniques result in smaller values for the relative error in comparison to the other methods, with D-coh-opt and A-coh-opt to perform slightly better than E-coh-opt. This superiority is signified in Figs. 4a-b where the oversampling ratio is small. These plots suggest that for the case of high-order PC expansions, the alphabetic-coherence-optimal techniques are the most suitable sampling methods, in particular, when the number of samples is slightly larger than the number of unknown coefficients $P$. For the applications where evaluation of the QoI is computationally expensive, we prefer to approximate the PCE with a minimum number of solution realizations, usually with $N$ slightly larger than $P$.

The standard MC and LH techniques provide the least accurate results, specially, for the case of $N=242$ where the oversampling ratio is low. In this case, these two methods fail to accurately approximate the solution as the error features large mean and variance (Fig. 4b). As the number of samples is increased, MC and LH methods perform considerably better such that with $N=440$ samples, both techniques result in acceptable errors.

Similar to the previous example, the errors obtained by the coherence-optimal sampling falls somewhere between those of MC (or LH) and the alphabetic-coherence-optimal sampling methods. By comparing Figs. 4a, 4c, and 4e, it can be seen that, as $N$ is increased, the coherence-optimal results converge to the alphabetic-coherence-optimal results at a higher rate in comparison to MC and LH methods.

5.3 Prediction of remaining useful life of batteries

The last problem of interest in this study is the model-based estimation of the remaining useful life (RUL) of a Lithium-ion battery (LIB). RUL of an LIB is defined as the amount of time before the battery health falls below a defined threshold, e.g., reaching the cut-off potential at the end of discharge [102], and is studied under the context of prognostics and health management [118]. For the sake of completeness, a brief introduction to model-based prognostics of LIB’s with an emphasis on RUL predictions is presented in the following.

5.3.1 Model-based prognostics

In general, prognosis approaches may be categorized as data-driven and model-based [105, 7]. A mathematical representation of the system based on the physics of the problem is used in model-based methods, while data-driven techniques use system monitoring data instead [105]. Both techniques have been successfully employed to predict the RUL of batteries [110, 24, 68, 104, 103, 102, 69]. Here, we focus on model-based RUL estimation for LIB’s which is suitable for online health monitoring.

Let us assume that the system model is given by the following set of equations, [24, 104, 103, 102],

[TABLE]

where $\bm{z}(t)$ is the state vector, $\bm{f}$ the state equation, $\bm{\theta}(t)$ the model parameter vector, $\bm{\nu}(t)$ the input vector, $\bm{v}_{p}(t)$ the process noise vector, $\bm{y}(t)$ the output vector, $\bm{h}$ the output equation, and $\bm{v}_{m}(t)$ the measurement noise vector.

As illustrated in [102, Fig. 1], the computational framework for prognostics and RUL predictions includes three main steps: estimation of the state $\bm{z}$ at time $t_{p}$, prediction of the future states, and finally RUL computation. Kalman or particle filtering approaches may be employed within a Bayesian framework to estimate the state vector $\bm{z}$ using the output data measured up to time $t_{p}$. Both Eqs. (31) and (32) should be used in the estimation step. Following [102], in this study we assume that the estimated state vector $\bm{z}$ at time $t_{p}$ is given, and instead we focus on the prediction and RUL computation steps.

In the prediction step, only the state space model in (31) is used since no measured data is available for $t>t_{p}$, and the goal is to predict an event $E$ (depending on the states, parameters, and inputs) indicating the LIB’s end of life. In the present work, $E$ occurs when the LIB’s terminal voltage $V$ is smaller than a user-specified cut-off voltage $V_{\text{cutoff}}$. More specifically, for each $t_{p}$, we define

[TABLE]

and RUL as

[TABLE]

5.3.2 An electrical equivalent circuit model for LIB’s

In order to facilitate the application of the model-based RUL estimation for the online health monitoring of LIB’s, the mathematical model described by Eqs. (31) and (32) is desired to be computationally inexpensive. Following closely [104, 103, 102], we consider an empirical model of the LIB based on the electrical equivalent circuit model shown in Fig. 5. The charge of the LIB, denoted by $q_{b}$, is held by the nonlinear capacitance $C_{b}$. The voltage drop due to the surface over-potential is given by $R_{sp}-C_{sp}$; see Fig. 5. Additionally, the resistances $R_{s}$ and $R_{p}$ are used to model the Ohmic drop and the parasitic self-discharging, respectively, [104, 103, 102].

Denoting the LIB maximum capacity and charge by $C_{max}$ and $q_{max}$, respectively, the state of charge ($SoC$), is given by

[TABLE]

[103]. The concentration polarization resistance $R_{sp}$ is a non-linear function of $SoC$,

[TABLE]

where $R_{sp0}$, $R_{sp1}$, and $R_{sp2}$ are empirical parameters [104, 103, 102]. $C_{b}$ is also given by the nonlinear (in $SoC$) relation

[TABLE]

with $C_{bi}$, $i=0,1,2,3$, empirically determined constants [104, 103, 102]. The model parameters of the considered cell are taken from [104, 103, 102] and shown in Table 2.

Voltage drops across the individual circuit elements and their associated currents and charges are governed by the set of equations presented in Table 3, where $I$ denotes the discharge current. The terminal voltage of the cell $V$ is given by $V=V_{b}-V_{sp}-V_{s}$. We assume that in this example $V_{\text{cutoff}}=16\text{V}$.

Different uncertainty sources such as those in the future inputs, state, model and process noise have been identified for the stochastic RUL estimations [105, 24, 100, 23, 104, 103, 102]. We refer the readers to the provided references for more details on the sources of uncertainty for this problem. In this example, we assume that the battery is discharged at a constant current (in Amps) represented by a random variable following a Beta distribution with shape parameters $\alpha=21.2$ and $\beta=31.8$. As mentioned earlier, we assume that the state estimates are already obtained. For this example, we consider an estimation of the state variables, $q_{b}$, $q_{sp}$, and $q_{s}$, given in [102, Figs. 5-7], respectively, with a coefficient of variation (COV) of 0.1. Whenever the estimated state is equal to zero, we assume its standard deviation is equal to 0.1. In addition, three process noise terms with zero mean normal distributions and variances of 0.1, $10^{-4}$, and $10^{-6}$ are considered in the state equations describing $q_{b}$, $q_{sp}$, and $q_{s}$, respectively. Assuming that there is no model uncertainty, the total stochastic dimension of this problem is therefore $d=7$. We note that the univariate PC basis functions $\psi_{{j}_{k}}(\Xi_{k})$ for this problem consist of Jacobi and Hermite polynomials which correspond to Beta and normal random variables, respectively.

Fig. 6 shows the relative error of the sampling techniques in estimating the RUL of LIB at $t=0,200,400$, and 600 seconds. Similar to the previous example, 10000 independent samples are used to compute the relative error. Since model-based RUL predictions are mainly used in online applications where fast computations are critical, we limit the number of sample points in this example to be $N=P+1$, which corresponds to 36 and 121 samples for the second and third order PCE’s, respectively. Similar to the previous example, we assumed $N_{c}=\lfloor 1.5P\log(P)\rfloor$. Fig. 6 demonstrates that similar to the high-order PC expansions, alphabetic-coherence-optimal techniques perform better that other methods for the high-dimensional problems as well when the oversampling ratio is low. Moreover, it can be seen that for the case of $d>p$ with low oversampling ratios, coherence-optimal method is not necessarily advantageous over the standard MC and LH samplings.

It is shown in Fig. 6c-d that the D-coh-opt approach results in the smallest errors. We next compare the pdf of the RUL predicted with a $p=3$ PCE, obtained using the D-coh-opt approach, against a MC-based pdf in Fig. 7. Although the pdf’s obtained by these two different approaches coincide, we note that $20\times 10^{6}$ solution realizations are used in constructing the MC pdf, while D-coh-opt-based PCE only requires 121 samples. This highlights the efficacy of the PCE approach, equipped with optimal sampling methods, for the online RUL prediction of LIB’s.

6 Conclusion and summary

We provided an extensive review of the sampling techniques available in the literature for PCE’s via least squares approximation. We also proposed a hybrid sampling approach referred to as alphabetic-coherence-optimal in which, given the same sample size, a combination of alphabetic-optimal designs, such as D-, A-, and E-optimal, and the recently developed coherence-optimal technique results in more accurate approximations of the PC coefficients in comparison to other sampling methods.

In order to provide insight for practitioners seeking the best sampling method in building PCE’s via least squares approximation, we compared the empirical performance of several sampling techniques through their application to three problems involving uniform, Gaussian, and Beta random inputs, and for both low and high order (and dimension) expansions. The highlights of our observations, backed up by theoretical results, are summarized below:

High-order PCE with $p>d$: We observed that for problems in which high-order expansions are required to meet accuracy requirements, the alphabetic-coherence-optimal techniques deliver the most accurate approximations. Among the alphabetic-coherence-optimal techniques, D-coh-opt results were slightly better than those obtained by A-coh-opt. In our experiments, we observed that the construction of the D-coh-opt sample sets is considerably faster than that of the A-coh-opt samples; hence, we recommend D-coh-opt for high-order PCE’s. We note that generating input samples according to the alphabetic-coherence-optimal design is more expensive than sampling the inputs, for instance, randomly according to their joint pdf. However, in most practical engineering applications, generating realizations of the output QoI is considerably expensive, so the extra computational cost for sampling the inputs may be ignored.

• 2.

High-dimensional PCE with $d>p$: Although the alphabetic-coherence-optimal techniques still perform better than other methods considered here for the case of $d>p$, the standard MC and LH sampling methods may also deliver accurate approximations when low-order expansions are sufficient. In such cases, one may employ standard MC, LH, or other advanced space-filling sampling techniques discussed in this study as they may be less computationally expensive to generate than the alphabetic-coherence-optimal approach.

• 3.

The oversampling ratio $N/P$: In the numerical examples provided in this study, it was observed that when $N$ is 3-4 times larger than the number of PC coefficients $P$, one can expect that the standard MC and LH sampling methods exhibit a fairly acceptable performance as long as the order $p$ of the expansion is small. When high order expansions are required, even with a fairly large oversampling ratio, these methods may fail to provide an accurate PC approximation. For the cases with small values of $N/P$, alphabetic-coherence-optimal techniques outperform other methods considerably, specially when high-order PCE’s are employed.

Finally, we note that existing algorithms for solving alphabetic-optimal design problems are computationally expensive for large values of $P$, thus limiting their applicability to relatively low-order and/or low-dimensional PCE’s. Developing scalable optimization techniques to address this difficulty remains an open problem in ODE.

Acknowledgements

MH’s work was supported by the National Science Foundation grant CMMI-1201207. This material is based upon work of AD supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, under Award Number DE-SC0006402 and NSF grant CMMI-1454601. AD is thankful of Prof. Paul Constantine from Colorado School of Mines for fruitful discussions on least squares approximation and regression.

\appendixpage

Appendix A Derivation of bounds in Corollary 1

Using the assumption $\|\bm{M}-\bm{I}\|\leq\delta$, $0\leq\delta<1$, and the symmetry of $\bm{M}-\bm{I}$, we have that

[TABLE]

This implies $|\lambda(\bm{M}-\bm{I})|\leq\delta$, which in turn gives $\lambda_{\max}(\bm{M})-1\leq\delta$ and $-\delta\leq\lambda_{\min}(\bm{M})-1$; hence, the inequality in (12). A direct application of (12) gives (13).

To see (14), we first note that $\|\hat{u}\|_{L_{2}(\Omega,f)}^{2}=\mathbb{E}[\hat{u}^{2}]=\|\hat{\bm{c}}\|^{2}$ due to the orthonormality of the PC basis. From (5) and following [18], we have

[TABLE]

To bound $\|\left(\bm{W}\bm{\Psi}\right)^{T}\|$, we first observe that $\|\left(\bm{W}\bm{\Psi}\right)^{T}\|=\|\bm{W}\bm{\Psi}\|=\|(\bm{W}\bm{\Psi})^{T}(\bm{W}\bm{\Psi})\|^{1/2}$. Notice that $\|1/N(\bm{W}\bm{\Psi})^{T}(\bm{W}\bm{\Psi})\|\leq\|1/N(\bm{W}\bm{\Psi})^{T}(\bm{W}\bm{\Psi})-\bm{I}\|+1\leq(1+\delta)$. Therefore,

[TABLE]

To bound $\|\left((\bm{W}\bm{\Psi})^{T}(\bm{W}\bm{\Psi})\right)^{-1}\|$, we follow a similar analysis as in [97, Section 7.2]. Specifically,

[TABLE]

which implies

[TABLE]

Putting together (38) and (39), we arrive at the inequality (13).

Notice that (15) is a standard result for which we refer the reader to [47, Section 3.5].

Appendix B Matrix determinant update formula

Let $\bm{A}$ be an invertible matrix of size $P\times P$, and let $\bm{b}$ and $\bm{a}$ be two vectors of size $P$. Then

[TABLE]

Proof.

When $\bm{A}=\bm{I}$, we have

[TABLE]

Taking the determinant of both sides results in $|\bm{I}\pm\bm{b}\bm{a}^{T}|=1\pm\bm{a}^{T}\bm{b}$. Finally, one may write $|\bm{A}\pm\bm{b}\bm{a}^{T}|=|\bm{A}(\bm{I}\pm\bm{A}^{-1}\bm{y}\bm{a}^{T})|=|\bm{A}|(1\pm\bm{a}^{T}{\bm{A}}^{-1}\bm{b})$. ∎

Appendix C Matrix trace update formula

Let $\bm{A}$ be an invertible matrix of size $P\times P$, and let $\bm{b}$ and $\bm{a}$ be two vectors of size $P$. Additionally, assume $1\pm\bm{a}^{T}\bm{A}^{-1}\bm{b}\neq 0$. Then

[TABLE]

Proof.

Sherman-Morrison rank-one update formula states

[TABLE]

which can be easily verified by substituting (42) in identity $(\bm{A}\pm\bm{b}\bm{a}^{T})(\bm{A}\pm\bm{b}\bm{a}^{T})^{-1}=\bm{I}$. Taking the trace of both sides of (42) completes the proof. ∎