Approximation and Uncertainty Quantification of Stochastic Systems with Arbitrary Input Distributions using Weighted Leja Interpolation

Approximation and uncertainty quantification methods based on Lagrange interpolation are typically abandoned in cases where the probability distributions of a stochastic system’s input parameters are not normal, uniform, or closely related ones, due to the lack of suitable interpolation nodes. This paper suggests the use of weighted Leja node sequences as a remedy to this situation. We present the recently introduced weighted Leja interpolation rules, 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 the results of three numerical experiments, where the respective models feature extreme value and truncated normal input distributions. The suggested approach is also compared against a well-established polynomial chaos method and found to be either comparable or superior in terms of approximation and statistical moment estimation accuracy. keywords– adaptive algorithms, arbitrary probability distributions, sparse interpolation, uncertainty quantification, weighted Leja sequences.


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 non-intrusive UQ methods, such that the complex and often proprietary simulation softwares which model physical systems can be incorporated into UQ studies without any modifications.
Sampling methods such as (quasi-) Monte Carlo (MC) [12] or latin hypercube sampling (LHS) [28] continue to be the most common approaches in the context of non-intrusive UQ, and for good reasons; they are very simple to implement, embarrassingly 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 sampling methods unattractive or even intractable in cases where high accuracies are sought, or if single simulations are highly time-demanding.
Computationally efficient alternatives have been found since more than 25 years in the use of spectral UQ methods [24,48], which approximate the functional dependence between a stochastic system's random inputs and outputs. The most popular black-box methods of this category employ Lagrange or generalized polynomial chaos (gPC) [49] global polynomial approximations, where the former rely on interpolation schemes [3,4] and the latter are most commonly based on least squares (LS) regression [8,9,32] or pseudo-spectral projection [25] methods. Assuming a smooth inputoutput dependence, spectral UQ methods provide fast convergence rates, in some cases even exponential ones. However, their performance scales badly with the number of input parameters, an effect known as the curse of dimensionality [7]. To overcome this bottleneck, state-of-the-art approaches rely on algorithms for the adaptive construction of sparse approximations [8,9,14,32,36,42].
In the case of Lagrange interpolation methods, adaptive sparse approximation algorithms rely upon the existence of nested sequences of univariate interpolation nodes. Such nested node sequences are readily available in the case of input parameters characterized by uniform [15] or normal [20] probability density functions (PDFs), but not for more general probability measures. This shortcoming is also obvious from the available implementations of several free softwares, e.g., the Sparse-Grid-Matlab-Kit [44], Spinterp [22], Tasmanian [41], or Dakota [1]. On the other hand, the main requirement in gPC approximations, either full or sparse, are polynomials which are orthogonal with respect to the input PDFs. In the case of arbitrary PDFs, such polynomials can be numerically constructed [39,46,47] and corresponding implementations are available in a number of free softwares, e.g., UQLab [31], Chaospy [19], OpenTURNS [5] and Dakota [1]. Unsurprisingly and despite the fact that, according to a number of studies [17,18,34], Lagrange interpolation methods enjoy an advantage over gPC-based ones with respect to the error-cost ratio, they are abandoned whenever the input PDFs are not uniform, normal, or closely related ones, e.g., log-uniform or log-normal. As a natural consequence, the ongoing research work on testing and establishing the suitability of various candidate methods in different settings lacks comparisons between Lagrange interpolations and gPC-based approximations for more "exotic" probability measures.
In this work, we suggest using the recently introduced weighted Leja node sequences [33] as a way to enable the use of sparse interpolation methods for UQ studies in the case of arbitrary input PDFs. We note that we consider here continuous PDFs of arbitrary shapes and not data-driven approaches which only make use of moments without further assumptions regarding the PDFs [2,37]. 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 [14,35,38]. Nonetheless, with the exception of [30] where skewed beta distributions are considered, the use of Leja nodes has been limited to uniform and normal input distributions. The contribution of this work is exactly on the application of weighted Leja interpolation to enablepossibly high-dimensional -approximation and UQ for systems with more arbitrary input PDFs.
To that end, the present paper presents 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. In Section 2 we first present a general overview of the UQ setting in which the suggested approach shall be applied. In Section 3 we recall the basics of univariate Lagrange interpolation, both standard and hierarchical. Original and weighted Leja rules along with their properties are presented in Section 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 in Section 5. Numerical experiments which verify the suitability of the suggested approach for physical systems with arbitrarily distributed uncertain input parameters are made available in Section 6. In the same section, the adaptive weighted Leja interpolation is compared against a well-known adaptive gPC algorithm based on least angle regression (LAR) and numerically constructed orthogonal polynomials. Finally, a discussion on the findings of this work and on possibilities for further research can be found in Section 7.

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 where D = D(u, y) is a differential operator, y ∈ R N an N -dimensional parameter vector, and u = u (y) the parameter-dependent solution of (1). We are interested in a specific model output, commonly called a quantity of interest (QoI), which is computed by post-processing the solution u (y). For simplicity, we assume that the QoI is a real scalar, however, the methodology discussed in Section 5 can be applied to complex or vector-valued QoIs with only minor modifications. We further assume that the functional dependence between the input parameter vector and the QoI is given by the map g : y → g (y), which is assumed to be deterministic, i.e., the exact same output g (y) is observed each and every time the QoI is evaluated for the same parameter vector y. Parametric problems in the form of (1) arise in both deterministic and stochastic settings. An example of the former case would be an optimization study, while an example of the latter case would be a UQ study as the ones considered in this work. In the stochastic setting, the parameter vector y corresponds to a realization of an N -dimensional random variable (RV) Y = (Y 1 , Y 2 , . . . , Y N ), alternatively called a random vector. We introduce a complete probability space (Θ, Σ, P ), where Θ is the sample (outcome) space, Σ the σ-algebra (set of events) and P the probability measure, i.e., a map P : Σ → [0, 1] from events to probabilities. The random vector is defined as the map Y : Θ → Ξ from the sample space to a measurable space Ξ ⊆ R N , called the image space. Accordingly, the realizations y = Y (θ), θ ∈ Θ, take values in Ξ. Finally, the random vector Y is characterized by a PDF ̺ : Ξ → R ≥0 . The notation is summarized in Table 1. Further assuming that the RVs Y n , n = 1, . . . , N , are mutually independent, it holds that While the case of dependent RVs is not addressed in this work, we mention the option of using generalized Nataf transformations [26] to transform the dependent RVs into independent ones. The sparse interpolation method presented in Section 5 can then be applied to the transformed RVs. Due to its dependence on random 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 as well. The main goal of UQ studies is to quantify this output uncertainty by computing statistical measures, such as moments, event probabilities, or sensitivity metrics. Denoting the expectation operator with E [·], any statistical measure can be written in the general form where the functional φ is chosen according to the sought statistical measure, e.g., φ (g) = g corresponds to the expected value of the QoI, φ (g) = (g − E [g]) 2 to its variance, and so on. A secondary Table 1: Notation summary regarding stochastic parametric models.
parameter-dependent PDE solution g : R N → R map from the input parameters to the QoI g (y) parameter-dependent QoI (Θ, Σ, P ) complete probability space Θ sample/outcome space Σ σ-algebra/set of events probability density function (PDF) goal is to approximate the input-output map g with an inexpensive surrogate model g ≈ g, which could replace the original model in computationally expensive tasks, e.g., sampling-based UQ, optimization under uncertainty, or others. Both of these goals are pursued here by means of the sparse interpolation method discussed in Section 5.

