K Konakli

B Sudret

Reliability analysis of high-dimensional models using low-rank tensor approximations

Keywords:

Engineering and applied sciences use models of increasing complexity to simulate the behavior of manufactured and physical systems. Propagation of uncertainties from the input to a response quantity of interest through such models may become intractable in cases when a single simulation is time demanding. Particularly challenging is the reliability analysis of systems represented by computationally costly models, because of the large number of model evaluations that are typically required to estimate small probabilities of failure. In this paper, we demonstrate the potential of a newly emerged meta-modeling technique known as low-rank tensor approximations to address this limitation. This technique is especially promising for high-dimensional problems because: (i) the number of unknowns in the generic functional form of the meta-model grows only linearly with the input dimension and (ii) such approximations can be constructed by relying on a series of minimization problems of small size independent of the input dimension. In example applications involving finite-element models pertinent to structural mechanics and heat conduction, low-rank tensor approximations built with polynomial bases are found to outperform the popular sparse polynomial chaos expansions in the estimation of tail probabilities when small experimental designs are used. It should be emphasized that contrary to methods particularly targeted to reliability analysis, the meta-modeling approach also provides a full probabilistic description of the model response, which can be used to estimate any statistical measure of interest.

INTRODUCTION

Analysis of the response of engineered and physical systems under uncertainties is of key importance in risk assessment and decision making in a wide range of fields. To this end, it is typical to use a computer model to represent the behavior of a system and perform repeated simulations to propagate uncertainties from the input to a response quantity of interest. However, because of the growing complexity of the computer models used across engineering and sciences, in many practical situations, a single simulation is time consuming, thus rendering uncertainty propagation non-affordable. Such situations are often encountered in reliability analysis due to the large number of model evaluations required to compute small failure probabilities. As a result, meta-modeling techniques are gaining increasing popularity. The key idea thereof is to substitute a computationally expensive model with a statistically equivalent one, so-called metamodel, which can be easily evaluated. Using the meta-model, the analyst can perform statistical analysis of a response quantity of interest at low cost.

Of interest herein is non-intrusive meta-modeling, in which the original model is treated as a "black box". Building a meta-model in an non-intrusive manner relies on the evaluation of the original model at a set of points in the input space, called experimental design. The efficiency of a meta-modeling technique depends on its ability to provide sufficiently accurate representations of the exact model response over the entire input range by using relatively small experimental designs. This can be particularly challenging in cases when determining the tails of the response distribution with high accuracy is important, as in the estimation of small exceedence probabilities required in reliability analysis.

In this paper, we demonstrate the potential of the newly-emerged approach called low-rank tensor approximations to provide meta-models appropriate for reliability analysis. Although different types of tensor decompositions may be used (see e.g. [START_REF] Kolda | Tensor decompositions and applications[END_REF]; [START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF]; [START_REF] Grasedyck | A literature survey of low-rank tensor approximation techniques[END_REF]), we confine our attention to canonical tensor formats.

In this context, low-rank approximations express the model response as a sum of a small number or rank-one tensors, where a rank-one tensor is a product of univariate functions in each of the input parameters. The idea of such decompositions originates in the work of Hitchcock [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF] in the first half of the 20th century and has been employed within the last 50 years in a wide range of fields, including -but not limited to -psychometrics [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "Eckart-Young" decomposition[END_REF]; [START_REF] Harshman | Foundations of the PARAFAC procedure: Models and conditions for an "explanatory" multi-modal factor analysis[END_REF], chemometrics [START_REF] Appellof | Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents[END_REF]; [START_REF] Bro | PARAFAC. tutorial and applications[END_REF], neuroscience [START_REF] Mocks | Topographic components model for event-related potentials and some biophysical considerations[END_REF]; [START_REF] Andersen | Structure-seeking multilinear methods for the analysis of fMRI data[END_REF], fluid mechanics [START_REF] Felippa | Mixed variational formulation of finite element analysis of acoustoelastic/slosh fluid-structure interaction[END_REF]; [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF], signal processing [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF]; De [START_REF] De Lathauwer | Tensor-based techniques for the blind separation of DS-CDMA signals[END_REF], image analysis [START_REF] Shashua | Linear image coding for regression and classification using the tensor-rank principle[END_REF]; [START_REF] Furukawa | Appearance based object modeling using texture database: Acquisition compression and rendering[END_REF] and data mining [START_REF] Acar | Collective sampling and analysis of high order tensors for chatroom communications[END_REF]; [START_REF] Beylkin | Multivariate regression and machine learning with sums of separable functions[END_REF]. This technique is now attracting an expanding interest from the more recently established community of uncertainty quantification [START_REF] Nouy | Proper generalized decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF]; [START_REF] Khoromskij | Tensor-structured Galerkin approximation of parametric and stochastic elliptic PDEs[END_REF]; Chevreuil et al. (2013,?); [START_REF] Doostan | Non-intrusive low-rank separated approximation of high-dimensional stochastic models[END_REF]; [START_REF] Hadigol | Partitioned treatment of uncertainty in coupled domain problems: A separated representation approach[END_REF]; [START_REF] Rai | Sparse Low Rank Approximation of Multivariate Functions -Applications in uncertainty quantification[END_REF]; [START_REF] Validi | Low-rank separated representation surrogates of high-dimensional stochastic functions: Application in Bayesian inference[END_REF]; Konakli and Sudret (2015b).

The focus of the present study is on low-rank tensor approximations that are made of polynomial functions due to the simplicity and versatility these offer. The considered meta-models therefore constitute an alternative to the widely used polynomial chaos expansions. In the latter, the number of unknown coefficients grows exponentially with the input dimension, requiring prohibitively large experimental designs when high-dimensional models are considered. Conversely, in low-rank tensor approximations, the number of unknown coefficients grows only linearly with the input dimension, which makes the approach particularly promising for dealing with high dimensionality. Such approximations can be constructed through a series of minimization problems of small size that is independent of the input dimension.

The paper is organized as follows: In Section 2, we describe the mathematical setup of nonintrusive meta-modeling and review basic concepts of reliability analysis with meta-models. In Section 3, we present the formulation of canonical low-rank approximations with polynomial bases and detail an algorithm for their construction; in the same section, we present, in a comparative way, the polynomial chaos expansions approach. In Section 4, we confront canonical low-rank approximations to state-of-art polynomial chaos expansions in the estimation of response probability density functions and of failure probabilities in reliability analysis. The paper concludes with a summary of the main findings and respective outlooks.

META-MODELS FOR UNCERTAINTY PROPAGATION AND RELIABILITY ANALYSIS

Non-intrusive meta-modeling and error estimation

We consider a computational model M that represents the behavior of a physical or engineering system of interest. We denote by X = {X 1 , . . . , X M } the M -dimensional input vector of M and by Y = {Y 1 , . . . , Y N } the N -dimensional response vector. To account for the uncertainty in the input and the resulting uncertainty in the response, the elements of X and Y are described by random variables. For the sake of simplicity, the case of a scalar model response (N = 1) is considered hereafter. Therefore, M is considered as the mapping:

X ∈ D X ⊂ R M -→ Y = M(X) ∈ R, (1)
where D X denotes the support of X. Note that the case of a vector model response can be addressed by separately treating each element of Y as in the case of a scalar response.

In general, the map in Eq. ( 1) is not known in a closed analytical form and may represent a computationally intensive process. A meta-model M is an analytical function that mimics the behavior of M; in other words, M(X) possesses similar statistical properties with M(X), while maintaining an easy-to-evaluate form. Replacing a complex computational model with a metamodel allows efficient uncertainty propagation from the random input to a response quantity of interest in cases when this is practically infeasible by using the original model due to the high computational cost.

In non-intrusive meta-modeling, which is of interest herein, the original computational model is treated as a "black box". Thus, in order to develop a meta-model in a non-intrusive manner, one only needs a set of N realizations of the input vector E = {χ (1) , . . . , χ (N ) }, called experimental design (ED), and the corresponding set of model responses Y = {M(χ (1) ), . . . , M(χ (N ) )}.

We underline that non-intrusive approaches do not require any further knowledge of the original model, which is used therein without any modification.

To define measures of accuracy of the meta-model response Y = M(X), we first introduce the discrete L 2 semi-norm:

a X = 1 n n i=1 a 2 (x i ) 1/2 , (2)
where a represents a function:

x ∈ D X -→ a(x) ∈ R and X = {x 1 , . . . , x n } ⊂ D X denotes
a set of n realizations of X. A good measure of accuracy of the meta-model response is the generalization error Err G , which is defined as the mean-square error of the difference (Y -Y )

and can be estimated by:

Err G = M -M 2 X val , (3)
where X val = {x 1 , . . . , x n val } is a sufficiently large set of realizations of X, called validation set.

