Approximation and Uncertainty Quantification of Systems with Arbitrary Parameter Distributions Using Weighted Leja Interpolation

Abstract

Approximation and uncertainty quantification methods based on Lagrange interpolation are typically abandoned in cases where the probability distributions of one or more system parameters are not normal, uniform, or closely related distributions, due to the computational issues that arise when one wishes to define interpolation nodes for general distributions. This paper examines the use of the recently introduced weighted Leja nodes for that purpose. Weighted Leja interpolation rules are presented, along with a dimension-adaptive sparse interpolation algorithm, to be employed in the case of high-dimensional input uncertainty. The performance and reliability of the suggested approach is verified by four numerical experiments, where the respective models feature extreme value and truncated normal parameter distributions. Furthermore, the suggested approach is compared with a well-established polynomial chaos method and found to be either comparable or superior in terms of approximation and statistics estimation accuracy.

1. Introduction

Studying physical systems often involves assessing the influence of random parameters upon a system’s behavior and outputs. Such uncertainty quantification (UQ) studies are becoming increasingly popular in a variety of technical domains. Of particular importance to industrial scientists and engineers are so-called non-intrusive UQ methods, where the term “non-intrusive” signifies that the often complex and typical proprietary simulation software applied to model physical systems are incorporated into UQ studies without any modifications.

Sampling methods such as (quasi-) Monte Carlo (MC) [1] or latin hypercube sampling (LHS) [2] are the most common approaches applied for non-intrusive UQ, as they are easy to implement, straightforwardly parallelizable, and can be applied to problems with a large number of parameters under relatively mild assumptions regarding the regularity of the underlying mathematical model. However, their slow convergence rates render them unattractive or even intractable in cases where high accuracy is demanded and single simulations are highly time-consuming.

Computationally efficient alternatives are spectral UQ methods [3,4], which approximate the functional dependence of a system’s quantity of interest (QoI) on its input parameters. The most popular black-box methods of this category employ global polynomial approximations by aid of either Lagrange interpolation schemes [5,6] or generalized polynomial chaos (gPC) [7], where the latter approach is typically based on least squares (LS) regression [8,9,10] or pseudo-spectral projection methods [11,12,13], or, less commonly, on interpolation [14]. Assuming a smooth input-output dependence, spectral UQ methods provide fast convergence rates, in some cases even of exponential order. However, their performance decreases rapidly for an increasing number of input parameters due to the curse of dimensionality [15]. To overcome this bottleneck, state-of-the-art approaches rely on algorithms for the adaptive construction of sparse approximations [8,9,10,16,17,18].

Adaptive sparse approximation algorithms for Lagrange interpolation methods typically employ nested sequences of univariate interpolation nodes. While not strictly necessary [19], nested node sequences are helpful for the efficient construction of sparse grid interpolation algorithms. Nested node sequences are readily available for the well studied cases of uniformly or normally distributed random variables (RVs) [20,21], but not for more general probability measures. In principle, one could derive nested interpolation (or quadrature) rules tailored to an arbitrary probability density function (PDF) [22,23], however, this is a rather cumbersome task for which dedicated analyses are necessary each and every time a new PDF is considered.

Diversely, the main requirement in gPC approximations, either full or sparse, is to find a complete set of polynomials that are orthogonal to each other with respect to the input parameter PDFs. In the case of arbitrary PDFs, such polynomials can be numerically constructed [24,25,26]. As a result, despite the fact that according to a number of studies Lagrange interpolation methods enjoy an advantage over gPC-based methods with respect to the error-cost ratio [27,28,29], they are abandoned whenever the input parameter PDFs are not uniform, or normal, or closely related distributions, e.g., log-uniform or log-normal. Consequently, the literature lacks comparisons between Lagrange interpolations and gPC-based approximations for other probability measures.

This work investigates the use of Lagrange interpolation methods based on the recently introduced weighted Leja node sequences [30] for the approximation of systems with arbitrarily distributed random parameters. Weighted Leja sequences provide interpolation and quadrature nodes tailored to an arbitrary PDF in a seamless way. Moreover, being nested by definition, they allow for the efficient construction of sparse grids, thus enabling high-dimensional UQ studies. We note that this paper focuses on continuous PDFs of arbitrary shapes and does not consider data-driven approaches which define a PDF through its moments without any further assumptions about it [31,32].

Due to a number of desirable properties, e.g., granularity, interpolation, and quadrature stability, and nestedness, weighted Leja nodes have been getting increasing attention in the context of approximation and UQ [16,30,33,34,35,36]. Nonetheless, with the exceptions of [34,37,38] where beta distributions are considered for the parameters, the use of Leja nodes has been limited to uniform and log-normal parameter distributions. The contribution of this paper is exactly to start filling this gap with the application of weighted Leja interpolation for possibly high-dimensional systems with arbitrary parameter PDFs.

To that end, we present here the use of weighted Leja rules in the context of sparse Lagrange interpolation for approximation and UQ purposes in a comprehensive and self-consistent way. Section 2 summarizes the main theory related to weighted Leja interpolation. In Section 2.1 we first present a general overview of the UQ setting in which Leja interpolation will be applied. Then, in Section 2.2 we recall the basics of univariate Lagrange interpolation, both standard and hierarchical. Unweighted and weighted Leja rules along with their properties are presented in Section 2.3, while in in Section 2.4 the extension of Lagrange interpolation to higher dimensions, along with a dimension-adaptive, Leja-based scheme for the construction of sparse approximations and its use for UQ purposes, is presented. Numerical experiments which verify the suitability of the Leja interpolation scheme for physical systems with arbitrarily distributed stochastic input parameters are presented in Section 3. In the same section, the adaptive weighted Leja interpolation UQ method is compared against a well-known adaptive gPC algorithm based on least angle regression (LAR). We conclude with a discussion on the overall findings of this work and on possibilities for further research in Section 4.

2. Theory and Methodology

2.1. Stochastic Parametric Models

Without loss of generality, let us assume that the behavior of the physical system under consideration can be modeled by a set of parametric partial differential equations (PDEs). The general form of the corresponding mathematical model is given by:

$${\mathcal{D}}_{\mathbf{y}}\left(u\right)=0,$$

(1)

where $\mathcal{D}={\mathcal{D}}_{\mathbf{y}}\left(u\right)$ is a differential operator, $\mathbf{y}\in \mathrm{\Xi }\subset {\mathbb{R}}^N$ a N-dimensional parameter vector, $u=u_{\mathbf{y}}$ the parameter-dependent solution of Equation (1), and $\mathrm{\Xi }$ a measurable set called the image space. Typically, we are interested in a certain quantity, that is the QoI, which is computed by post-processing the solution $u=u_{\mathbf{y}}$. For simplicity, we assume that the QoI is a real scalar, however, the methodology discussed in Section 2.4 can be applied to complex or vector-valued QoIs with only minor modifications. We further assume that the functional dependence between the parameter vector and the QoI is given by the map $g:\mathbf{y}\mapsto g\left(\mathbf{y}\right)$, which is assumed to be deterministic, i.e., the exact same output $g\left(\mathbf{y}\right)$ is observed each and every time the QoI is evaluated for the same parameter vector $\mathbf{y}$.