Univariate Interpolation Schemes
Let us assume a univariate input-output map g (y), where y = Y (θ) denotes a realization of a single RV Y . We let the non-negative integer i ∈ Z ≥0 denote the interpolation level, j=1 the level-i grid of interpolation nodes, and m : Z ≥0 → Z + , m (0) = 1, a monotonically increasing "levelto-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.

Univariate Lagrange Interpolation
The level-i univariate Lagrange interpolation reads where l i,j are univariate Lagrange polynomials, defined on the interpolation grid Z i as Given a smooth map g, the interpolation coefficients are expected to decay fast [3,4]. Hence, the magnitudes of those coefficients can be used to appropriately truncate the expansion (5), e.g., omit terms corresponding to coefficients with magnitudes below a certain tolerance. In that way, the interpolation can be truncated according to the desired accuracy. 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 I i [g] is related to the best polynomial approximation error ǫ best through [4,14] g where P i denotes the space of Lagrange polynomials with degree up to m (i) − 1 and L i the Lebesgue constant 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.

Hierarchical Univariate Interpolation
We further demand that the interpolation nodes form sequences of nested grids for increasing interpolation levels, such that In that case, we may define a hierarchical interpolation scheme as follows. We first define the operator ∆ i as the difference between two consecutive univariate interpolants, i.e., where I −1 is the null operator. The interpolation (5) can now be given as a sum of interpolant differences, such that The hierarchy of formula (10) is obvious if we consider two consecutive interpolation levels, i − 1 and i, in which case the level-i interpolation reads where the interpolation coefficients s i,j , given by are called hierarchical surpluses. The obvious advantage of using nested grids and the hierarchical formulas (10) and (11) is that at each new level i the QoI must be evaluated only at the new nodes Moreover, univariate hierarchical interpolation rules constitute the backbone of the sparse adaptive multivariate interpolation algorithms discussed in Section 5.