The estimate of the relative generalization error err G , is obtained by normalizing Err G with the empirical variance of Y val = {M(x 1 ), . . . , M(x n val )}, which denotes the set of model responses at the validation set. Unfortunately, a validation set is not available in typical meta-modeling applications, where a large number of model evaluations is not affordable. An alternative estimate that relies solely on the ED is the empirical error Err E , which is given by:

Err E = M -M 2 E . (4)
In the above equation, the subscript E emphasizes that the semi-norm is evaluated at the points of the ED. The relative empirical error err E , is obtained by normalizing Err E with the empirical variance of Y = {M(χ (1) ), . . . , M(χ (N ) )}, which denotes the set of model responses at the ED.

Although the empirical error reuses the points of the ED, it has the strong drawback that it tends to underestimate the actual generalization error, which might be severe in cases of overfitting.

By using the information contained in the ED only, one can obtain fair approximations of the generalization error by means of cross-validation (CV) techniques (see e.g. [START_REF] Viana | Multiple surrogates: how cross-validation errors can help us to obtain the best predictor[END_REF]; [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF]). In brief, the basic idea of k-fold CV is to randomly partition the ED into k sets of approximately equal size, build the meta-model by relying on all but one of the partitions, and use the excluded set to evaluate the generalization error. By alternating through the k sets, one obtains k meta-models of which the average generalization error serves as the error estimate of the meta-model built with the full ED.

Reliability analysis

Reliability analysis aims at computing the probability that the system under consideration fails to satisfy prescribed criteria. A failure criterion is mathematically represented by the so-called limit-state function g(Z); in a general case, Z = {X, M(X), X } may depend on the input parameters of the model describing the system, response quantities obtained from the model and additional random parameters gathered in X . Conventionally, the limit-state function is formulated so that failure corresponds to g(z) ≤ 0; the set of points that satisfy this condition comprise the failure domain D f with respect to g, i.e. D f = {z : g(z) ≤ 0}. The associated probability of failure is therefore given by:

P f = D f f Z (z)dz, (5)
where f Z denotes the probability density function (PDF) of Z. Note that in a general case, failure of a system may be defined in terms of multiple limit-state functions representing different failure criteria.

A universal method for computing the integral in Eq. ( 5) is Monte Carlo simulation (MCS).

The MCS approach involves generating a sufficiently large sample of realizations of Z, say {z 1 , . . . , z n }, and then, estimating P f as the empirical mean:

P f, MCS = 1 n n k=1 I D f (z k ), (6)
where I D f denotes the indicator function of the failure domain. Obviously, the MCS-based estimator in Eq. ( 6) is unbiased, which means E P f, MCS = P f . Typically, of interest in reliability analysis are failure events with small probabilities of occurrence, i.e. P f << 1, leading to the following approximation of the coefficient of variation (CoV) of P f, MCS :

δ P f, MCS ≈ 1/ nP f . (7)
Eq. ( 7) indicates that estimation of a failure probability with magnitude of the order of 10 -k with CoV < 10% requires that a number of samples larger than 10 k+2 is used. Clearly, MCS is impractical for computing small failure probabilities in cases when a single evaluation of the model response is computationally costly. This limitation is overcome when the original model M is substituted by a meta-model M. Accordingly, the actual failure domain D f is approximated by D f = {z : g(z) ≤ 0}, where g(Z) = g(X, M(X), X ). Once the meta-model is available, evaluation of P f, MCS in Eq. ( 6) by using large Monte Carlo samples becomes essentially costless.

Various techniques have been devised with the purpose of computing efficiently the small failure probabilities that are of interest in reliability analysis; a thorough review listing the advantages and drawbacks of different methods can be found in [START_REF] Morio | A survey of rare event simulation methods for static input-output models[END_REF]. We herein briefly describe three widely-used methods, which are considered in the application section of the present study:

• First-order reliability method (FORM) [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF]; [START_REF] Rackwitz | Structural reliability under combined load sequences[END_REF]: FORM relies on determining the design point P * , i.e. the point of the failure domain that is closest to the origin in the standard normal space. The failure domain is then approximated by the half space defined as the hyperplane that is tangent to the limit-state surface {z : g(z) = 0} at P * , leading to the first-order approximation of P f .

• Second-order reliability method (SORM) [START_REF] Breitung | Asymptotic approximations for probability integrals[END_REF]; Der Kiureghian and de Stefano (1991): SORM provides a correction to the FORM solution by approximating the limitstate surface at the design point by a second-order surface.

• Importance sampling (IS) [START_REF] Melchers | Importance sampling in structural systems[END_REF]; [START_REF] Au | Important sampling in high dimensions[END_REF]: The basic idea in IS is to recast the definition of P f in Eq. ( 5) by means of an auxiliary PDF that is more efficient in generating samples within the failure domain; appropriate weights are introduced in the computation of the integral in order to account for the change in the PDF.