Parametric problems in the form of Equation (1) arise in both deterministic and stochastic settings. An example of the former case is, e.g., an optimization study, while an example of the latter is, e.g., a UQ study as considered in this work. In the stochastic setting, the parameter vector $\mathbf{y}$ corresponds to a realization of a N-dimensional RV $\mathbf{Y}=\left(Y_1,Y_2,\cdots ,Y_N\right)$, alternatively called a random vector. Assuming that the RVs $Y_n$, $n=1,\cdots ,N$, are mutually independent, the joint PDF is given as a product of univariate PDFs, such that $\varrho \left(\mathbf{y}\right)={\prod }_{n=1}^N{\varrho }_n\left(y_n\right).$ While the case of dependent RVs is not rigorously addressed here, we mention the option of using generalized Nataf or Rosenblatt transformations to transform the dependent RVs into independent RVs [39,40]. The sparse interpolation method presented in Section 2.4 can then be applied to the transformed RVs.

Due to its dependence on stochastic input parameters, the deterministic map g results now in a random output. In other words, the input uncertainty propagates through the deterministic model and renders the QoI uncertain. The main goal of UQ studies is to quantify this output uncertainty by computing statistics, such as moments, event probabilities, or sensitivity metrics. Denoting the expectation operator by $\mathbb{E}\left[·\right]$, any statistic can be written in the general form:

$$\mathbb{E}\left[\varphi \left(g\right)\right]={\int }_{\mathrm{\Xi }}\varphi \left(g\left(\mathbf{y}\right)\right)\varrho \left(\mathbf{y}\right)\mathrm{d}\mathbf{y},$$

(2)

where the functional $\varphi$ corresponds to the statistic of interest, e.g. $\varphi \left(g\right)=g$ for the expected value of the QoI and $\varphi \left(g\right)={\left(g-\mathbb{E}\left[g\right]\right)}^2$ for its variance.

The integral of Equation (2) cannot in general be computed analytically. Therefore, one relies on numerical integration methods, such as MC integration or quadrature. Considering the case of MC integration, if the original input-output map g corresponds to an expensive numerical model, then it is desirable to replace it by an inexpensive albeit accurate surrogate model $\tilde{g}\approx g$. In this work, this surrogate is obtained by interpolating g. To capture accurately the behavior of $g\left(\mathbf{y}\right)$, given that $\mathbf{y}$ are realizations of a random vector, the interpolation must be based on nodes tailored to the PDF $\varrho \left(\mathbf{y}\right)$. The surrogate model can also be used for tasks other than computing statistics, e.g., for optimization. Moreover, by integrating the interpolant, one can derive a quadrature scheme for the computation of Equation (2). Similar to the interpolation, nodes tailored to the PDF are crucial for the convergence of the quadrature. Finally, in the case of a high-dimensional random vector $\mathbf{Y}$, sparse quadrature and interpolation schemes must be used in order to mitigate the effect of the curse of dimensionality. Such sparse schemes most commonly rely on nested univariate node sequences (see Section 2.2).

As pointed out in Section 1, there exist interpolation nodes that meet the aforementioned requirements for the well studied cases of uniform and normal PDFs, e.g., Clenshaw–Curtis and Genz–Keister nodes, respectively [20,21]. We note that only a finite number of Genz–Keister nodes are available, thus imposing an additional limitation on its achievable accuracy. However, constructing suitable interpolation and quadrature rules for more exotic or even arbitrary PDFs, while in principle possible [22,23], remains a computational challenge. This problem is here circumvented with the use of weighted Leja nodes, discussed in Section 2.3, which are easy to compute, by definition nested, and tailored to a continuous PDF.

2.2. Univariate Interpolation Schemes

Let us assume a univariate input-output map $g\left(y\right)$, where y denotes a realization of a single RV Y. We let the non-negative integer $i\in {\mathbb{Z}}^{\ge 0}$ denote the interpolation level, $Z_i={\left\{y_{i,j}\right\}}_{j=1}^{m\left(i\right)}$ the level-i grid of interpolation nodes, and $m:{\mathbb{Z}}^{\ge 0}\to {\mathbb{Z}}^{+}$, $m\left(0\right)=1$, a strictly monotonically increasing “level-to-nodes” function, relating the interpolation level to the grid size. Based on these three elements, in the following subsections we define Lagrange and hierarchical interpolation rules, as well as related quadrature schemes.

2.2.1. Univariate Lagrange Interpolation

The level-i univariate Lagrange interpolation is defined as:

$${\mathcal{I}}_i\left[g\right]\left(y\right)=\sum _{j=1}^{m\left(i\right)}g\left(y_{i,j}\right)l_{i,j}\left(y\right),$$

(3)

where $l_{i,j}$ are univariate Lagrange polynomials, defined on the interpolation grid $Z_i$ as:

$$l_{i,j}\left(y\right)=\prod _{k=1,k\ne j}^{m\left(i\right)}\frac{y-y_{i,k}}{y_{i,j}-y_{i,k}}.$$

The accuracy of the interpolation depends crucially on the choice of nodes, e.g., Gauss–Legendre and Gauss–Hermite quadrature nodes are popular choices for RVs following uniform or normal distributions, respectively. In particular, the error between g and ${\mathcal{I}}_i\left[g\right]$ is related to the best polynomial approximation error ${ϵ}_{\mathrm{best}}$ through [6,16]:

$${∥g-{\mathcal{I}}_i\left[g\right]∥}_{L^{\infty }\left({\mathbb{P}}_i\right)}\le \left(1+{\mathbb{L}}_i\right){ϵ}_{\mathrm{best}},$$

where ${\mathbb{P}}_i$ denotes the space of Lagrange polynomials with degree up to $m\left(i\right)-1$ and ${\mathbb{L}}_i$ the Lebesgue constant:

$${\mathbb{L}}_i=\sup \frac{{∥{\mathcal{I}}_i\left[g\right]∥}_{L^{\infty }\left({\mathbb{P}}_i\right)}}{{∥g∥}_{L^{\infty }\left({\mathbb{P}}_i\right)}}.$$

Therefore, assuming that the interpolation error decreases for increasing interpolation levels, the interpolation nodes should be associated to a Lebesgue constant which grows at a slower rate.

2.2.2. Hierarchical Univariate Interpolation

We further demand that the interpolation nodes form sequences of nested grids for increasing interpolation levels, such that $Z_{i-1}\subset Z_i$, $\forall i>0$. In that case, we may define a hierarchical interpolation scheme as follows. We first define the operator ${\mathsf{\Delta }}_i$ as the difference between two consecutive univariate interpolants, i.e.,

$${\mathsf{\Delta }}_i={\mathcal{I}}_i-{\mathcal{I}}_{i-1},$$

(4)

where ${\mathcal{I}}_{-1}\left[g\right]\left(y\right)=0$. The interpolation in Equation (3) can now be given as a sum of interpolant differences, such that:

$${\mathcal{I}}_i\left[g\right]\left(y\right)=\sum _{k=0}^i{\mathsf{\Delta }}_k\left[g\right]\left(y\right).$$

(5)

The hierarchy of formula in Equation (5) is obvious if we consider two consecutive interpolation levels, $i-1$ and i, in which case the level-i interpolation reads:

$$\begin{array}{}\mathrm{(6a)}& {\mathcal{I}}_i\left[g\right]\left(y\right)\hfill & ={\mathcal{I}}_{i-1}\left[g\right]\left(y\right)+\sum _{j:y_{i,j}\in Z_i{Z}_{i-1}}\left(g\left(y_{i,j}\right)-{\mathcal{I}}_{i-1}\left[g\right]\left(y_{i,j}\right)\right)l_{i,j}\left(y\right)\hfill \mathrm{(6b)}& & ={\mathcal{I}}_{i-1}\left[g\right]\left(y\right)+\sum _{j:y_{i,j}\in Z_i{Z}_{i-1}}{s}_{i,j}{l}_{i,j}\left(y\right),\hfill \end{array}$$

where the interpolation coefficients $s_{i,j}$, given by:

$$s_{i,j}=g\left(y_{i,j}\right)-{\mathcal{I}}_{i-1}\left[g\right]\left(y_{i,j}\right),y_{i,j}\in Z_i{Z}_{i-1},$$

(7)

are called hierarchical surpluses. In the particular case of $m\left(i\right)=i+1$, which can in fact be accomplished using the Leja sequences discussed in Section 2.3, a single grid point is added for each new interpolation level. Then, the index j can be dropped and Equation (6) is simplified to ${\mathcal{I}}_i\left[g\right]\left(y\right)={\mathcal{I}}_{i-1}\left[g\right]\left(y\right)+s_i{l}_i\left(y\right)$.

The obvious advantage of using nested grids and the hierarchical Equations (5) or (6), is that at each new level i the QoI must be evaluated only at the new nodes $y_{i,j}\in Z_i{Z}_{i-1}$. Moreover, univariate hierarchical interpolation rules constitute the backbone of the sparse adaptive multivariate interpolation algorithms discussed in Section 2.4.

2.2.3. Interpolatory Univariate Quadrature

Considering the UQ goal of estimating specific statistics given by Equation (2), we may use an available interpolation ${\mathcal{I}}_i\left[g\right]$ to derive an appropriate quadrature scheme. For example, assuming an univariate interpolation given by:

$${\mathcal{I}}_i\left[g\right]\left(y\right)=\sum _{j=1}^{m\left(i\right)}{s}_{i,j}{l}_{i,j}\left(y\right),$$

the expected value of the QoI can be estimated as:

$$\mathbb{E}\left[g\right]\approx \mathbb{E}\left[{\mathcal{I}}_i\left[g\right]\right]=\sum _{j=1}^{m\left(i\right)}{s}_{i,j}{w}_{i,j},$$

(8)

where $w_{i,j}$ are quadrature weights given by:

$$w_{i,j}=\mathbb{E}\left[l_{i,j}\right]={\int }_{\mathrm{\Xi }}{l}_{i,j}\left(y\right)\varrho \left(y\right)\mathrm{d}y.$$

Quadrature schemes for the estimation of statistics other than the expected value can be derived in a similar way, by applying the expectation operator along with the associated functional $\varphi$.

2.3. Leja Interpolation Nodes

As pointed out in Section 2.1 and Section 2.2, the univariate interpolation nodes employed in Lagrange interpolation schemes must satisfy certain requirements. First, for the interpolation to become more accurate with increasing levels, the Lebesgue constant associated with the interpolation nodes must increase at a rate lower than the decrease in the polynomial approximation error [6,16]. Second, the choice of nodes must result in accurate quadrature schemes, similar to Equation (8), to be used in the computation of statistics. Third, interpolation nodes which form nested grids can be used to derive hierarchical interpolation schemes, which result in reduced computations per added interpolation level. Additionally, as discussed in Section 2.4, sparse interpolations depend crucially on hierarchical univariate interpolation rules. Finally and most importantly in the context of this work, the interpolation must be constructed with respect to arbitrary weight functions, i.e., input PDFs. All requirements are addressed here by using weighted Leja sequences [30], discussed in the following.

2.3.1. Unweighted Leja Nodes

A sequence of standard, unweighted Leja nodes [41] ${\left\{y_j\right\}}_{j\ge 0}$, $y_j\in \left[-1.1\right]$, is defined by solving the optimization problem:

$$y_j=\underset{y\in \left[-1,1\right]}{argmax}\prod _{k=0}^{j-1}\left|y-y_k\right|,$$

(9)

for each node $y_j$, $j\ge 1$. The initial node $y_0\in \left[-1.1\right]$ can be chosen arbitrarily, therefore Leja sequences are not unique. Using simple scaling rules, the unweighted Leja sequences defined in Equation (9) are directly applicable to Lagrange interpolation involving constant weight functions, equivalently, uniform PDFs in the UQ context. The extension to non-uniform PDFs is discussed in Section 2.3.2. With respect to the remaining node requirements, it has been shown that the Lebesgue constant associated with unweighted Leja nodes increases subexponentially [42]. Regarding Leja-based quadrature rules, an extensive discussion can be found in [30]. Considering uniform PDFs, a number of works [30,36,37] show that Leja-based quadrature schemes are sufficiently accurate, albeit suboptimal compared to more standard choices, e.g., based on Clenshaw–Curtis nodes. Finally, the optimization problem of Equation (9) results in nested node sequences irrespective of the employed level-to-nodes function $m\left(i\right)$. This is an additional advantage of using Leja nodes, since the nested node sequences can be as granular as the user wishes. In this work, we use $m\left(i\right)=i+1$ to get the minimum of one extra node per interpolation level, i.e., $\#\left(Z_i{Z}_{i-1}\right)=1$. Another popular choice is to use $m\left(i\right)=2i+1$ [36], which is typically employed for the construction of symmetric Leja sequences, where $y_0=0$, $y_1=1$, $y_j$ is given as in Equation (9) for odd $j>1$, and $y_j=-y_{j-1}$ for even j. In all cases, Leja sequences are more granular compared to other nested node sequences, e.g., the level-to-nodes function $m\left(i\right)=2^i+1$ must be used for nested Clenshaw–Curtis nodes.

2.3.2. Weighted Leja Nodes

In more general UQ settings, e.g., for the arbitrary input distributions considered in this work, the definition of Leja sequences must be adjusted according to the input PDF $\varrho \left(y\right)$, which acts as a weight function. To that end, we employ the definition of weighted Leja sequences given in [30]. Using this formulation, weighted Leja sequences are constructed by incorporating the PDF $\varrho \left(y\right)$ into the optimization problem of Equation (9), which is transformed into

$$y_j=\underset{y\in \mathrm{\Xi }}{argmax}\sqrt{\varrho \left(y\right)}\prod _{k=0}^{j-1}\left|y-y_k\right|.$$

(10)

We note that other formulations exist as well [43,44], however, the form given in Equation (10) is preferred due to the fact that the resulting Leja nodes are asymptotically distributed as weighted Gauss quadrature nodes [30]. While this is a promising property, it is not necessarily sufficient to produce an accurate interpolation rule. In some cases, e.g., for weight functions resembling the Gaussian PDF [44], it can be shown that weighted Leja nodes are also associated with a subexponentially growing Lebesgue constant. While more general theoretical results are currently not available, weighted Leja interpolation has been found to perform very well in practice [33,37], as also demonstrated in this paper by the results of the numerical experiments presented in Section 3. The other beneficial properties of Leja sequences, i.e., granularity and nestedness, are preserved in the weighted case. Moreover, since the interpolation rule is now tailored to the specific weight function that is the PDF, the corresponding quadrature rules will be tailored to that PDF as well. For a more detailed discussion on the properties of weighted Leja sequences, see [30].