Interpolatory Univariate Quadrature
Considering the UQ goal of estimating specific statistical measures given by (4), we may use an available interpolation I i [g] to derive an appropriate quadrature scheme. For example, assuming a univariate interpolation given by the expected value of the QoI can be estimated as where w i,j are quadrature weights which can be analytically computed by Quadrature schemes for the estimation of statistical measures other than the expected value can be derived in a similar way, by applying the expectation operator along with the functional φ which corresponds to the sought statistical measure, e.g., analytical formulas for higher order moments can be found in [1].

Leja Interpolation Nodes
As already pointed out in Section 3, 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 to the interpolation nodes must increase at a rate lower than the polynomial approximation error decreases [4,14]. Second, the choice of nodes must result in accurate quadrature schemes, similar to (14), to be used in the computation of statistical measures. 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 will be shown in Section 5, 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. The first three points are easily addressed in the case of a uniform or normal PDF by choosing Clenshaw-Curtis [15] and Genz-Keister [20] nodes, respectively. However, satisfying the final requirement poses significant difficulties. In principle, one could derive specific interpolation and quadrature rules for each given PDF [10,40]. However, this is a cumbersome task, for which dedicated analyses are necessary each and every time a new PDF is considered. This can be avoided by using weighted Leja sequences [33], discussed in the following.

Unweighted Leja Nodes
A sequence of standard, unweighted Leja nodes [27] {y j } j≥0 , y j ∈ [−1, 1], is defined by solving the optimization problem for each node y j , j ≥ 1. The initial node y 0 ∈ [−1, 1] can be chosen arbitrarily, therefore Leja sequences are not unique. Using simple scaling rules, the unweighted Leja sequences defined in (16) 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 4.2.
With respect to the remaining node requirements, it has been shown that the Lebesgue constant associated with unweighted Leja nodes increases subexponentially [45], where exponential convergence is expected for sufficiently smooth QoIs. Regarding Leja-based quadrature rules, an extensive discussion can be found in [33]. Considering uniform PDFs, a number of works [30,33,38] show that Leja-based quadrature schemes are sufficiently accurate, albeit suboptimal compared to more standard choices, e.g., Clenshaw-Curtis-based ones. Finally, the optimization problem (16) results in nested node sequences irrespective of the employed level-to-nodes function m (i). 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 (i) = i + 1 to get the minimum of one extra node per interpolation level, i.e., # (Z i \ Z i−1 ) = 1. Another popular choice is to use m(i) = 2i + 1 [38], which is typically employed for the construction of symmetric Leja sequences, defined as with y 0 = 0 and y 1 = 1, as well as for R-Leja sequences [13], defined as projections on [−1, 1] of Leja sequences computed on the complex unit disc, such that where y 0 = 0, y 1 = 1, and R (z) denotes the real part of a complex number z. In all cases, Leja sequences are more granular compared to other nested node sequences, e.g., the level-to-nodes function m(i) = 2 i + 1 must be used for nested Clenshaw-Curtis nodes.

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 ̺ (y), which acts as a weight function. To that end, we employ the definition of weighted Leja sequences given in [33].
Using that formulation, weighted Leja sequences are constructed by incorporating the PDF ̺ (y) into the optimization problem (16), which is transformed into We note that other formulations exist as well [16], however, the one given in (19) is preferred due to the fact that the resulting Leja nodes are asymptotically distributed as weighted Gauss quadrature nodes. While this is a promising property, it is not necessarily sufficient to produce an accurate interpolation rule. In particular, unlike the unweighted case, weighted Leja rules are not associated with a subexponentially growing Lebesgue constant, hence, an improved accuracy cannot be guaranteed for interpolation grids of increasing sizes. Nonetheless, weighted Leja interpolation has been found to perform very well in practice [30], as also demonstrated by the results of the numerical experiments presented in Section 6. The remaining good 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, equivalently, 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 [33].