We emphasize that the aforementioned methods are particularly targeted to reliability analysis; on the other hand, meta-modeling comprises a more general tool for uncertainty propagation, which may be used to conduct any type of statistical analysis of the model response, e.g. PDF estimation, evaluation of statistical moments and confidence intervals, analysis of variance, and so forth (see e.g. [START_REF] Xiu | Modelling uncertainty in steady state diffusion problems via generalized polynomial chaos[END_REF]; [START_REF] Acharjee | Uncertainty propagation in finite deformations -A spectral stochastic Lagrangian approach[END_REF] 3 LOW-RANK TENSOR APPROXIMATIONS

Formulation using polynomial bases

We consider the map in Eq. ( 1) assuming that the components of X are independent, with the marginal PDF of X i denoted by f Xi for {i = 1, . . . , M }. Let M LRA denote a meta-model of M belonging to the class of low-rank approximations (LRA); as mentioned in the Introduction, the term "rank" herein refers to the so-called "canonical rank". The corresponding approximation of Y = M(X) has the general form:

Y LRA = M LRA (X) = R l=1 b l w l (X), (8)
in which b l is a normalizing constant and w l is a rank-one function of X. The rank-one function w l is a product of univariate functions of the components of X:

w l (X) = M i=1 v (i) l (X i ), (9)
where v

(i) l
denotes a univariate function of X i . Accordingly, R in Eq. ( 8) represents the number of rank-one components retained in the approximation. Naturally, representations with a small number of rank-one components are of interest, thus named low-rank.

In order to obtain a representation of Y = M(X) in terms of polynomial functions, we expand v

(i) l onto a polynomial basis {P (i)
k , k ∈ N} that is orthonormal with respect to f Xi , i.e. satisfies:

< P (i) j , P (i) k > = DX i P (i) j (x i )P (i) k (x i )f Xi (x i )dx i = δ jk , (10)
where D Xi denotes the support of X i and δ jk is the Kronecker delta symbol, equal to one if j = k and zero otherwise. Accordingly, the univariate function of X i takes the form:

v (i) l (X i ) = pi k=0 z (i) k,l P (i) k (X i ), (11)
where

P (i) k
is the k-th degree univariate polynomial in the i-th input variable of maximum degree p i and z

(i) k,l is the coefficient of P (i)
k in the l-th rank-one term. By substituting Eq. (11) into Eq. (8), we obtain:

Y LRA = M LRA (X) = R l=1 b l M i=1 pi k=0 z (i) k,l P (i) k (X i ) . ( 12
)
Disregarding the redundant parameterization arising from the normalizing constants, the number of unknowns in Eq. ( 12) is R • M i=1 (p i + 1), which grows only linearly with the input dimension M . We will see later that this is a key factor for the higher efficiency of LRA as compared to polynomial chaos expansions when dealing with high-dimensional problems.

Classical algebra allows one to build a family of polynomials satisfying Eq. ( 10) [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]. For standard distributions, the associated families of orthonormal polynomials are well known; for instance, a uniform variable with support [-1, 1] is associated with the family of Legendre polynomials, whereas a standard normal variable is associated with the family of Hermite polynomials [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. However, it is common in practical situations that the input variables do not follow standard distributions. In such cases, the random vector X is first transformed into a basic random vector U (e.g. a standard normal or standard uniform random vector) through an isoprobabilistic transformation X = T -1 (U ) and then, the model response M(T -1 (U )) is expanded onto the polynomial basis associated with U .

Cases with mutually dependent input variables can also be treated through an isoprobabilistic transformation into a vector of independent variables, e.g. the Nataf transformation in the case of a joint PDF with Gaussian copula Lebrun and Dutfoy (2009b,a). We underline that although the focus of the present work is on LRA developed with polynomial functions, the use of such functions is not a constraint in a general case.

Construction with greedy approaches

Different non-intrusive algorithms have been proposed recently for developing LRA in the form of Eq. ( 12); see e.g. Chevreuil et al. (2013,?); [START_REF] Doostan | Non-intrusive low-rank separated approximation of high-dimensional stochastic models[END_REF]; [START_REF] Rai | Sparse Low Rank Approximation of Multivariate Functions -Applications in uncertainty quantification[END_REF]; [START_REF] Validi | Low-rank separated representation surrogates of high-dimensional stochastic functions: Application in Bayesian inference[END_REF].

A common attribute of these algorithms is that the computation of the polynomial coefficients relies on an alternated least-squares (ALS) minimization approach. The ALS technique consists in solving a series of small-size least-squares minimization problems, where each minimization is performed along a single dimension. Chevreuil et al. Chevreuil et al. (2013) proposed to construct LRA in a greedy manner by successively adding rank-one components and updating the entire set of normalizing constants following each increase of the rank. Aspects of this algorithm were further investigated by Konakli and Sudret Konakli and Sudret (2015a). This greedy approach is employed in the present study and described analytically below.

Let Y LRA r = M LRA r (X) denote the rank-r approximation of Y = M(X): Y LRA r = M LRA r (X) = r l=1 b l w l (X), (13)
where:

w l (X) = M i=1 pi k=0 z (i) k,l P (i) k (X i ) . ( 14
)
The employed algorithm comprises a sequence of pairs of a correction step and an updating step, so that the r-th correction step yields the rank-one component w r and the r-th updating step yields the set of coefficients {b 1 , . . . , b r }. Details on these steps are given next.

Correction step: Let R r (X) denote the residual after the completion of the r-th iteration:

R r (X) = M(X) -M LRA r (X). ( 15
)
The sequence is initiated by setting M LRA 0 (X) = 0 leading to R 0 (X) = M(X). In the r-th correction step, the new rank-one tensor w r is determined by minimizing the empirical error with respect to the current residual:

w r = arg min ω∈W R r-1 -ω 2 E , (16)
where W represents the space of rank-one tensors. Eq. ( 16) is solved by means of an ALS scheme that involves successive minimizations along the dimensions {1, . . . , M }. In the minimization along dimension j, the polynomial coefficients in all other dimensions are "frozen" at their current values; the coefficients z

(j) r = {z (j) 1,r . . . z (j)
pj ,r } are therefore obtained as:

z (j) r = arg min ζ∈R p j +1 R r-1 -C (j) • pj k=0 ζ k P (j) k 2 E , (17)
where C (j) represents the "frozen" component:

C (j) (X 1 , . . . , X j-1 , X j+1 , . . . , X M ) = i =j v (i) r (X i ) = i =j pi k=0 z (i) k,r P (i) k (X i ) . ( 18
)
Because Eq. ( 17) involves only (p j + 1) unknowns (p j < 20 in typical applications), it can be easily solved using the ordinary least squares (OLS) method.

The correction step is initiated by assigning arbitrary values to the unknowns and may involve several iterations over the set of dimensions. Note that assigning initial values to the functions

v (i)
r (Eq. ( 11)) is sufficient; for instance, unity values may be used. Konakli and Sudret Konakli and Sudret (2015a) investigated the effect of the number of iterations performed in a correction step on the accuracy of LRA. They proposed a stopping criterion combining the number of iterations I r with the decrease in the relative empirical error ∆ err r in two successive iterations.

The relative empirical error err r is obtained by normalizing the error measure:

Err r = R r-1 -w r 2 E (19)
with the empirical variance of Y = {M(χ (1) ), . . . , M(χ (N ) )}, the latter denoting the set of model responses at the ED. Accordingly, the algorithm exits the r-th correction step if either I r reaches a maximum allowable value I max or ∆ err r becomes smaller than a prescribed threshold ∆ err min . Based on numerical investigations in different case studies, Konakli and Sudret Konakli and Sudret (2015a) proposed to use I max = 50 and ∆ err min = 10 -6 .

Updating step: After the completion of a correction step, the algorithm moves to an updating step, in which the set of coefficients b = {b 1 . . . b r } is obtained by minimizing the empirical error with respect to the response of the original model:

b = arg min β∈R r M - r l=1 β l w l 2 E . ( 20
)
Note that in each updating step, the size of vector b is increased by one. In the r-th updating step, the value of the new element b r is determined for the first time, whereas the values of the existing elements {b 1 , . . . , b r-1 } are updated (recomputed). Because Eq. ( 20) involves only r unknowns (recall that small ranks are of interest in LRA), it can be easily solved using OLS.

Construction of a rank-R representation in the form of Eq. ( 12) requires repeating pairs of a correction and an updating step for r = 1, . . . , R. The algorithm is summarized below.

Algorithm 1: Construction of a rank-R representation of Y = M(X) with polynomial bases, using a set of input samples E = {χ (1) , . . . , χ (N ) } and the corresponding model responses

Y = {M(χ (1) ), . . . , M(χ (N ) )}:
1. Set R 0 (χ (q) ) = M(χ (q) ), q = 1, . . . , N . (c) While ∆ err r > ∆ err min and I r < I max , repeat steps i-iv:

For

i. Set I r ← I r + 1. ii. Determine z (i) r = {z (i) 0,r . . . z (i) pi,r }, i = 1, .
. . , M , using Eq. ( 17).

iii. Use the current values of z

(i) r , i = 1, . . . , M , to update w r .
iv. Compute Err r using Eq. ( 19) and update ∆ err r .

(d) Determine b = {b 1 . . . b r } using Eq. ( 20).

(e) Evaluate M LRA r (χ (q) ), q = 1, . . . , N , using Eq. ( 13).

(f) Evaluate R r (χ (q) ), q = 1, . . . , N , using Eq. ( 15).

Algorithm 1 describes the construction of LRA for a given rank R. However, in a typical application, the optimal rank is not known a priori. Because Algorithm 1 yields a set of LRA of progressively increasing rank r = 1, . . . , R, the optimal among those can be selected using error-based criteria. In the present study, the optimal LRA is identified by means of 3-fold CV, as proposed by Chevreuil et al. Chevreuil et al. (2013) (see Section 2 for details on kfold CV). Thus, we set r = 1, . . . , r max in Step 2 of Algorithm 1, where r max is a maximum allowable candidate rank, and at the end, select the optimal rank R ∈ {1, . . . , r max } as the one yielding the minimum 3-fold CV error estimate. Konakli and Sudret Konakli and Sudret (2015a) investigated the accuracy of rank selection based on 3-fold CV in different case studies and found that it leads to optimal or nearly optimal LRA in terms of the relative generalization errors, with the latter estimated using large validation sets. We note that the 3-fold CV error estimate may also be used to select the optimal polynomial degrees (see Konakli and Sudret (2015a) for an investigation of the accuracy of this approach).

Comparison to polynomial chaos expansions

A popular method for developing meta-models with polynomial bases is the use of polynomial chaos expansions (PCE). In this section, we provide a brief description of the PCE technique, noting its similarities with LRA.

We consider again the map in Eq. ( 1) assuming that the components of X are independent with a joint PDF f X . Analogously to LRA, the case of dependent input variables can be herein treated with an appropriate isoprobabilistic transformation (see Section 3.1). A PCE approximation of Y = M(X) has the form [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]; [START_REF] Soize | Physical systems with random uncertainties: chaos representations with arbitrary probability measure[END_REF]:

Y PCE = M PCE (X) = α∈A y α Ψ α (X), ( 21
)
where A is a set of multi-indices α = (α 1 , . . . , α M ), {Ψ α , α ∈ A} is a set of multivariate polynomials that are orthonormal with respect to f X and {y α , α ∈ A} is the set of polynomial coefficients. The orthonormal polynomial bases in Eq. ( 21) can be obtained by tensorization of univariate polynomials that are orthonormal with respect to the marginals f Xi :

Ψ α (X) = M i=1 P (i) αi (X i ), (22)
where

P (i)
αi is a univariate polynomial of degree α i in the i-th input variable belonging to an appropriate family. Obviously, the families of the univariate polynomials used to formulate the multivariate PCE basis are the same as the families of polynomials that form the bases of the univariate functions in LRA (see Section 3.1). However, as seen in Eq. ( 12), LRA retain the tensor-product form of Eq. ( 22), whereas the expanded form is considered in PCE. Thus, LRA with polynomial bases can be seen as equivalent compressed representations of PCE.

Different truncation schemes may be employed to determine the set of multi-indices {α ∈ A} in Eq. ( 21). When the maximum degree of P (i) αi is set to p i , i.e. A = {α ∈ N M : α i ≤ p i , i = 1, . . . , M }, the expansion in Eq. ( 21) relies on exactly the same polynomial functions with those used in Eq. ( 12). For this case, let us compare the number of unknowns in LRA and PCE considering a common maximum polynomial degree in all dimensions, i.e. p i = p for i = 1, . . . , M . One has (p + 1) M unknowns in the PCE representation versus (p + 1) • M • R in LRA when redundant parameters are disregarded. Note that the number of unknowns grows exponentially with M in PCE, but only linearly in LRA. For a typical engineering problem with dimensionality M = 10, considering polynomials of low degree p = 3 and an example low rank R = 10, the aforementioned formulas yield 1, 048, 576 PCE coefficients versus a mere 400 unknowns in LRA.

A more efficient truncation scheme is the hyperbolic scheme proposed by Blatman and Sudret [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. This is defined by the condition that the q-norm of any multi-index does not exceed a value p t , i.e. A = {α ∈ N M : α q ≤ p t } with:

α q = M i=1 α i q 1/q , 0 < q ≤ 1. ( 23
)
When q = 1, multivariate polynomials of maximum total degree p t are retained in the expansion.

The corresponding number of terms in the truncated series is:

cardA = M + p t p t = (M + p t )! M !p t ! , (24)
which grows polynomially with M . Smaller values of q impose limitations to the number of terms that include interactions between two or more input variables. Optimal values of p t and q in the hyperbolic truncation scheme can be determined by means of error-based criteria (e.g. the leave-one-out error described later).

Once the basis has been specified, the set of coefficients y = {y α , α ∈ A} may be computed by minimizing the empirical error of the approximation:

y = arg min υ∈R card A M - α∈A υ α Ψ α 2 E . ( 25
)
Even by employing a hyperbolic truncation scheme, the number of unknowns in Eq. ( 25) can be very large in high-dimensional problems, requiring EDs of non-affordable size. Note that contrary to LRA, where the computation of the polynomial coefficients in each dimension is performed separately, the entire set of PCE coefficients is determined from a single minimization problem.

To improve efficiency in the latter, one may substitute Eq. ( 25) with a respective regularized problem. By penalizing the L 1 norm of y, insignificant terms are disregarded from the set of predictors, leading to sparse PCE. An efficient method to solve L 1 -regularized problems is the least angle regression (LAR) method [START_REF] Efron | Least angle regression[END_REF]. A variation proposed by Blatman and Sudret [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] under the name hybrid LAR consists in using the LAR method to determine the best set of predictors and then, computing the PCE coefficients with OLS.

The

EXAMPLE APPLICATIONS

In this section, we confront LRA to sparse PCE in uncertainty propagation through four models with different characteristics and dimensionality. In the first example, we consider a structural- In all applications, the EDs used to build the LRA and PCE meta-models are obtained using Sobol pseudo-random sequences [START_REF] Niederreiter | Random number generation and quasi-Monte Carlo methods[END_REF]. The LRA meta-models are built by implementing Algorithm 1 in Section 3.2. A common maximum polynomial degree p 1 = . . . = p M = p is considered in all dimensions, with its optimal value selected by means of 3-fold CV.

The involved minimization problems are solved using the OLS method. In building the PCE meta-models, a candidate basis is first determined by employing a hyperbolic truncation scheme and then, a sparse expansion is obtained by evaluating the PCE coefficients with the hybrid LAR method, as described in Section 3.3. The optimal combination of the maximum total polynomial degree p t and the parameter q controlling the truncation, where q ∈ {0.25, 0.50, 0.75, 1.0}, is selected as the one leading to the minimum corrected leave-one-out error (see Section 3.3).

The PCE meta-models are built using the UQLab software Marelli andSudret (2014, 2015); an implementation of the algorithm for developing LRA in the same software in currently underway.

Beam deflection

In the first example, we perform reliability analysis of a simply supported beam subjected to a concentrated load at the midspan. The beam has a rectangular cross-section of width b and height h, length L and material Young's modulus E. The magnitude of the concentrated load is denoted by P . The aforementioned parameters are modeled as independent random variables following the distributions listed in Table 1. The response quantity of interest is the midspan deflection, which is obtained through basic structural mechanics as:

U = P L 3 4E bh 3 . ( 26
)
Because U is a product of lognormal random variables, the response PDF can be herein obtained analytically. In particular, U follows a lognormal distribution with parameters (mean and standard deviation of the corresponding normal variable log U ) given by:

λ U = -log(4) + λ P + 3λ L -λ E -λ b -3λ h (27)
and

ζ U = ζ 2 P + 9ζ 2 L + ζ 2 E + ζ 2 b + 9ζ 2 h 1/2 , ( 28
)
where λ Xi and ζ Xi respectively denote the mean and standard deviation of log X i . and3,respectively. In particular, Table 2 lists the rank R and polynomial degree p of the LRA meta-model, the 3-fold CV error estimate err CV3 and the generalization error err G .

Table 3 lists the parameter q controlling the truncation scheme and the total polynomial degree p t of the PCE meta-model, the corrected leave-one-out error err * LOO and the generalization error err G . The generalization errors are estimated using a validations set of size n val = 10 6 sampled with MCS. The ED-based error estimates err CV3 and err * LOO are fairly close to the corresponding generalization errors except for the LRA meta-model when N = 50. In the latter case, err CV3 underestimates err G by approximately one order of magnitude; however, err G is sufficiently small. Note that for each ED, the generalization error of LRA is 2-3 orders of magnitude smaller than that of sparse PCE, which can be justified by the rank-one structure of the herein considered model. N q p t err * LOO err G 30 1 2 3.13 • 10 -2 1.47 • 10 -2 50 1 2 1.56 • 10 -3 1.81 • 10 -3

In Figure 1, we compare the analytical response PDF f U to the respective kernel density estimates (KDEs) obtained with the LRA and sparse PCE meta-models. The KDEs are based on the meta-model responses at a set of n = 10 7 points in the input space sampled with MCS. In Figure 2, we show a similar comparison but using a logarithmic scale in the vertical axis in order to highlight the behavior at the tails of the PDF. Clearly, for the considered EDs of relatively small size, LRA yield superior estimates of the PDF as compared to PCE. It is remarkable that with the LRA approach, an ED of size as small as N = 30 is sufficient to obtain an excellent approximation of f U in the normal scale and a fairly good approximation of the tails.

In the sequel, we use the LRA and sparse PCE meta-models to estimate the failure probability of the actual model in Eq. ( 26). It is worth mentioning that the current implementation of PCE and LRA allows one to sample 10 7 values in a matter of a few seconds using a standard desktop.

P f = P(u lim -U ≤ 0) = P(U ≥ u lim ), i.

Truss deflection

In the second example, we conduct reliability analysis of the truss structure shown in Figure 5 with respect to the midspan deflection u. The random input comprises the six vertical loads, denoted by P 1 , . . . , P 6 , the cross-sectional area and Young's modulus of the horizontal bars, respectively denoted by A 1 and E 1 , and the cross-sectional area and Young's modulus of the vertical bars, respectively denoted by A 2 and E 2 . The distributions of the input random variables are listed in Table 4. The deflection is computed with an in-house finite-element analysis code developed in the Matlab environment.

We develop LRA and sparse PCE meta-models of U = M(A 1 , A 2 , E 1 , E 2 , P 1 , . . . , P 6 ) using two EDs of size N = 50 and N = 100. For both types of meta-models, we use Hermite polynomials to build the basis functions, after an isoprobabilistic transformation of the input variables to standard normal variables. In the LRA algorithm, we define the stopping criterion in the correction step by setting I max = 50 and ∆ err min = 10 -6 . Parameters and error estimates of the LRA and PCE meta-models are listed in Tables 5 and6, respectively. The generalization errors err G are estimated using a validation set of size n val = 10 6 sampled with MCS. For N = 50, the generalization error of LRA is nearly an order of magnitude smaller than that of sparse PCE. For N = 100, the generalization errors of the two types of meta-models are fairly close.

In Figure 6, we compare the KDEs of the response PDF f U obtained with the LRA and sparse PCE meta-models with that obtained with the actual model, which is considered the reference solution for f U . All aforementioned KDEs are based on the evaluation of the different models at a MCS sample of n = 10 6 points in the input space. In Figure 7, the same KDEs are shown in the logarithmic scale in order to emphasize the behavior at the tails. When the LRA approach is employed, the ED of size N = 50 is sufficient to approximate the response PDF with high accuracy in its entire range including the tails. The PCE solution converges more slowly to the reference solution. For N = 50, the discrepancy between the PCE-based KDE from the reference one is obvious even in the normal scale; for N = 100, the PCE-based KDE remains inaccurate for u > 0.11 m (99.1-th percentile).

Table 5: Truss-deflection problem: Parameters and error estimates of LRA meta-models.

N R p err CV3

err G 50 1 2 6.26 • 10 -3 2.85 • 10 -3 100 1 2 3.33 • 10 -3 2.10 • 10 -3

Next, we assess the comparative accuracy of LRA and sparse PCE in estimating the failure probabilities P f = P(u lim -U ≤ 0) = P(U ≥ u lim ), with the deflection threshold u lim varying in [10,15] cm. The LRA-and PCE-based estimates are compared with respective reference values obtained with (i) SORM and (ii) IS (see Section 2.2 for a brief description of these methods). In the IS approach, we utilize the results of a previous analysis with FORM and sequentially add samples of size N IS = 100 until the coefficient of variation of the estimated failure probability becomes smaller than 0.10. The SORM-and IS-based failure probabilities are computed with the software UQLab [START_REF] Marelli | UQLab: a framework for uncertainty quantification in MATLAB[END_REF]; Marelli et al. (2015). The IS-based estimates, considered the reference solution, vary in the range [4.13 • 10 -2 , 3.90 • 10 -6 ]. The LRA-and PCE-based estimates are obtained using a MCS approach with an input sample of size n = 3 • 10 7 , which is sufficient to estimate the smallest failure probability with CoV < 0.10, i.e. with a coefficient of variation similar to that of the reference IS solution. The results are shown in Figure 8. It is remarkable that with the LRA approach, an ED of size as small as N = 50 proves sufficient to evaluate failure probabilities of the order of 10 -6 . The PCE-based estimates converge to the reference solution with increasing N , but at a slower rate than LRA. Figure 9 shows the corresponding ratios of the LRA-and PCE-based reliability indices, β LRA and β PCE

Err C G = M -M 2 X C val . ( 29
)
The conditional generalization error is computed similarly to the generalization error in Eq. ( 3), but by considering only a subset X C val of the validation set X val , defined by an appropriate condition. The corresponding relative error is obtained after normalization with the empirical variance of Y C val , which denotes the set of model responses at X C val . In reliability analysis, we are interested in conditional errors evaluated at subsets of the validation set defined as:

X C val = {x ∈ X val : u = M(x) ≥ u lim }. ( 30
)
In Table 7, we list the relative conditional generalization errors of the LRA and sparse PCE metamodels obtained with N = 100, considering the same values of u lim as in the above reliability analysis. These errors are significantly smaller for LRA than for PCE (about one order of magnitude smaller for the larger response thresholds), which is consistent with Figure 10 and with the results of the reliability analysis. To further highlight the efficiency of the LRA meta-modeling approach, in Table 8, we list the number of evaluations of the actual (finite-element) model required by SORM and IS to compute the failure probability for each threshold. The given number for IS includes the model evaluations required to obtain the FORM estimate. To limit the number of model evaluations in the sequential analyses performed for the increasing thresholds, we start a new FORM analysis from the previous design point. Note that by using the LRA meta-model, we can obtain P f values similar to those computed with SORM and IS, while relying on a much smaller number of model evaluations. We underline that once a meta-model is built, the failure probability for any threshold can be estimated without any additional model evaluations, whereas a new set of model evaluations for each threshold is required by the FORM, SORM and IS techniques.

Heat conduction with spatially varying diffusion coefficient

The present example, inspired by [START_REF] Nouy | Proper generalized decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF], concerns two-dimensional stationary heat-conduction defined on the square domain D = (-0.5, 0.5) m × (0.5, 0.5) m shown in Figure 11. The temper-ature field T (z), z ∈ D, is described by the partial differential equation:

-∇(κ(z)∇T (z)) = QI A (z), (31)
with boundary conditions T = 0 on the top boundary and ∇T •n = 0 on the left, right and bottom boundaries, where n denotes the vector normal to the boundary. In Eq. ( 31), Q = 2 • 10 3 W/m 3 , A = (0.2, 0.3) m×(0.2, 0.3) m is a square domain within D (see Figure 11) and I A is the indicator function equal to 1 if z ∈ A and 0 otherwise. The diffusion coefficient κ(z) is a lognormal random field defined as:

κ(z) = exp[a κ + b κ g(z)], (32)
where g(z) denotes a standard Gaussian random field with autocorrelation function:

ρ(z, z ) = exp (-z -z 2 / 2 ). ( 33
)
In Eq. ( 32), the parameters a κ and b κ are such that the mean and standard deviation of κ are

µ κ = 1 W/C • m and σ κ = 0.3 W/C
• m, respectively, while in Eq. ( 33), = 0.2 m.

To solve Eq. ( 31), the Gaussian random field g(z) in Eq. ( 32) is first discretized using the expansion optimal linear estimation (EOLE) method [START_REF] Li | Optimal discretization of random fields[END_REF]. Let {ζ 1 , . . . , ζ n } denote the points of an appropriately defined grid in D. By retaining the first M terms in the EOLE series, g(z) is approximated by:

g(z) = M i=1 ξ i √ l i φ T i C zζ (z), ( 34
)
where {ξ 1 , . . . , ξ M } are independent standard normal variables; C zζ is a vector with elements

C (k) zζ = ρ(z, ζ k ) for k = 1, .
. . , n; and (l i , φ i ) are the eigenvalues and eigenvectors of the correlation matrix C ζζ with elements [START_REF] Sudret | Stochastic finite element methods and reliability: a state-of-the-art report[END_REF], it is recommended that for a square-exponential autocorrelation function, the size of the element in the EOLE grid must be 1/2 -1/3 of . Accordingly, in the present numerical application, we use a square grid with element size 0.01 m, thus comprising n = 121 points.

C (k,l) ζζ = ρ(ζ k , ζ l ) for k, l = 1, . . . , n. In
The number of terms in the EOLE series is determined according to the rule:

M i=1 l i / n i=1 l i ≥ 0.99, ( 35
)
herein leading to M = 53. The shapes of the first 20 basis functions {φ T i C zζ (z), i = 1, . . . , 20} are shown in Figure 12.

The response quantity of interest is the average temperature in the square domain B = (-0.3, -0.2) m × (-0.3, -0.2) m (see Figure 11), denoted by T :

T = 1 |B| z∈B T (z)dz. ( 36
)
For a given realization of {ξ 1 , . . . , ξ M }, the "exact" model response is obtained with an in-house finite-element analysis code developed in the Matlab environment. The employed finite-element discretization in 16,000 triangular T3 elements is depicted in Figure 11; this discretization is obtained using software Gmsh Geuzaine and Remacle (2009). Figure 13 shows the temperature field T (z) for two example realizations of the conductivity random field. We develop LRA and sparse PCE meta-models of T = M(ξ 1 , , . . . , , ξ 53 ) using EDs of size varying from 50 to 2, 000. Because the random input herein comprises standard normal variables, we build the basis functions by relying on the associated family of Hermite polynomials. In the LRA algorithm, we define the stopping criterion in the correction step by setting I max = 50 and ∆ err min = 10 -6 . Parameters and error estimates of the resulting LRA and PCE metamodels are listed in Tables 9 and10, respectively. The generalization errors err G are estimated using a validation set of size n val = 10 4 sampled with MCS. The two types of meta-models are characterized by generalization errors of the same order of magnitude when a certain ED is considered, except for N = 2, 000. Note that the LRA meta-models exhibit smaller err G than sparse PCE only when N ≤ 200. For both types of meta-models, the ED-based error estimates are well approximated by the respective generalization errors, particularly for the larger EDs.

A B D ∇T • n = 0 ∇T • n = 0 ∇T • n = 0 T = 0 ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡
We first compare the KDEs of the response PDF f T obtained with the LRA and sparse PCE meta-models with that obtained with the actual model, which is considered the reference 01230123the behavior at the tails of the PDF. With only N = 200, the LRA approach yields a good approximation of the reference PDF in both the normal and the logarithmic scales, which is clearly superior to the PCE approximation. For N = 500, the PCE estimate becomes fairly accurate, except for the upper tail of the PDF. We underline that by using the meta-models, we can easily sample larger sets of responses and thus obtain KDEs with smooth tails; for instance, the current implementation of PCE and LRA allows sampling 10 6 values in only a few seconds with a standard desktop. On the other hand, obtaining such large sets of responses by using the actual finite-element model requires prohibitively high computational times (a single evaluation takes approximately 16 sec). We herein consider the same set comprising n = 10 4 points for the evaluation of all KDEs for the sake of comparison. Next, we assess the accuracy of the LRA and sparse PCE meta-models in estimating tail probabilities of the form P f = P( t lim -T ≤ 0) = P( T ≥ t lim ) for the thresholds t lim ∈ {6.0, 6.5} • C.

The reference failure probabilities are obtained with a MCS approach using n = 10 4 evaluations of the actual model, which leads to P f = 4.53 • 10 -2 (β = 1.69) for t lim = 6.0 • C and P f = 1.43 • 10 -2 (β = 2.19) for t lim = 6.5 • C; the CoVs of these estimates are < 0.10. The same input sample is used to estimate the failure probabilities with the two types of meta-models. The above analysis demonstrates that although sparse PCE are characterized by smaller generalization errors for N ≥ 500, the LRA-based estimates of the tail probabilities remain superior. This is because, contrary to PCE, the LRA responses tend to be unbiased at the tails, even in cases when they exhibit a larger dispersion around the actual responses than the PCE ones. This is illustrated in Figure 17 for the case with N = 1, 000. This figure depicts the responses of the LRA and sparse PCE meta-models, denoted by t LRA and t PCE respectively, versus the actual model responses, denoted by t, at the points x of the validation set satisfying the condition t = M(x) ≥ 5.0 • C (27.7-th upper percentile). The LRA responses are characterized by an overall larger dispersion around the actual model responses, which leads to a larger generalization error. However, the LRA responses tend to be unbiased, whereas PCE systematically underestimate the actual model responses at the upper tail. In this case, the conditional generalization errors evaluated at X C val = {x ∈ X val : t = M(x) ≥ t lim } for t lim ∈ {6.0, 6.5} • C (the same thresholds considered in the reliability analysis) remain larger for LRA. By accounting only for the absolute differences between the estimated responses from the actual ones and disregarding their signs, the generalization errors do not reflect herein the superior performance of LRA over PCE in the estimation of tail probabilities. As mentioned earlier, by using the meta-models, we can easily sample larger sets of responses and thus, estimate lower failure probabilities than the above. In Figure 18, we plot the LRAand PCE-based estimates of P f versus the ED size for t lim ∈ {8.0, 8.5} • C. A MCS sample of size n = 2 • 10 6 is used in order to achieve a CoV< 0.10. For both temperature thresholds, the LRA estimates practically reach convergence at N = 500; the corresponding failure probabilities are P f = 3.22 • 10 -4 for t lim = 8.0 • C and P f = 7.35 • 10 -5 for t lim = 8.5 • C. The PCE estimate appears to converge to a similar value with increasing ED size, but at a much slower rate. It is remarkable that an ED of size as small as N = 50 is sufficient to obtain a reasonable preliminary estimate of the failure probability with LRA even for the higher temperature threshold, whereas at least N = 500 points are required to obtain such an estimate with sparse PCE. The values of the LRA-and PCE-based failure probabilities and respective reliability indices for the cases with N = 500 and N = 2, 000 are listed in the Appendix. We emphasize that for the considered thresholds, reference solutions that rely on the evaluation of the actual finite-element model cannot be obtained at an affordable computational time. Note that because of the high dimensionality of the problem, the typically low-cost FORM and SORM approaches are herein inefficient.

Frame displacement

In the last example, we consider the frame structure shown in Figure 19, also studied in [START_REF] Liu | Optimization algorithms for structural reliability[END_REF]; [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. We conduct reliability analysis with respect to the horizontal displacement u at the top right corner of the top floor under the depicted horizontal loads acting at the floor levels. The random input comprises the load values P 1 , P 2 and P 3 , the Young's moduli of the column and beam elements, respectively denoted by E C and E B , the moments of inertia of the column and beam elements, respectively denoted by I Bi and I Ci , i = 1, . . . , 4, and the cross-sectional areas of the column and beam elements, respectively denoted by I Bi and I Ci , i = 1, . . . , 4. The distributions of the aforementioned variables are listed in Table 11. Contrary to the previous examples where the input variables were independent, the input of the herein considered model has a dependence structure described by means of a Gaussian copula (for further information on the modeling of probabilistic dependence with copulas, the interested reader is referred to [START_REF] Nelsen | An introduction to copulas[END_REF]). The non-zero elements of the associated linear correlation matrix are defined as follows: the correlation coefficient between the two Young's moduli is ρ EC ,EB = 0.90; the correlation coefficient between the cross-sectional area A i and the moment of inertia I i of a certain element i is ρ Ai,Ii = 0.95; the correlation coefficient between the geometric properties of two distinct elements i and j are ρ Ai,Ij = ρ Ii,Ij = ρ Ai,Aj = 0.13. In the original example, the above values represent the corresponding linear correlation coefficients in the standard normal space; however, Blatman and Sudret [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] note that the differences between the two are insignificant. For a given realization of the input random vector, the frame displacement is computed with an in-house finite-element analysis code developed in the Matlab environment. We develop LRA and sparse PCE meta-models of X = {P 1 , P 2 , P 3 , E C , E B , I C1 , . . . , I B4 , A C1 , . . . , A B4 } using two EDs of size N = 500 and N = 1, 000. For both types of meta-models, we use Hermite polynomials to build the basis functions, after an isoprobabilistic transformation of the input variables to independent standard normal variables. In the LRA algorithm, we define the stopping criterion in the correction step by setting I max = 50 and ∆ err min = 10 -6 . Parameters and error estimates of the LRA and PCE meta-models are listed in Tables 12 and13, respectively.

P 1 P 1 P 1 P 2 P 3 B 1 B 1 B 1 B 1 B 2 B 2 B 2 B 2 B 2 B 2 B 3 B 3 B 3 B 3 B 4 C 3 C 3 C 3 C 3 C 4 C 4 C 1 C 1 C 1 C 1 C 2 C 2 C 2 C 2 C 2 C 2 C 2
The generalization errors err G are estimated using a validation set of size n val = 10 6 sampled with MCS. For N = 500, the two types of meta-models exhibit similar generalization errors, while for N = 1, 000, the PCE error is slightly smaller. Note that the generalization errors are approximated fairly well by the corresponding ED-based error estimates.

In Figure 20, we assess the accuracy of the LRA and sparse PCE meta-models in estimating the failure probability P f = P(u lim -U ≤ 0) = P(U ≥ u lim ), with the displacement threshold u lim varying in the range [3,8] -6 ]; for the largest displacement thresholds, these tend to be underestimated when the actual model is replaced by the metamodels. The LRA estimates are characterized by higher accuracy than the PCE ones and remain within the order of magnitude of the reference value, even for the smallest ED and the largest displacement threshold considered. Note again the superior performance of LRA in the prediction of extreme responses, despite the smaller generalization error of PCE. Figure 21 shows the corresponding ratios of the reliability indices β LRA and β PCE , based on the LRA and sparse PCE meta-models, to the reliability index β, obtained in terms of the reference P f . The latter varies in the range [1.36, 4.33]. When the LRA approach is employed, the relative error in β remains smaller than 5% for N = 500 and smaller than 4% for N = 1, 000. The PCE errors are slightly higher; for the largest displacement threshold, they exceed 10% and 5% when N = 500 To highlight the computational gain achieved by using a meta-modeling approach, in Table 14, we list the number of evaluations of the actual model (for each threshold) required to compute the reference failure probability. The given numbers include the model evaluations used to obtain the FORM estimate; as in the truss-deflection problem, the FORM analysis for each threshold, except for the smallest, is initiated from the design point corresponding to the previous threshold. The total computational cost of the reference solution comprises the sum of the listed model evaluations for all considered thresholds. We underline that this cost can increase fast with decreasing target coefficient of variation (note the orders-of-magnitude larger number of model evaluations in the present example, requiring CoV< 0.01, as compared to that in the truss-deflection problem, requiring CoV< 0.10). Conversely, when the analysis relies on the meta-models, IS is essentially costless, independently of the target coefficient of variation, and the computational effort is confined to the evaluations of the actual model at the ED. As seen above, LRA herein provides an estimate of a reliability index higher than 4 with an error smaller than 5% based on an ED as small as N = 500.

E B [KN/m 2 ] Truncated Gaussian over [0, ∞) 2.1738 • 10 7 1.9152 • 10 6 I C 1 [m 4 ] Truncated Gaussian over [0, ∞) 8.1344 • 10 -3 1.0834 • 10 -3 I C 2 [m 4 ] Truncated Gaussian over [0, ∞) 1.1509 • 10 -2 1.2980 • 10 -3 I C 3 [m 4 ] Truncated Gaussian over [0, ∞) 2.1375 • 10 -2 2.5961e -03 I C 4 [m 4 ] Truncated Gaussian over [0, ∞) 2.5961 • 10 -2 3.0288 • 10 -3 I B 1 [m 4 ] Truncated Gaussian over [0, ∞) 1.0811 • 10 -2 2.5961 • 10 -3 I B 2 [m 4 ] Truncated Gaussian over [0, ∞) 1.4105 • 10 -2 3.4615 • 10 -3 I B 3 [m 4 ] Truncated Gaussian over [0, ∞) 2.3279 • 10 -2 5.6249 • 10 -3 I B 4 [m 4 ] Truncated Gaussian over [0, ∞) 2.5961 • 10 -2 6.4902 • 10 -3 A C 1 [m 2 ] Truncated Gaussian over [0, ∞) 3.1256 • 10 -1 5.5815 • 10 -2 A C 2 [m 2 ] Truncated Gaussian over [0, ∞) 3.7210 • 10 -1 7.4420 • 10 -2 A C 3 [m 2 ] Truncated Gaussian over [0, ∞) 5.0606 • 10 -1 9.3025 • 10 -2 A C 4 [m 2 ] Truncated Gaussian over [0, ∞) 5.5815 • 10 -1 1.1163 • 10 -1 A B 1 [m 2 ] Truncated Gaussian over [0, ∞) 2.5302 • 10 -1 9.3025 • 10 -2 A B 2 [m 2 ] Truncated Gaussian over [0, ∞) 2.9117 • 10 -1 1.0232 • 10 -1 A B 3 [m 2 ] Truncated Gaussian over [0, ∞) 3.7303 • 10 -1 1.2093 • 10 -1 A B 4 [m 2 ] Truncated Gaussian over [0, ∞) 4.1860 • 10 -1 1.9537 • 10 -1

CONCLUSIONS

Reliability analysis faces challenges in cases when the systems under consideration are represented by complex high-dimensional computational models. In this paper, we demonstrate that meta-models belonging to the class of canonical low-rank approximations (LRA) can provide an accurate representation of the probability density function (PDF) of the model response at the tails, thus leading to efficient estimation of the small exceedence probabilities required in reliability analysis. By replacing a complex model by a meta-model that possesses similar statistical properties, evaluation of a response quantity of interest becomes essentially costless from a computational viewpoint. The LRA approach can be particularly efficient in high-dimensional problems because: (i) the number of unknowns grows only linearly with the input dimension and (ii) their construction relies on a series of least-square minimization problems of small size that is independent of the input dimension.

In this paper, canonical LRA developed with polynomial bases are of interest because of the simplicity and versatility characterizing the use of polynomial functions. The formulation and construction of such LRA meta-models in a non-intrusive manner is detailed. Furthermore, the links between canonical LRA and the popular meta-modeling technique of polynomial chaos expansions (PCE) are explained. Canonical LRA are confronted to sparse PCE in reliability applications involving a rank-one model (for which an analytical solution is available) and three finite-element models pertinent to structural mechanics and heat conduction. In all aforementioned applications, canonical LRA are found to outperform sparse PCE for cases when the size of the experimental design is relatively small with respect to the input dimension. By providing unbiased representations of the model responses at the tails, canonical LRA provide superior estimates of small exceedence probabilities compared to sparse PCE, even in cases when the latter exhibit smaller generalization errors. In the examined applications, failure probabilities of the order of 10 -5 are predicted with sufficient accuracy by canonical LRA based on 5M -25M evaluations of the original model, where M denotes the dimension of the random input. The LRA approach is also found to outperform methods particularly targeted to reliability analysis in terms of the required number of model evaluations. We underline that contrary to such methods, LRA provide a full probabilistic description of the model response, which can be used to estimate any statistical measure of interest beyond the probabilities of exceedence.

Having introduced canonical LRA in reliability analysis and demonstrated its strong potential for dealing with high dimensionality, we underline the need for further studies that will establish the efficacy of the approach in diverse reliability applications. The construction and use of tensor approximation of different formats, beyond the herein considered canonical formulation, is an active research topic in the field of uncertainty quantification. We hope that the present study will motivate further investigations into the capacities of such meta-models to accurately predict extreme responses of high-dimensional models, thus opening new paths to the risk assessment of complex systems.

  r = 1, . . . , R, repeat steps (a)-(f): (a) Assign initial values to v (i) r , i = 1, . . . , M (e.g. unity values). (b) Set I r = 0 and ∆ err r = > ∆ err min .

  mechanics model of dimension M = 5 having an analytical rank-one structure. The following three examples involve finite-element models; in particular, we consider a truss model with independent input of dimension M = 10, a heat-conduction model with thermal conductivity described by a random field, which is approximated by a series expansion of dimension M = 53, and a frame model with correlated input of dimension M = 21. For the aforementioned models, we investigate the comparative accuracy of LRA and PCE in the estimation of small failure probabilities P f and of the corresponding reliability indices β = -Φ -1 ( P f ), where Φ denotes the standard normal cumulative distribution function (CDF).

Figure 1 :Figure 2 :

 12 Figure 1: Beam-deflection problem: Probability density function of the response (normal scale).

Figure 3 Figure 3 :

 33 Figure 3 depicts the estimates of the failure probability for the two EDs of size N = 30 and N = 50 together with the respective analytical solutions. For N = 30, the LRA-based estimates are close to the analytical solutions, especially for the smaller deflection thresholds, whereas the PCE-based estimates are overall highly inaccurate. For N = 50, the LRA-based estimates are excellent in the entire range examined, whereas the PCE-based estimates remain poor for the larger deflection thresholds. Figure 4 shows the corresponding ratios of the reliability index estimates based on the LRA and sparse PCE meta-models, respectively denoted by β LRA and β PCE , to the reliability index based on the analytical solution, denoted by β. For the considered deflection thresholds, β varies in the range [1.51, 4.25]. Note that with the LRA approach, we estimate the largest β with a relative error smaller than 2% using only N = 30 evaluations of the actual model. The values of the failure probabilities and corresponding reliability indices

Figure 4 :

 4 Figure 4: Beam-deflection problem: Ratios of meta-model-based to reference reliability indices.

Figure 5 :

 5 Figure 5: Truss structure.

Figure 6 :Figure 7 :

 67 Figure 6: Truss-deflection problem: Probability density function of the response (normal scale).

Figure 8 :

 8 Figure 8: Truss-deflection problem: Failure probabilities.

Figure 9 :

 9 Figure 9: Truss-deflection problem: Ratios of meta-model-based to reference reliability indices.

Figure 10 :

 10 Figure 10: Truss-deflection problem: Meta-model versus actual model responses at a subset of the validation set corresponding to the upper tail of the response distribution.

  Figure 11: Heat-conduction problem: Domain and boundary conditions (left); finite-element mesh (right).

Figure 12 :

 12 Figure 12: Heat-conduction problem: Shapes of the first 20 basis functions in the EOLE discretization (from left-top to bottom-right row-wise).

Figure 13 :

 13 Figure 13: Heat-conduction problem: Example realizations of the temperature field. Table 9: Heat-conduction problem: Parameters and error estimates of LRA meta-models. N R p err CV3 err G 50 3 1 6.51 • 10 -2 1.08 • 10 -1 100 1 1 3.56 • 10 -2 2.46 • 10 -2 200 1 1 1.79 • 10 -2 1.38 • 10 -2 500 1 2 1.13 • 10 -2 9.68 • 10 -3 1,000 1 2 8.40 • 10 -3 8.19 • 10 -3 2,000 1 2 7.81 • 10 -3 7.72 • 10 -3

Figure 14 :

 14 Figure 14: Heat-conduction problem: Probability density function of the response (normal scale).

Figure 15 :Figure 16 :

 1516 Figure 15: Heat-conduction problem: Probability density function of the response (log-scale).

Figure 17 :

 17 Figure 17: Heat-conduction problem: Meta-model versus actual model responses at a subset of the validation set corresponding to the upper tail of the response distribution.

Figure 18 :

 18 Figure 18: Heat-conduction problem: Failure probabilities.

Figure 19 :

 19 Figure 19: Frame structure.

  cm. All failure probabilities are herein computed with the IS technique. The reference values P f are obtained by employing IS in conjunction with the actual finite-element model; the LRA and PCE values are computed with exactly the same algorithm but using the respective meta-models in lieu of the original model. As in the truss-deflection problem, the auxiliary PDF in IS is defined in terms of the design point indicated by a previous FORM analysis. Samples of size N IS = 1, 000 are then sequentially added until the coefficient

Figure 20 :

 20 Figure 20: Frame-displacement problem: Failure probabilities.

Figure 21 :

 21 Figure 21: Frame-displacement problem: Ratios of meta-model-based to reference reliability indices.

t

  lim [ • C] Reference LRA PCE N = 500 N = 2, 000 N = 500 N = 2, 000 6 4.53 • 10 -2 4.57 • 10 -2 4.59 • 10 -2 4.13 • 10 -2 4.57 • 10 -2 6.5 1.43 • 10 -2 1.42 • 10 -2 1.44 • 10 -2 1.15 • 10 -2 1.36 • 10 -2 8 -3.22 • 10 -4 3.23 • 10 -4 1.23 • 10 -4 2.27 • 10 -4 8.5 -7.35 • 10 -5 7.40 • 10 -5 1.65 • 10 -5 4.45 • 10 -5

  PCE accuracy can be assessed by means of the leave-one-out error Err LOO , corresponding to the CV error for the extreme case k = N[START_REF] Allen | The prediction sum of squares as a criterion for selecting predictor variables[END_REF]. Using algebraic manipulations, this error can be computed based on a single PCE that is built with the full ED (see[START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] for details). The corresponding relative error, denoted by err LOO , is obtained after normalizing Err LOO with the empirical variance of Y = {M(χ (1) ), . . . , M(χ (N ) )}. Be-

cause err LOO can be too optimistic,

[START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]

proposed the use of the corrected leave-one-out error, which includes a multiplication factor derived by

[START_REF] Chapelle | Model selection for small sample regression[END_REF]

.

Table 1 :

 1 Beam-deflection problem: Distributions of input variables.

	Variable Distribution mean CoV
	b [m]	Lognormal	0.15	0.05
	h [m]	Lognormal	0.3	0.05
	L [m]	Lognormal	5	0.01
	E [MPa] Lognormal 30,000 0.15
	P [KN]	Lognormal	10	0.20

Table 2 :

 2 Beam-deflection problem: Parameters and error estimates of LRA meta-models.

	N R p	err CV3	err G
	30 1 2 1.21 • 10 -4 2.32 • 10 -4 50 1 3 3.14 • 10 -7 2.63 • 10 -6

Table 3 :

 3 Beam-deflection problem: Parameters and error estimates of PCE meta-models.

Table 4 :

 4 Truss-deflection problem: Distributions of input random variables.

	Variable	Distribution	mean	CoV
	A 1 [m]	Lognormal	0.002	0.10
	A 2 [m]	Lognormal	0.001	0.10
	E 1 , E 2 [MPa]	Lognormal 210,000 0.10
	P 1 , . . . , P 6 [KN]	Gumbel	50	0.15

Table 6 :

 6 Truss-deflection problem: Parameters and error estimates of PCE meta-models.

	N 50 0.25 2 4.25 • 10 -2 1.24 • 10 -2 q p t err * err G LOO 100 1 2 3.50 • 10 -3 2.56 • 10 -3

Table 7 :

 7 Truss-deflection problem: Relative conditional generalization errors of the LRA and PCE meta-models obtained with N = 100 for different response thresholds.

	u lim [cm]	LRA	PCE
	10 11 12 13 14 15	2.89 • 10 -2 8.68 • 10 -2 5.70 • 10 -2 2.58 • 10 -1 1.04 • 10 -1 6.45 • 10 -1 2.00 • 10 -1 1.45 • 10 0 3.94 • 10 -1 3.55 • 10 0 2.67 • 10 0 2.62 • 10 1
		21

Table 8 :

 8 Truss-deflection problem: Number of model evaluations required in the computation of failure probabilities with SORM and IS.

	u lim [cm] SORM IS
	10	387	475
	11	285	473
	12	297	585
	13	309	597
	14	321	709
	15	333	921

Table 10 :

 10 Heat-conduction problem: Parameters and error estimates of PCE meta-models.

	N 50	q 0.25 2 1.71 • 10 -1 2.53 • 10 -1 p t err * err G LOO

Table 11 :

 11 Frame-displacement problem: Distributions of input random variables.

	Variable	Distribution	Mean	Standard deviation
	P 1 [KN]	Lognormal	133.45	40.04
	P 2 [KN]	Lognormal	88.97	35.59
	Lognormal E C [KN/m 2 ] Truncated Gaussian over [0, ∞) 2.3796 • 10 7 71.17 P 3 [KN]	28.47 1.9152 • 10 6

Table 12 :

 12 Frame-displacement problem: Parameters and error estimates of LRA meta-models.

	N	R p	err CV3	err G
	500 1,000 1 3 2.73 • 10 -3 3.19 • 10 -3 1 3 2.93 • 10 -3 3.35 • 10 -3

Table 13 :

 13 Frame-displacement problem: Parameters and error estimates of PCE meta-models.

	N 500 0.50 7 1.67 • 10 -3 3.06 • 10 -3 q p t err * err G LOO 1,000 1 3 9.37 • 10 -4 1.42 • 10 -3
	of variation of the estimated probability becomes smaller than 0.01. Values of the so-obtained
	reference failure probability vary in the range [8.72•10 -2 , 7.38•10

Table 14 :

 14 Frame-displacement problem: Number of model evaluations required in the computation of the reference failure probabilities.

	u lim [cm]	IS
	3	18,340
	4	29,116
	5	41,162
	6	58,185
	7	79,231
	8	132,300

Table A2 :

 A2 Beam-deflection problem: Reliability indices (see Figure4).

	u lim [mm]	Analytical	LRA N = 30 N = 50 N = 30 N = 50 PCE
	4	1.51	1.51	1.51	1.58	1.52
	5	2.26	2.27	2.26	2.52	2.36
	6	2.88	2.91	2.88	3.35	3.09
	7	3.40	3.45	3.41	4.10	3.74
	8	3.85	3.91	3.86	4.89	4.35
	9	4.25	4.32	4.26	-	4.99

Table A3 :

 A3 Truss-deflection problem: Failure probabilities (see Figure8). 10 -2 4.13 • 10 -2 4.35 • 10 -2 4.23 • 10 -2 3.17 • 10 -2 3.83 • 10 -2 11 9.28 • 10 -3 9.86 • 10 -3 8.88 • 10 -3 8.47 • 10 -3 3.66 • 10 -3 6.45 • 10 -3 12 1.63 • 10 -3 1.34 • 10 -3 1.52 • 10 -3 1.41 • 10 -3 2.58 • 10 -4 8.19 • 10 -4 13 2.57 • 10 -4 2.14 • 10 -4 2.28 • 10 -4 2.09 • 10 -4 1.20 • 10 -5 8.03 • 10 -5 14 3.77 • 10 -5 3.46 • 10 -5 3.23 • 10 -5 2.73 • 10 -5 4.00 • 10 -7 6.03 • 10 -6 15 5.31 • 10 -6 3.90 • 10 -6 4.73 • 10 -6 3.10 • 10 -6 -4.33 • 10 -7

	u lim [cm]	Reference SORM	IS	N = 50	LRA N = 100	N = 50	PCE N = 100
	10	4.47 •

Table A4 :

 A4 Truss-deflection problem: Reliability indices (see Figure9).

	u lim [cm]	Reference SORM IS N = 50 N = 100 N = 50 N = 100 LRA PCE
	10	1.70	1.74	1.71	1.72	1.86	1.77
	11	2.35	2.33	2.37	2.39	2.68	2.49
	12	2.94	3.00	2.96	2.99	3.47	3.15
	13	3.47	3.52	3.51	3.53	4.22	3.77
	14	3.96	3.98	4.00	4.04	4.94	4.38
	15	4.40	4.47	4.43	4.52	-	4.92

Table A5 :

 A5 Heat-conduction problem: Failure probabilities (see Figures16 and 18).

Table A6 :

 A6 Heat-conduction problem: Reliability indices (see Figures16 and 18).

	t lim [ • C]	Reference	LRA N = 500 N = 2, 000 N = 500 N = 2, 000 PCE
	6	1.69	1.69	1.69	1.74	1.69
	6.5	2.19	2.19	2.19	2.27	2.21
	8	-	3.41	3.41	3.67	3.51
	8.5	-	3.80	3.79	4.15	3.92
			40

Table A7 :

 A7 Frame displacement: Failure probabilities (see Figure20). = 1, 000 N = 500 N = 1, 000 3 8.72 • 10 -2 8.51 • 10 -2 8.69• 10 -2 8.54 • 10 -2 8.75 • 10 -2 4 1.15 • 10 -2 1.10 • 10 -2 1.12 • 10 -2 1.04 • 10 -2 1.08 • 10 -2 5 1.53 • 10 -3 1.32 • 10 -3 1.35 • 10 -3 1.07 • 10 -3 1.22 • 10 -3 6 2.27 • 10 -4 1.61 • 10 -4 1.70 • 10 -4 1.04 • 10 -4 1.31 • 10 -4 7 3.78 • 10 -5 2.10 • 10 -5 2.32 • 10 -5 9.16 • 10 -6 1.45 • 10 -5 8 7.38 • 10 -6 2.94 • 10 -6 3.32 • 10 -6 7.36 • 10 -7 1.58 • 10 -6

	u lim	Reference	LRA	PCE
	[cm]		N = 500 N

Table A8 :

 A8 Frame displacement: Reliability indices (see Figure21).

	u lim	Reference	LRA		PCE
	[cm]		N = 500 N = 1, 000 N = 500 N = 1, 000
	3	1.36	1.37	1.36	1.37	1.36
	4	2.27	2.29	2.28	2.31	2.30
	5	2.96	3.01	3.00	3.07	3.03
	6	3.51	3.60	3.58	3.71	3.65
	7	3.96	4.10	4.07	4.28	4.18
	8	4.33	4.53	4.51	4.82	4.66

We develop LRA and sparse PCE meta-models of U = M(P, L, E, b, h) using two EDs of size N = 30 and N = 50. For both types of meta-models, we use Hermite polynomials to build the basis functions, after an isoprobabilistic transformation of the input variables to standard normal variables. In the LRA algorithm, we define the stopping criterion in the correction step by setting I max = 50 and ∆ err min = 10 -8 (it was shown inKonakli and Sudret (2015a) that in the considered problem, selecting a small value for ∆ err min is critical for the LRA accuracy). Parameters and error estimates of the LRA and PCE meta-models are listed in

Appendix

Table A1: Beam-deflection problem: Failure probabilities (see Figure 3).

2.00 • 10 -3 1.81 • 10 -3 1.96 • 10 -3 4.00 • 10 -4 1.00 • 10 -3 7 3.37 • 10 -4 2.80 • 10 -4 3.25 • 10 -4 2.07 • 10 -5 9.08 • 10 -5 8 5.86 • 10 -5 4.57 • 10 -5 5.57 • 10 -5 5.00 • 10 -7 6.70 • 10 -6 9 1.07 • 10 -5 7.90 • 10 -6 1.01 • 10 -5 -3.00 • 10 -7