2.4. Sparse Adaptive Leja Interpolation

The multivariate interpolation is based on suitable combinations of the univariate interpolation rules discussed in Section 2.2. We denote with $\mathbf{i}=\left(i_1,i_2,\cdots ,i_N\right)$ a multi-index of interpolation levels per input RV, with $\mathsf{\Lambda }$ a multi-index set, and with ${\mathsf{\Delta }}_{\mathbf{i}}$ the tensor-product difference operator given by:

$${\mathsf{\Delta }}_{\mathbf{i}}={\mathsf{\Delta }}_{1,i_1}\otimes {\mathsf{\Delta }}_{2,i_2}\otimes \cdots \otimes {\mathsf{\Delta }}_{N,i_N},$$

where the univariate difference operators ${\mathsf{\Delta }}_{n,i_n}$ are defined as in Equation (4), to be recalled in the following subsections. In Section 2.4.1, we discuss interpolation on generalized sparse grids. A dimension-adaptive interpolation scheme based on Leja sequences is presented in Section 2.4.2. Post-processing multivariate interpolations for UQ purposes is discussed in Section 2.4.3.

2.4.1. Generalized Sparse Grid Interpolation

Given a multi-index set $\mathsf{\Lambda }$, the interpolation-based approximation of the multivariate input-output map $g\left(\mathbf{y}\right)$ reads:

$${\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\left(\mathbf{y}\right)=\sum _{\mathbf{i}\in \mathsf{\Lambda }}{\mathsf{\Delta }}_{\mathbf{i}}\left[g\right]\left(\mathbf{y}\right).$$

(11)

We note that Equation (11) is not in general interpolatory, i.e., the interpolation property ${\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\left({\mathbf{y}}_{*}\right)=g\left({\mathbf{y}}_{*}\right)$, where ${\mathbf{y}}_{*}$ is a multivariate interpolation node, does not necessarily hold. The interpolation property holds only if $\mathsf{\Lambda }$ is a downward-closed set and the underlying univariate interpolation rules are based on nested nodes [6]. Downward-closed sets, also known as monotone or lower sets, satisfy the property:

$$\forall \mathbf{i}\in \mathsf{\Lambda }\Rightarrow \mathbf{i}-{\mathbf{e}}_n\in \mathsf{\Lambda },\forall n=1,2,\cdots ,N,\mathrm{with}{i}_n>0,$$

(12)

where ${\mathbf{e}}_n$ denotes the unit vector in the nth dimension. If, additionally, the multi-index set $\mathsf{\Lambda }$ satisfies some sparsity constraint, then Equation (11) is a sparse interpolation formula, while the corresponding sparse grid is given by:

$$Z_{\mathsf{\Lambda }}=\bigcup _{\mathbf{i}\in \mathsf{\Lambda }}{Z}_{\mathbf{i}}=\bigcup _{\mathbf{i}\in \mathsf{\Lambda }}\left(Z_{1,i_1}\times Z_{2,i_2}\times \cdots \times Z_{N,i_N}\right).$$

A commonly used sparsity constraint is:

$$\mathsf{\Lambda }=\left\{\mathbf{i}:\sum _{n=1}^N{i}_n\le i_{max}\right\},$$

resulting in the so-called isotropic sparse grids [5,6]. Compared to full isotropic tensor grids, i.e., those given by:

$$\mathsf{\Lambda }=\left\{\mathbf{i}:\underset{n=1,\cdots ,N}{max}{i}_n\le i_{max}\right\},$$

the complexity of which is $\mathcal{O}\left(m{\left(i_{max}\right)}^N\right)$, isotropic sparse grids delay significantly the curse of dimensionality, since their complexity is reduced to $\mathcal{O}\left(m\left(i_{max}\right){\left(logm\left(i_{max}\right)\right)}^{N-1}\right)$ [45]. Nonetheless, isotropic sparse grids still grow exponentially with respect to the number of parameters.

In most practical cases, the influence of the input parameters upon the QoI is anisotropic. For example, the QoI may be more sensitive to some parameters compared to others, in the sense that variations in the first set of parameters result in comparably greater variations in the QoI. Accordingly, the functional dependence of the QoI on a certain parameter might be significantly more difficult to approximate compared to another, e.g., a highly nonlinear versus a linear relation.

The interpolation of Equation (11) can be constructed such that this anisotropy is represented in its terms, equivalently in the multi-indices comprising the set $\mathsf{\Lambda }$ [16,17]. This is accomplished via an anisotropic refinement of the parameter space according to the contribution of each parameter to the interpolation accuracy. Compared to isotropic schemes, anisotropic interpolations typically result in great computational savings for the same approximation accuracy, or, equivalently, in greater interpolation accuracy for an equal cost. Since, in most cases, the parameter anisotropy cannot be a priori estimated in an accurate way [46], anisotropic interpolations are usually constructed using greedy, adaptive algorithms [18,47,48]. For the particular case of Leja nodes, such an algorithm is discussed in Section 2.4.2. Leja-based, adaptive, anisotropic sparse-grid interpolation algorithms are also available in [16,30].

2.4.2. Adaptive Anisotropic Leja Interpolation

We will base the adaptive construction of anisotropic sparse interpolations on the dimension-adaptive algorithm first presented in [47] for quadrature purposes. Variants of this algorithm appear in a number of later works [16,30,36,37,48]. In this work, we assume the all underlying univariate interpolation rules employ Leja nodes, presented in Section 2.3, along with the level-to-nodes function $m\left(i\right)=i+1$. Our approach is depicted in Algorithm 1, while a detailed presentation follows.

Let us assume that a sparse interpolation ${\mathcal{I}}_{\mathsf{\Lambda }}$, based on a downward-closed multi-index set $\mathsf{\Lambda }$ and on univariate Leja sequences, is given as in Equation (11). We define the set of admissible multi-indices, i.e., those that satisfy the downward-closedness property of Equation (12) if added to $\mathsf{\Lambda }$, as:

$${\mathsf{\Lambda }}^{\mathrm{adm}}=\left\{\mathbf{i}:\mathbf{i}\notin \mathsf{\Lambda }\mathrm{and}\mathsf{\Lambda }\cup \left\{\mathbf{i}\right\}\mathrm{is}\mathrm{downward}-\mathrm{closed}\right\}.$$

Due to the fact that all univariate Leja rules employ the level-to-nodes function $m\left(i\right)=i+1$, each multi-index $\mathbf{i}\in {\mathsf{\Lambda }}^{\mathrm{adm}}$ is associated with a single new interpolation node ${\mathbf{y}}_{\mathbf{i}}=Z_{\mathbf{i}}{Z}_{\mathsf{\Lambda }}$. Therefore, should an admissible multi-index $\mathbf{i}\in {\mathsf{\Lambda }}^{\mathrm{adm}}$ be added to $\mathsf{\Lambda }$, we obtain the hierarchical interpolation scheme:

$${\mathcal{I}}_{\mathsf{\Lambda }\cup \left\{\mathbf{i}\right\}}\left[g\right]\left(\mathbf{y}\right)={\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\left(\mathbf{y}\right)+s_{\mathbf{i}}{L}_{\mathbf{i}}\left(\mathbf{y}\right),$$

where, similarly to Equation (7), the hierarchical surpluses $s_{\mathbf{i}}$ are given by:

$$s_{\mathbf{i}}=g\left({\mathbf{y}}_{\mathbf{i}}\right)-{\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\left({\mathbf{y}}_{\mathbf{i}}\right),\mathbf{i}\in {\mathsf{\Lambda }}^{\mathrm{adm}},$$

and $L_{\mathbf{i}}$ are multivariate Lagrange polynomials given by:

$$L_{\mathbf{i}}\left(\mathbf{y}\right)=\prod _{n=1}^N{l}_{n,i_n}\left(y_n\right),$$

where, due to the use of nested Leja sequences with the level-to-nodes function $m\left(i\right)=i+1$, the univariate Lagrange polynomials can now be defined hierarchically [16], such that,

$$l_{n,i_n}\left(y_n\right)=\prod _{k=0}^{i_n-1}\frac{y_n-y_{n,k}}{y_{n,i_n}-y_{n,k}}.$$

Similar to [47], we expect that coefficients of small magnitudes correspond to terms with relatively small contributions to the available interpolation, which can thus be omitted. Accordingly, terms whose coefficients have large magnitudes are expected to enhance the accuracy of the interpolation and should therefore be added. Based on this argument, we define for each multi-index $\mathbf{i}\in {\mathsf{\Lambda }}^{\mathrm{adm}}$ the contribution indicator:

$${\eta }_{\mathbf{i}}=\left|s_{\mathbf{i}}\right|.$$

Naturally, the multi-index set $\mathsf{\Lambda }$ should be enriched with the admissible multi-index:

$${\mathbf{i}}^{*}=\underset{\mathbf{i}\in {\mathsf{\Lambda }}^{\mathrm{adm}}}{argmax}{\eta }_{\mathbf{i}},$$

that is the index corresponding to the maximum contribution. This procedure can be continued iteratively until a simulation budget B is reached, or until the total contribution of the admissible set is below a tolerance $ϵ$. At any given step of the algorithm, the total number of simulations, equivalently, model evaluations or function calls, is equal to $\#Z_{\mathsf{\Lambda }\cup {\mathsf{\Lambda }}^{\mathrm{adm}}}$. After the termination of the iterations, the final interpolation is constructed using all multi-indices in $\mathsf{\Lambda }\cup {\mathsf{\Lambda }}^{\mathrm{adm}}$, since the hierarchical surpluses corresponding to the admissible nodes have already been computed.

2.4.3. Post-Processing

Once available, the sparse interpolation can be used as an inexpensive surrogate model that can replace the original model in computationally demanding tasks, such as sampling-based estimation of statistics. Specific statistics can be computed directly after the interpolation terms. One possibility is to first transform the Lagrange basis into a gPC basis [49] and then use the well-known formulas of gPC approximations for the estimation of moments or sensitivity indices [8,9]. Another approach is to apply the expectation operator along with the functional corresponding to the statistic directly upon the interpolation, as to derive an appropriate quadrature scheme, similar to the univariate case presented in Section 2.2.3. For example, assuming a Leja-based interpolation given by:

$${\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\left(\mathbf{y}\right)=\sum _{\mathbf{i}\in \mathsf{\Lambda }}{s}_{\mathbf{i}}{L}_{\mathbf{i}}\left(\mathbf{y}\right),$$

the expected value of the QoI can be estimated as:

$$\mathbb{E}\left[g\right]\approx \mathbb{E}\left[{\mathcal{I}}_{\mathsf{\Lambda }}\left[g\right]\right]=\sum _{\mathbf{i}\in \mathsf{\Lambda }}{s}_{\mathbf{i}}{w}_{\mathbf{i}},$$

(13)

where the quadrature weights $w_{\mathbf{i}}$ are given as products of the respective univariate weights, i.e.,

$$w_{\mathbf{i}}=\mathbb{E}\left[L_{\mathbf{i}}\right]=\mathbb{E}\left[\prod _{n=1}^N{l}_{n,i_n}\right]=\prod _{n=1}^N\mathbb{E}\left[l_{n,i_n}\right]=\prod _{n=1}^N{w}_{n,i_n}.$$

Clearly, analytical formulas similar to Equation (13) can be derived for other statistics as well.

3. Results

In this section, we consider models with relatively atypical input distributions, for which the application of interpolation-based UQ methods is not straightforward. In particular, we consider input RVs following truncated normal and Gumbel distributions, the latter also known as the generalized extreme value distribution of type 1 or the log-Weibull distribution. We denote the truncated normal distribution with $\mathcal{TN}\left(\mu ,{\sigma }^2,l,u\right)$, where $\mu$ and ${\sigma }^2$ refer to the mean value and the variance of a normal distribution, while l and u are the lower and upper truncation limits. The Gumbel distribution is denoted with $\mathcal{G}\left(\ell ,\beta \right)$, where ℓ is the location parameter and $\beta$ the scaling parameter. Dedicated nested interpolation nodes are not readily available with respect to the corresponding input PDFs, therefore, interpolation-based UQ methods have not been considered so far in the literature for such random inputs. Here, we compute Leja nodes tailored to truncated normal and Gumbel distributions simply by applying Equation (10) with weight functions that are equal to the corresponding PDFs. Exemplarily, Figure 1 shows the first 10 computed Leja points with respect to two truncated normal and two Gumbel distributions, each with different distribution parameters.

Based on the numerical results presented in this section, we demonstrate that weighted Leja interpolation can reliably be used for such atypical input PDFs, in terms of both approximation and UQ. For further verification, we compare the performance of the weighted Leja interpolation against gPC approximations with numerically constructed orthogonal polynomials [24]. All models feature multiple input parameters, hence, we use adaptive algorithms resulting in sparse approximations. The dimension-adaptive Algorithm 1 based on weighted Leja nodes is implemented in our in-house developed DALI (Dimension-Adaptive Leja Interpolation) software [14,37]. The sparse gPC approximations are constructed with a well-established, degree-adaptive algorithm based on least angle regression (LAR) [9] and implemented in the UQLab [50] software. Quasi-random experimental designs based on Sobol sequences are used in the LAR-gPC approach, while numerically constructed orthogonal polynomials are used to tackle the input PDFs [24].

3.1. Error Metrics

Let us assume that an interpolation or gPC-based polynomial approximation $\tilde{g}\approx g$ is available. The following errors are used to estimate the performance of $\tilde{g}$ for approximation and UQ purposes.

First, we want to estimate the performance of $\tilde{g}$ in terms of approximation accuracy, i.e., its suitability to act as a surrogate model and reliably replace the input-output map g, equivalently, the original computational model. To that end, we use a validation sample ${\left\{{\mathbf{y}}_q\right\}}_{q=1}^Q$, which is randomly drawn from the joint input PDF $\varrho \left(\mathbf{y}\right)$ and compute the root-mean-square (RMS) validation error:

$${ϵ}_{\mathrm{RMS}}=\sqrt{\frac{1}{Q}\sum _{q=1}^Q{\left(\tilde{g}\left({\mathbf{y}}_q\right)-g\left({\mathbf{y}}_q\right)\right)}^2}.$$

(14)

From an interpretation point of view, ${ϵ}_{\mathrm{RMS}}$ indicates the average (expected) approximation accuracy of $\tilde{g}$, while at the same time being sensitive to outliers by penalizing large errors. In all numerical experiments, the size of the validation sample is set to $Q={10}^5$.

Secondly, we want to estimate the performance of $\tilde{g}$ in terms of UQ quality, i.e., estimate the accuracy of the statistics values provided by post-processing the approximation. For that purpose, we use the expected value of the QoI as a representative statistic and compute the relative error:

$${ϵ}_{\mathrm{rel},\mathbb{E}}=\left|\frac{\mathbb{E}\left[g\right]-\widetilde{\mathbb{E}\left[g\right]}}{\mathbb{E}\left[g\right]}\right|,$$

(15)

where $\mathbb{E}\left[g\right]$ is a reference value and $\widetilde{\mathbb{E}\left[g\right]}$ is an estimate computed by directly post-processing the terms of $\tilde{g}$. In all numerical experiments, the reference expected value is computed with a quasi-MC integration scheme based on Sobol sequences, where the size of the quasi-MC sample is equal to ${10}^8$.

3.2. Accuracy versus Costs

In the following numerical experiments, the performance of both polynomial approximations is measured by the respective error-cost relation. In most engineering applications, the cost of running simulations for different sets of parameter values typically outweighs any other numerical operation necessary for the construction of the approximation. In this context, the cost of a method refers solely to the number of simulations, equivalently, solver calls or model evaluations, that are needed until the sought approximation accuracy is reached. This notion of cost is also used in this section, despite the fact that computationally inexpensive models are employed. In particular, in each numerical experiment and for both methods, we compute 19 polynomial approximations corresponding to increasing simulation budgets, i.e., for $B=10,20,\cdots ,100,200,\cdots ,1000$, and compare the errors of Equations (14) and (15) for the same costs.

We should note that regression-based gPC methods, such as the LAR-gPC [9] method employed here for comparison purposes, are significantly more expensive than interpolation-based methods in terms of the neglected costs. This is due to the fact that the solution of possibly expensive LS problems is necessary. On the other hand, regression-based gPC methods allow to pre-compute the QoI for different parameter values, i.e., run the required simulations offline before the algorithm for the construction of the gPC approximation is used. Moreover, this offline task is embarrassingly parallelizable. On the contrary, as can easily be observed from Algorithm 1, dimension-adaptive sparse interpolation is based on sequential computations and therefore does not allow the use of parallelization to that extent. Despite it not being addressed here, those differences must be taken into account by practitioners, considering the problem and the computational resources at hand.

3.3. Borehole Model

We consider the 8-dimensional parametric function:

$$g\left(\mathbf{y}\right)=\frac{2\pi T_{\mathrm{u}}\left(H_{\mathrm{u}}-H_{\mathrm{l}}\right)}{ln\left(\frac{r}{r_{\mathrm{w}}}\right)\left(1+\frac{T_{\mathrm{u}}}{T_{\mathrm{l}}}+\frac{2L{T}_{\mathrm{u}}}{ln\left(\frac{r}{r_{\mathrm{w}}}\right)r_{\mathrm{w}}^2{K}_{\mathrm{w}}}\right)},$$

(16)

which models the water flow through a borehole [51]. The water flow is given in m³/y and the parameter vector is $\mathbf{y}=\left(r_{\mathrm{w}},r,T_{\mathrm{u}},H_{\mathrm{u}},T_{\mathrm{l}},H_{\mathrm{l}},L,K_{\mathrm{w}}\right)$, where $r_{\mathrm{w}}$ and r respectively denote the radius of the borehole and the radius of influence (in m), $T_{\mathrm{u}}$ and $T_{\mathrm{l}}$ the transmissivity of the upper and lower aquifer (in m²/y), $H_{\mathrm{u}}$ and $H_{\mathrm{l}}$ the potentiometric head of the upper and lower aquifer (in m), L the length of the borehole (in m), and $K_{\mathrm{w}}$ the hydraulic conductivity of the borehole (in m/y).

Regarding the parameter distributions, we follow the setting given in [51], where specific value ranges and probability distributions are given for each parameter and use this information to recast all parameters as to follow truncated normal distributions. In this numerical experiment, the input parameters are modeled as follows.

• The parameter $r_{\mathrm{w}}$ originally follows the normal distribution $\mathcal{N}\left(\mu =0.1,\sigma ={0.0161812}^2\right)$. The distribution is now truncated to the range $\left[\mathrm{0.05.0.15}\right]$;
• The parameter r originally follows the log-normal distribution $\mathcal{LN}\left({\mu }_{\mathrm{LN}}=7.71,{\sigma }_{\mathrm{LN}}=1.0056\right)$. Therefore, the parameter r has a mean value equal to $exp\left({\mu }_{\mathrm{LN}}+\frac{{\sigma }_{\mathrm{LN}}^2}{2}\right)$ and a variance equal to $\left(exp\left({\sigma }_{\mathrm{LN}}^2\right)-1\right)exp\left(2{\mu }_{\mathrm{LN}}+{\sigma }_{\mathrm{LN}}^2\right)$. The corresponding truncated normal distribution is defined by these mean and variance values, as well as by the truncation range $\left[100.50000\right]$l
• The remaining parameters are originally uniformly distributed, such that $T_{\mathrm{u}}\in \left[63070,115600\right]$, $H_{\mathrm{u}}\in \left[990,1110\right]$, $T_{\mathrm{l}}\in \left[63.1,116\right]$, $H_{\mathrm{l}}\in \left[700,820\right]$, $L\in \left[1120,1680\right]$, and $K_{\mathrm{w}}\in \left[9855,12045\right]$. Assuming a uniform distribution with support in $\left[a,b\right]$, the corresponding truncated normal distribution is given as $\mathcal{TN}\left(\mu =\frac{a+b}{2},{\sigma }^2=\frac{{(b-a)}^2}{12},l=a,u=b\right)$, i.e., the mean value and variance of the normal distribution correspond to those of the original uniform distribution, while the truncation limits coincide with the uniform distribution’s support boundaries.

The borehole model of Equation (16) is approximated by both the dimension-adaptive weighted Leja interpolation and the degree-adaptive LAR-gPC expansion. Figure 2 shows the error-cost relation for both approximations, with respect to the errors in Equations (14) and (15). Regarding the ${ϵ}_{\mathrm{RMS}}$ error, Leja interpolation outperformed the LAR-gPC method by approximately one order of magnitude over the whole cost-range. Leja interpolation was also advantageous with respect to the ${ϵ}_{\mathrm{rel},\mathbb{E}}$ error, again always outperforming the LAR-gPC method and stagnating much sooner. We note that this error stagnation was observed due to the fact that both methods reached the accuracy of the quasi-MC estimate. Therefore, for the case of the borehole model, the weighted Leja interpolation approach was not only found to perform well for the considered truncated normal input PDFs, but was also superior to the LAR-gPC method for both considered error metrics.

3.4. Steel Column Limit State Function

We consider the 10-dimensional parametric limit state function:

$$g\left(\mathbf{y}\right)=F_{\mathrm{s}}-P_{\mathrm{t}}\left(\frac{1}{2BD}+\frac{F_0{E}_{\mathrm{b}}}{BDH\left(E_{\mathrm{b}}-P_{\mathrm{t}}\right)}\right),$$

(17)

which relates the reliability of a steel column to its cost [51,52]. On the right hand side of Equation (17), $P_{\mathrm{t}}=P_{\mathrm{d}}+P_1+P_2$ is the total load and $E_{\mathrm{b}}=\frac{{\pi }^{2}EBD{H}^2}{2L^2}$ is known as the Euler buckling load. Thereby, the parameter vector $\mathbf{y}$ consists of ten parameters, namely the yield stress $F_{\mathrm{s}}$ (in MPa), the dead weight load $P_{\mathrm{d}}$ (in N), the variable loads $P_1$ and $P_2$ (in N), the flange breadth B (in mm), the flange thickness D (in mm), the profile height H (in mm), the initial deflection $F_0$ (in mm), Young’s modulus E (in MPa), and the column length L (in mm).

In the original setting given in [51,52], $P_1$ and $P_2$ are modeled with Gumbel distributions, E follows a Weibull distribution, $F_{\mathrm{s}}$, B, D, and H are log-normally distributed, and $P_{\mathrm{d}}$ and $F_0$ are normally distributed. In this numerical experiment, we replace all log-normal and normal distributions by truncated normal distributions and use the Gumbel distribution to model the remaining three parameters. The corresponding distributions are given in

Table 1. Parameter distributions of the steel column function.

Parameter	Distribution
$F_{\mathrm{s}}$	$\mathcal{TN}\left(\mu =400,{\sigma }^2={35}^2,l=295,u=505\right)$
$P_{\mathrm{d}}$	$\mathcal{TN}\left(\mu =500000,{\sigma }^2={50000}^2,l=350000,u=650000\right)$
$P_1$	$\mathcal{G}\left(\ell =559495,\beta =70173\right)$
$P_2$	$\mathcal{G}\left(\ell =559495,\beta =70173\right)$
B	$\mathcal{TN}\left(\mu =300,{\sigma }^2=3^2,l=291,u=309\right)$
D	$\mathcal{TN}\left(\mu =20,{\sigma }^2=2^2,l=14,u=26\right)$
H	$\mathcal{TN}\left(\mu =300,{\sigma }^2=5^2,l=285,u=315\right)$
$F_0$	$\mathcal{TN}\left(\mu =30,{\sigma }^2={10}^2,l=0,u=60\right)$
E	$\mathcal{G}\left(\ell =208110,\beta =3275\right)$
L	$\mathcal{TN}\left(\mu =7500,{\sigma }^2={7.5}^2,l=7470,u=7530\right)$

. With the exception of the parameter L, for which the truncation limits are set to $\mu \pm 4\sigma$, all other truncated normal distributions have support in $\left[\mu -3\sigma ,\mu +3\sigma \right]$. For all parameters that follow a truncated normal distribution, $\mu$ and ${\sigma }^2$ coincide with the original mean and variance values given in [51,52]. Accordingly, the location and scale parameters of the Gumbel distributions are chosen such that the original mean values and variances given in [51,52] are preserved.

We approximate the steel column function Equation (17) using the dimension-adaptive weighted Leja interpolation and the degree-adaptive, LAR-based gPC method. For both approximations and for both error metrics in Equations (14) and (15), the error-cost relation is presented in Figure 3. Once more, the Leja interpolation method performs very well for the given input PDFs. However, contrary to the results of Section 3.3, no advantage is observed over the LAR-gPC method. With respect to the ${ϵ}_{\mathrm{RMS}}$ error, the performance of both methods is comparable. This can be probably attributed to a reduced regularity of the QoI defined in Equation (17), in comparison to the QoI given in Equation (16). The Leja interpolation method has a slight edge for costs greater than 200 simulations, however, the difference is minor. With the exception of costs greater than 800 simulations, the LAR-gPC method is superior to the Leja interpolation in terms of the ${ϵ}_{\mathrm{rel},\mathbb{E}}$ error. Overall, the weighted Leja interpolation is again able to yield accurate approximations and statistics and it may be regarded as competitive to the LAR-gPC method.

3.5. Meromorphic Function

We consider the 16-dimensional meromorphic function:

$$g\left(\mathbf{y}\right)=\frac{1}{1+\mathbf{w}·\mathbf{y}},$$

(18)

where $\mathbf{w}$ is a vector of positive weights. The weight vector in Equation (18) is given by $\mathbf{w}=\frac{\widehat{\mathbf{w}}}{2\parallel \widehat{\mathbf{w}}{\parallel }_1}$, where $\widehat{\mathbf{w}}=\left(1,0.5,0.1,0.05,0.01,\cdots ,5·{10}^{-8}\right)$ [10]. The input vector $\mathbf{y}$ takes values in the image space of the random vector $\mathbf{Y}=\left(Y_1,Y_2,\cdots ,Y_{16}\right)$, where the single RVs are given as:

$$\begin{array}{c}\hfill Y_n\sim \left\{\begin{array}{cc}\mathcal{TN}\left(\mu =0,{\sigma }^2=1,l=0,u=3\right),\hfill & n=1,3,\cdots ,15,\hfill \\ \mathcal{TN}\left(\mu =0,{\sigma }^2=1,l=-3,u=0\right),\hfill & n=2,4,\cdots ,16.\hfill \end{array}\right.\end{array}$$

While a normal distribution truncated exactly at its (untruncated) mean value is not very likely to be encountered in practical applications, the selected truncation limits result in very atypical PDFs which fit very well within the concept of arbitrary input PDFs.

Figure 4 shows the error-cost relations corresponding to the error metrics of Equations (14) and (15) for both approximations of the meromorphic function in Equation (18), i.e., using the dimension-adaptive weighted Leja interpolation algorithm and the degree-adaptive LAR-gPC algorithm. The results resemble those of Section 3.3, i.e., the Leja interpolation method had a clear advantage over the LAR-gPC approximation with respect to both error metrics. Regarding the ${ϵ}_{\mathrm{RMS}}$ error, the difference between the two approximations was typically greater than one order of magnitude, while it increased for an increasing computational budget. Regarding the ${ϵ}_{\mathrm{rel},\mathbb{E}}$ error, the two methods showed a comparable performance for costs up to 100 simulations, however, the Leja interpolation was again clearly superior as the simulation budget increased. Hence, in this numerical experiment, the weighted Leja interpolation was able to accurately approximate a model with very atypical input PDFs, and provide accurate estimations of statistics. Moreover, it was found to have a clear edge over the well-established LAR-gPC method.

3.6. Dielectric Inset Waveguide

We consider a rectangular waveguide with a dispersive dielectric inset, similar to the one investigated in [53]. A 2D illustration of the waveguide is depicted in Figure 5, where w denotes its width, ℓ is the length of the dielectric material, and d is a vacuum offset. The waveguide is further extended in the y direction by its height h, which is not shown in Figure 5. The dielectric material has permittivity $\epsilon ={\epsilon }_0{\epsilon }_{\mathrm{r}}$ and permeability $\mu ={\mu }_0{\mu }_{\mathrm{r}}$, where the subscripts “0” and “$\mathrm{r}$” refer to the absolute value of the material property in vacuum and to its relative value for the given dielectric material, respectively. The relative material values are given by Debye relaxation models of second order [54], such that:

$$\begin{array}{cc}\hfill {\epsilon }_{\mathrm{r}}& ={\epsilon }_{\infty }+\frac{{\epsilon }_{\mathrm{s},1}-{\epsilon }_{\infty }}{1+\left(ı\omega {\tau }_{\epsilon ,1}\right)}+\frac{{\epsilon }_{\mathrm{s},2}-{\epsilon }_{\infty }}{1+\left(ı\omega {\tau }_{\epsilon ,2}\right)},\hfill \\ \hfill {\mu }_{\mathrm{r}}& ={\mu }_{\infty }+\frac{{\mu }_{\mathrm{s},1}-{\mu }_{\infty }}{1+\left(ı\omega {\tau }_{\mu ,1}\right)}+\frac{{\mu }_{\mathrm{s},2}-{\mu }_{\infty }}{1+\left(ı\omega {\tau }_{\mu ,2}\right)},\hfill \end{array}$$

where ${\tau }_{\epsilon /\mu ,1/2}$ are relaxation time constants, the subscript “∞” refers to a very high frequency value of the relative material property, the subscript “s” to a static value of the relative material property, and ı denotes the imaginary unit. The waveguide’s input and output boundaries are denoted with ${\mathrm{\Gamma }}_{\mathrm{in}}$ and ${\mathrm{\Gamma }}_{\mathrm{out}}$, respectively. The remaining waveguide walls are assumed to be perfect electric conductors (PECs) and the corresponding boundary is denoted with ${\mathrm{\Gamma }}_{\mathrm{PEC}}$.

We assume that the input port boundary ${\mathrm{\Gamma }}_{\mathrm{in}}$ is excited by an incoming plane wave. We further assume that the excitation coincides with the fundamental transverse electric (TE) mode and that all other propagation modes can be neglected, e.g., are quickly attenuated. Under these assumptions, the electric field $\mathbf{E}$ can be computed by solving the boundary value problem:

$$\begin{array}{}\mathrm{(19a)}& \hfill \nabla \times \left({\mu }^{-1}\nabla \times \mathbf{E}\right)-{\omega }^2\epsilon \mathbf{E}& =0,\hfill & \hfill & \mathrm{in}D,\hfill \mathrm{(19b)}& \hfill {\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{PEC}}}\times \mathbf{E}& =0,\hfill & \hfill & \mathrm{on}{\mathrm{\Gamma }}_{\mathrm{PEC}},\hfill \mathrm{(19c)}& \hfill {\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{in}}}\times \left(\nabla \times \mathbf{E}\right)+ık_z{\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{in}}}\times \left({\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{in}}}\times \mathbf{E}\right)& ={\mathbf{U}}^{\mathrm{i}},\hfill & \hfill & \mathrm{on}{\mathrm{\Gamma }}_{\mathrm{in}},\hfill \mathrm{(19d)}& \hfill {\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{out}}}\times \left(\nabla \times \mathbf{E}\right)+ık_z{\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{out}}}\times \left({\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{out}}}\times \mathbf{E}\right)& =0,\hfill & \hfill & \mathrm{on}{\mathrm{\Gamma }}_{\mathrm{out}},\hfill \end{array}$$