Sparse Adaptive Leja Interpolation
The multivariate interpolation is based on suitable combinations of the univariate interpolation rules discussed in Section 3. We denote with i = (i 1 , i 2 , . . . , i N ) a multi-index of interpolation levels per input RV, with Λ a multi-index set, and with ∆ i the tensor-product difference operator given by where the univariate difference operators ∆ n,in are defined as in (9), to be recalled in the following subsections. In Section 5.1 we discuss interpolation on generalized sparse grids. A dimensionadaptive interpolation scheme based on Leja sequences is presented in Section 5.2. Post-processing multivariate interpolations for UQ purposes is discussed in Section 5.3.

Generalized Sparse Grid Interpolation
Given a multi-index set Λ, the interpolation-based approximation of the multivariate input-output map g(y) reads while the corresponding grid of multivariate interpolation nodes Z Λ is given by We note that the formula (21) is not in general interpolatory. The interpolation property holds only if Λ is a downward-closed set and the underlying univariate interpolation rules are based on nested nodes [4]. Downward-closed sets, also known as monotone or lower sets, satisfy the property ∀i ∈ Λ ⇒ i − e n ∈ Λ, ∀n = 1, 2, . . . , N, with i n > 0, where e n denotes the unit vector in the nth dimension. If, additionally, the multi-index set Λ satisfies some sparsity constraint, then (21) is a sparse interpolation formula and the grid in (22) is called a sparse grid. A commonly used sparsity constraint is resulting in the so-called isotropic sparse grids [3,4]. Compared to full isotropic tensor grids, i.e., those given by the complexity of which is O m (i max ) N , isotropic sparse grids delay significantly the curse of dimensionality, since their complexity is reduced to O m (i max ) (log m (i max )) N −1 . Nonetheless, isotropic sparse grids still grow exponentially fast with respect to the number of parameters.
In most practical cases, the influence of the input parameters upon the QoI is anisotropic, e.g., some parameters have greater impact or their dependence is more difficult to approximate. Hence, the interpolation (21) can be constructed such that this anisotropy is represented in its terms, equivalently in the multi-indices comprising the set Λ [14,36]. This anisotropic refinement of the parameter space according to the contribution of each parameter to the interpolation results in greater interpolation accuracy for the same number of terms, equivalently, for equal costs, compared to isotropic schemes. Since, in most cases, the parameter anisotropy cannot be a priori estimated in an accurate way [6], anisotropic interpolations are usually constructed using greedy, adaptive algorithms [21,22,42]. For the particular case of Leja nodes, such an algorithm is discussed in Section 5.2.

Adaptive Anisotropic Leja Interpolation
We will base the adaptive construction of anisotropic sparse interpolations on the dimension-adaptive algorithm first presented in [21] for quadrature purposes. Variants of this algorithm appear in a number of later works [14,22,30,33,38]. In this work, we assume the all underlying univariate interpolation rules employ Leja nodes, presented in Section 4, along with the level-to-nodes function m(i) = i + 1. Our approach is depicted in Algorithm 1, while a detailed presentation follows.
Due to the fact that all univariate Leja rules employ the level-to-nodes function m(i) = i + 1, each multi-index i ∈ Λ adm is associated with a single new interpolation node y i = Z i \ Z Λ . Therefore, should an admissible multi-index i ∈ Λ adm be added to Λ, we obtain the hierarchical interpolation scheme where, similarly to (12), the hierarchical surpluses s i are given by and L i are multivariate Lagrange polynomials given by where, due to the use of nested Leja sequences with the level-to-nodes function m(i) = i + 1, the univariate Lagrange polynomials can now be defined hierarchically [14], such that As already discussed in Section 3.1 for the univariate case, interpolation 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 i ∈ Λ adm the contribution indicator Naturally, the multi-index set Λ should be enriched with the admissible multi-index i.e., the one 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 Λ∪Λ adm . After the termination of the iterations, the final interpolation is constructed using all multi-indices in Λ ∪ Λ adm , since the hierarchical surpluses corresponding to the admissible nodes have already been computed.

Post-Processing
Once available, the sparse interpolation can be used as an inexpensive surrogate model which can replace the original model in computationally demanding tasks. For example, any statistical measure E [φ (g)] can be estimated by sampling the interpolation I Λ [g] instead of g. While the slow convergence of sampling methods remains, the interpolation is typically computationally inexpensive compared to the full simulation model, thus enabling the use of much larger sample sizes.
Specific statistical measures can be computed directly after the interpolation terms. One possibility is to first transform the Lagrange basis into a gPC basis [11] 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 sought statistical measure directly upon the interpolation, as to derive an appropriate quadrature scheme, similar to the univariate case presented in Section 3.3. For example, assuming a Leja-based interpolation given by the expected value of the QoI can be estimated as where the quadrature weights w i are given as products of the respective univariate ones, i.e., Analytical formulas similar to (34) can be derived for other statistical measures as well, e.g., for the QoI's variance, skewness, and kurtosis [1].

Numerical Experiments
In this section, we consider stochastic 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. 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. 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 [39]. 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 [29]. 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 [31] 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 [39].

Error Metrics
Let us assume that an interpolation or gPC-based polynomial approximation g ≈ g is available. The following errors are used to estimate the performance of g for approximation and UQ purposes.
First, we want to estimate the performance of 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 cross-validation sample {y q } Q q=1 which is randomly drawn from the joint input PDF ̺(y) and compute the root-mean-square (RMS) crossvalidation error From an interpretation point of view, ǫ cv,RMS indicates the average (expected) approximation accuracy of g, while at the same time being sensitive to outliers by penalizing large errors. In all numerical experiments, the size of the cross-validation sample is set to Q = 10 5 . Secondly, we want to estimate the performance of g in terms of UQ quality, i.e., estimate the accuracy of the statistical measure values provided by post-processing the approximation. For that purpose, we use the expected value of the QoI as a representative statistical measure and compute the relative error where E [g] is a reference value and E [g] is an estimate computed by directly post-processing the terms of 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 .

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 only analytical 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, . . . , 100, 200, . . . , 1000, and compare the errors (36) and (37) 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 not being addressed here, those differences must be taken into account by practitioners, considering the problem and the computational resources at hand.

Borehole Model
We consider the 8-dimensional parametric function which models the water flow through a borehole [43]. The water flow is given in m 3 /yr and the parameter vector is y = (r w , r, T u , H u , T l , H l , L, K w ). All model parameters are presented in Table 2.
We follow the setting given in [43], 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 [10]. We denote a truncated normal distribution with T N µ, σ 2 , l, u , where µ and σ 2 refer to the mean value and the variance of a normal distribution, while l and u are the lower and upper truncation limits, i.e., the minimum and maximum values of the image space Ξ. In this numerical experiment, the input parameters are modeled as follows.
• The parameter r w originally follows the normal distribution N µ = 0.1, σ = 0.0161812 2 . The distribution is now truncated to the value range [0.05, 0.15].
Therefore, the parameter r has a mean value equal to exp µ LN + • The remaining parameters are originally uniformly distributed. Assuming a uniform distribution with support in [a, b], the corresponding truncated normal distribution is given as T N µ = a+b 2 , σ 2 = (b−a) 2 12 , l = a, u = b , 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 probability distribution of each parameter is given in Table 3.
The borehole model (38) is approximated by both the dimension-adaptive weighted Leja interpolation and the degree-adaptive LAR-gPC expansion. Figure 1 shows the error-cost relation for both approximations, with respect to the errors (36) and (37). Regarding the ǫ cv,RMS error, Leja interpolation outperforms the LAR-gPC method by approximately one order of magnitude over the  , σ 2 = e 1.0056 2 − 1 e 2·7.71+1.0056 2 , l = 100, u = 50000 whole cost-range. Leja interpolation is also advantageous with respect to the ǫ rel,E error, again always outperforming the LAR-gPC method and stagnating much sooner. We note that this error stagnation is observed due to the fact that both methods reach the accuracy of the quasi-MC estimate. Therefore, for the case of the borehole model, the weighted Leja interpolation approach is not only found to perform well for the considered truncated normal input PDFs, but is also superior to the LAR-gPC method for both considered error metrics.   [29,30] and with a degree-adaptive LAR-gPC algorithm [9,31]. The size of the cross-validation sample is Q = 10 5 . The reference expected value is computed with quasi-MC integration based on a Sobol sample of size equal to 10 8 .