where D is the computational domain, $\omega =2\pi f$ the angular frequency, f the frequency, ${\mathbf{U}}^{\mathrm{i}}$ the incoming plane wave, ${\mathbf{n}}_{{\mathrm{\Gamma }}_{\mathrm{in}}/{\mathrm{\Gamma }}_{\mathrm{out}}/{\mathrm{\Gamma }}_{\mathrm{PEC}}}$ are outwards-pointing normal vectors, and $\mathbf{k}=\left(0,0,k_z\right)$ the wavevector. The QoI is the reflection coefficient at the input port ${\mathrm{\Gamma }}_{\mathrm{in}}$, $r=\left|\frac{{\mathbf{E}}_{{\mathrm{\Gamma }}_{\mathrm{in}}^{-}}}{{\mathbf{E}}_{{\mathrm{\Gamma }}_{\mathrm{in}}^{+}}}\right|\in \left[0,1\right]$. Usually, Equation (19) is solved numerically, e.g., using the finite element method (FEM). For this simple model, a semi-analytical solution exists for the reflection coefficient r. Therefore, errors due to spatial discretization can be neglected.

Equation (19) is transformed into a stochastic problem by recasting all geometrical and Debye material model parameters, as well as the frequency, to RVs. In the nominal configuration, the parameter values are $\overline{f}=6\mathrm{G}\mathrm{Hz}$, $\overline{w}=30\mathrm{m}\mathrm{m}$, $\overline{h}=3\mathrm{m}\mathrm{m}$, $\overline{\ell }=7\mathrm{m}\mathrm{m}$, $\overline{d}=1\mathrm{m}\mathrm{m}$, ${\overline{\epsilon }}_{\mathrm{s},1}=2$, ${\overline{\epsilon }}_{\mathrm{s},2}=2.2$, ${\overline{\epsilon }}_{\infty }=1$, ${\overline{\mu }}_{\mathrm{s},1}=2$, ${\overline{\mu }}_{\mathrm{s},2}=3$, ${\overline{\mu }}_{\infty }=1$, ${\overline{\tau }}_{\epsilon ,1}=1$, ${\overline{\tau }}_{\epsilon ,2}=1.1$, ${\overline{\tau }}_{\mu ,1}=1$, and ${\overline{\tau }}_{\mu ,2}=2$. Each parameter is now assumed to follow a truncated normal distribution, such that $Y_n\sim \mathcal{TN}\left({\overline{y}}_n,{\left(0.01{\overline{y}}_n\right)}^2,0.95{\overline{y}}_n,1.05{\overline{y}}_n\right)$, $n=1,\cdots ,15$, i.e., the standard deviation of the Gaussian distribution is equal to $1\%$ of the nominal parameter value and a maximum deviation of $\pm 5\%$ from the nominal parameter value is imposed.