Steel Column Limit State Function
We consider the 10-dimensional parametric limit state function which relates the reliability of a steel column to its cost [23,43]. The parameter vector is y = (F, P d , P 1 , P 2 , B, D, H, F 0 , E), where P t = P d + P 1 + P 2 is the total load and E b = π 2 EBDH 2 2L 2 is known as the Euler buckling load. All parameters of function (39) are presented in Table 4.
Once more, we follow the original setting given in [43,23]. Originally, P 1 , P 2 , and E follow extreme value distributions, where P 1 and P 2 are modeled with Gumbel distributions and E with a Weibull distribution, F s , B, D and H follow log-normal distributions, and P d and F 0 follow normal distributions. In this numerical experiment, we transform all log-normal and normal distributions to truncated normal ones and use the Gumbel distribution to model the remaining three parameters. The corresponding distributions are given in Table 5, where the truncated normal distribution is denoted as in Section 6.3 and the Gumbel distribution is denoted with G (ℓ, β), where ℓ is the location parameter and β the scaling parameter. With the exception of the parameter L, for which the truncation limits are set to µ ± 4σ, all other truncated normal distributions have support in [µ − 3σ, µ + 3σ]. For all parameters which follow a truncated normal distribution, µ and σ 2 coincide with the original mean and variance values given in [43,23]. Accordingly, the parameters ℓ and β of the Gumbel distributions are chosen such that the original mean values and variances given in [43,23] are preserved.
We approximate the steel column function (39) using the dimension-adaptive weighted Leja interpolation and the degree-adaptive, LAR-based gPC method. For both approximations and for both error metrics (36) and (37), the error-cost relation is presented in Figure 2. Once more, the Leja interpolation method performs very well for the given input PDFs. However, contrary to the results of Section 6.3, no advantage is observed over the LAR-gPC method. With respect to the ǫ cv,RMS error, the performance of both methods is comparable. 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 ǫ rel,E error. Overall, the weighted Leja interpolation is again able to yield accurate approximations and statistical measures and it may be regarded as competitive to the LAR-gPC method. Table 5: Input distributions of the steel column function.

Meromorphic Function
We consider the 16-dimensional meromorphic function where w is a vector of positive weights which govern the parameter anisotropy of function (40) [32]. The weight vector is given by w =ŵ 2 ŵ 1 , whereŵ = 1, 0.   [29,30] and with a degree-adaptive LAR-gPC algorithm [9,31]. The size of the cross-validation sample is Q = 10 5 . The reference expected value is computed with quasi-MC integration based on a Sobol sample of size equal to 10 8 .
While a normal distribution truncated exactly at its 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 3 shows the error-cost relations corresponding to the error metrics (36) and (37) for both approximations of the meromorphic function (40), i.e., using the dimension-adaptive weighted Leja interpolation algorithm and the degree-adaptive LAR-gPC algorithm. The results resemble those of Section 6.3, i.e., the Leja interpolation method has a clear advantage over the LAR-gPC approximation with respect to both error metrics. Regarding the ǫ cv,RMS error, the difference between the two approximations is typically greater than one order of magnitude, while it increases for an increasing computational budget. Regarding the ǫ rel,E error, the two methods show comparable performance for costs up to 100 simulations, however, the Leja interpolation is again clearly superior as the simulation budget increases. Hence, in this numerical experiment, the weighted Leja interpolation is able to accurately approximate a model with very atypical input PDFs, and provide accurate estimations of statistical measures. Moreover, it is found to have a clear edge over to the well-established LAR-gPC method.

Summary and Conclusions
In this work, we have suggested and 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 stochastic models with random variables characterized by PDFs of arbitrary shapes. For that purpose, we applied a weighted Leja-based dimension-adaptive interpolation algorithm to three stochastic models featuring 8, 10, and 16 input 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 6 indicate that the weighted Leja interpolation performs very well in all three numerical experiments, both in terms of approximation as well as of statistical measure estimation accuracy. Moreover, the suggested approach shows either superior or equivalent performance to the gPC method used for comparison purposes. Those comparison results can be seen as an extension of those reported in [17,18], 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 poses a viable and reliable alternative to arbitrary gPC methods [39,46,47].
An obvious limitation of this work is the assumption that the input RVs can be described by continuous PDFs. Data-driven approaches, such as those presented in [2,37], 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.