The comparison between the dimension-adaptive Leja interpolation and the degree-adaptive LAR-gPC for both error metrics in Equations (14) and (15) is presented in Figure 6. Similarly to the results obtained in Section 3.3 and Section 3.5, the Leja interpolation method was significantly superior to the LAR-gPC approximation. For both error metrics, the difference was typically one order of magnitude, with an increasing tendency for larger simulation budgets.

4. Discussion

In this work, we investigated the use of a sparse interpolation method based on weighted Leja node sequences in order to address the problems of approximation and uncertainty quantification for models with input random variables characterized by PDFs of arbitrary shapes. For that purpose, we applied a weighted Leja-based, dimension-adaptive interpolation algorithm to four models featuring 8, 10, 16, and 15 parameters, respectively. Truncated normal and extreme value input distributions have been used to model the random input parameters. The suggested approach has also been compared to a well-established adaptive gPC method, which uses numerically constructed orthogonal polynomials to address arbitrary input PDFs.

The numerical results presented in Section 3 indicated that the weighted Leja interpolation performed very well in all four numerical experiments, both in terms of approximation as well as of statistics estimation accuracy. Moreover, the suggested approach showed either superior or equivalent performance to the gPC method used for comparison purposes. The comparison results could be seen as an extension of those reported in [27,28], where uniform, normal, and log-normal distributions were considered, to the case of arbitrary input PDFs. Overall, we may conclude that the approach presented in this work posed a viable and reliable alternative to arbitrary gPC methods [24,25,26].

An obvious limitation of this work is the assumption that the input RVs could be described by continuous PDFs. Data-driven approaches, such as those presented in [31,32], are able to construct gPC approximations based only on moment values computed after a dataset, without any other assumption regarding a PDF, which might rely on subjective interpretations of the data and result in the introduction of severe biases. The extension of interpolation-based UQ methods to problems where the input PDFs are not explicitly given in closed-form would be an interesting path for further research, to be pursued in a later study.