C Lataniotis

S Marelli

B Sudret

Extending classical surrogate modelling to ultrahigh dimensional problems through supervised dimensionality reduction: a data-driven approach

Keywords:

Thanks to their versatility, ease of deployment and high-performance, surrogate models have become staple tools in the arsenal of uncertainty quantification (UQ). From local interpolants to global spectral decompositions, surrogates are characterised by their ability to efficiently emulate complex computational models based on a small set of model runs used for training. An inherent limitation of many surrogate models is their susceptibility to the curse of dimensionality, which traditionally limits their applicability to a maximum of O(10 2 ) input dimensions. We present a novel approach at high-dimensional surrogate modelling that is model-, dimensionality reduction-and surrogate model-agnostic (black box), and can enable the solution of high dimensional (i.e. up to O(10 4 )) problems. After introducing the general algorithm, we demonstrate its performance by combining Kriging and polynomial chaos expansions surrogates and kernel principal component analysis. In particular, we compare the generalisation performance that the resulting surrogates achieve to the classical sequential application of dimensionality reduction followed by surrogate modelling on several benchmark applications, comprising an analytical function and two engineering applications of increasing dimensionality and complexity.

Introduction

It is nowadays a common practice to study the behaviour of physical and engineering systems through computer simulation. In a real-world setting, such systems are driven by input parameters, the values of which can be uncertain or even unknown. Uncertainty quantification (UQ) aims at identifying and quantifying the sources of uncertainty in the input parameters to assess the uncertainty they cause in the model predictions. In the context 1 of Monte Carlo simulation, such workflow typically entails the repeated evaluation of the computational model. However, it may become intractable when a single simulation is computationally demanding, as is often the case with modern computer codes. A remedy to this problem is to substitute the model with a surrogate that accurately mimics the model response within the chosen parameter bounds, but is computationally inexpensive. An additional benefit of surrogate models is that they are often non-intrusive, i.e. their construction only depends on a training set of model evaluations, without access to the model itself. This includes the case when the model is not available, but only a pre-existent data set is, as is typical in machine learning applications. The latter setting is the focus of this paper. Popular surrogate modelling techniques (SM) include Gaussian process modelling and regression [START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], polynomial chaos expansions [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF][START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF][START_REF] Xiu | Numerical methods for stochastic computations -A spectral method approach[END_REF], low-rank tensor approximations [START_REF] Chevreuil | A least-squares method for sparse low rank approximation of multivariate functions[END_REF]Konakli and Sudret, 2016b), and support vector regression [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. Parametrising and training a surrogate model, however, can become harder or even intractable as the number of input parameters increases, a well known problem often referred to as curse of dimensionality (see e.g. [START_REF] Verleysen | The curse of dimensionality in data mining and time series prediction[END_REF]).

For the sake of clarity, in the following we will classify high-dimensional inputs in two broad categories, depending on their characteristics: unstructured or structured. Unstructured inputs are characterised by the lack of an intrinsic ordering, and they are commonly identified with the so-called "model parameters", e.g. point loads on mechanical models, or resistance values in electrical circuit models. Structured inputs, on the other hand, are characterised by the existence of a natural ordering and/or a distance function (i.e. they show strong correlation across some physically meaningful set of coordinates), as it is typical for time-series or space-variant quantities represented by maps. Boundary conditions in complex simulations that rely on discretisation grids, e.g. time-dependent excitations at grid nodes, often belong to this second class. In most practical applications, unstructured inputs range in dimension in the order O(10 0-2 ), while structured inputs tend to be in the order O(10 2-6 ).

Several strategies have been explored in the literature to deal with high dimensional problems for surrogate modelling. A common approach in dealing with unstructured inputs is input variable selection, which consists in identifying the "most important" inputs according to some importance measure, see e.g. Saltelli et al. (2008); [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF], and simply ignoring the others (e.g. by setting them to their nominal value).

In the context of kernel-based emulators (e.g. Gaussian process modelling or support vector machines), some attention has been devoted to the use of simple isotropic kernels [START_REF] Djolonga | High-dimensional Gaussian process bandits[END_REF], or to the design of specific kernels for high-dimensional input vectors, sometimes including deep-learning techniques (e.g., [START_REF] Lawrence | Probabilistic non-linear principal component analysis with Gaussian process latent variable models[END_REF]; [START_REF] Durrande | Additive covariance kernels for high-dimensional Gaussian process modeling[END_REF]; [START_REF] Wilson | Deep kernel learning[END_REF]).

In more complex scenarios, the more general concept of dimensionality reduction (DR) is applied, which essentially consists in mapping the input space to a suitable lower dimensional space using an appropriate transformation prior to the surrogate modelling stage. The latter approach is considered in this work due to its applicability to cases for which variable selection seems inadequate or insufficient (e.g. in the presence of structured inputs).

In the current literature, a two-step approach is often followed for dealing with such problems: first, the input dimension is reduced; then, the surrogate model is constructed directly in the reduced (feature-) space. The dimensionality reduction step is based on an unsupervised objective, i.e. an objective that only takes into account the input observations.

Examples of unsupervised objectives include the minimisation of the input reconstruction error [START_REF] Vincent | Extracting and composing robust features with denoising autoencoders[END_REF], maximisation of the sample variance [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF], maximisation of statistical independence [START_REF] Hyvärinen | One-unit learning rules for independent component analysis[END_REF], and preservation of the distances between the observations [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Hinton | Stochastic neighbor embedding[END_REF]. While in principle attractive due to their straightforward implementation, unsupervised approaches for dimensionality reduction may be suboptimal in this context, because the input-output map of the reduced representation may exhibit a complex topology unsuitable for surrogate modelling [START_REF] Wahlström | Learning deep dynamical models from image pixels[END_REF][START_REF] Calandra | Manifold Gaussian processes for regression[END_REF].

To deal with this issue, various supervised techniques have been proposed, in the sense that the objective of the input compression takes into account the model outputs. One such approach that has received attention recently is based on the so-called active subspaces concept [START_REF] Constantine | Active subspace methods in theory and practice: applications to kriging surfaces[END_REF]. Various methods that belong to this category, provide a linear transformation of the high dimensional input space into a reduced space that is characterised by maximal variability w.r.t. the model output. However, active subspace methods often require the availability of the model gradient w.r.t. the input parameters, a limiting factor in data-driven scenarios where such information is not available and needs to be approximated [START_REF] Fornasier | Learning functions of few arbitrary linear parameters in high dimensions[END_REF]. Moreover, the numerical computation of the gradient may be infeasible in problems that involve structured inputs such as time series or 2D maps with O(10 2-6 ) components.

Other data-driven supervised DR techniques have been proposed in the literature, that are dependent on the properties of a specific combination of either DR or SM techniques. [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] employ multi-layer neural networks for both the DR and the SM steps. Specifically, an unsupervised objective based on the reconstruction error is followed by a generalisation performance objective that aims at fine tuning the network weights with respect to a measure of the surrogate modelling error. Similar approaches have been proposed with other combinations of methods. In [START_REF] Damianou | Deep Gaussian processes[END_REF], the same idea is extended by using stacked Gaussian processes instead of multilayer neural networks. In [START_REF] Huang | Scalable Gaussian process regression using deep neural networks[END_REF]; [START_REF] Calandra | Manifold Gaussian processes for regression[END_REF] this approach is extended by combining neural networks with Gaussian processes within a Bayesian framework.

All of these methods demonstrate that supervised methods yield a significant accuracy advantage over the unsupervised ones, as the final goal of the supervised learner (i.e. surrogate model accuracy) matches the final goal of high-dimensional surrogate modelling in the first place. However, this increased accuracy comes at the cost of restricting the applicability of such methods to specific combinations of DR and SM techniques.

In this paper, we propose a novel method of performing dimensionality reduction for surrogate modelling in a data-driven setting, which we name (perhaps with a lack of creative flair) DRSM. The aim of this method is to capitalise on the performance gains of supervised DR, while maintaining maximum flexibility in terms of both DR and SM methodologies.

Recognising that different communities, applications and researchers have in general access to one or two preferred techniques for either DR or SM, the proposed approach is fully non-intrusive, i.e. both the DR and the SM stages are considered as black boxes under very general conditions. The novelty lies in the way the two stages are coupled into a single problem, for which dedicated solvers are proposed. This paper is structured as follows: Section 2 introduces the main ingredients required by DRSM, namely dimensionality reduction and surrogate modelling. For the sake of clarity, some of the techniques that will be specifically used in the applications section are also introduced, i.e. kernel principal component analysis (KPCA) for DR, Gaussian process modelling, a.k.a. Kriging, and polynomial chaos expansions (PCE) for SM. The core framework underlying DRSM is then introduced. Finally, the effectiveness of DRSM is analysed on several benchmark applications including both unstructured and structured inputs, ranging from low-dimensional analytical functions to a complex engineering 2-dimensional heat-transfer problem.

Ingredients for surrogate modelling in high dimension

As the name implies, DRSM consists in the combination of two families of computational tools: dimensionality reduction and surrogate modelling. This section aims at highlighting the main features of each, and how they can be exploited without resorting to intrusive, dedicated algorithms.

Dimensionality reduction

Consider a set of high-dimensional samples X = x (i) ∈ R M , i = 1, . . . , N . In an abstract sense, dimensionality reduction (DR) refers to the parametric mapping g :

X ∈ R M → Z ∈ R m of the form: z = g(x; w) (1)
where z ∈ Z, x ∈ X , and w is the set of parameters associated with the mapping. Di-

mensionality reduction occurs if m M , i.e. if m = O 10 0-1 whereas M = O 10 2-4 .
The nature and number of the parameters w depends on the specific DR method under consideration.

Such transformations are motivated by the assumption that the samples in X lie on some manifold with dimensionality m that is embedded within the M -dimensional space. This specific value of m is in some applications referred to as the "intrinsic dimension" of X [START_REF] Fukunaga | Introduction to statistical pattern recognition[END_REF]. From an information theory perspective, the intrinsic dimension refers to the minimum number of scalars that is required to represent X without any loss w.r.t. an appropriate information measure. In practice it is a-priory unknown. In such cases DR is an ill-posed problem that can only be solved by assuming certain properties of X , such as its intrinsic dimension. Alternatively the later may be approximated and/or inferred from the available data by various approaches (see e.g. [START_REF] Camastra | Data dimensionality estimation methods: a survey[END_REF] for a comparative overview).

An important aspect of all parametric DR methods, regardless of their specificity, is that for each choice of dimension m the remaining parameters w are estimated by minimising a suitable error measure (sometimes referred to as loss function):

w = arg min Dw J(w; X ), (2)
where w denotes the estimated parameters, D w the feasible domain of w, J(•) the error measure and X the available data. The choice of the error measure depends on the specific application DR is used for. When the goal is direct compression of a high dimensional input without information loss (a common situation in telecommunication-related applications), a typical choice of J(•) is the so-called mean-squared reconstruction error, that reads:

J(w; X ) = 1 N N i=1 x (i) -x(i) 2 , (3)
where x = g -1 (z, w) denotes the reconstruction of the sample x, calculated through the inverse transform g -1 : Z ∈ R m → X ∈ R M . In the general case, additional parameters may be introduced in g -1 , or the inverse transform may not exist at all (see e.g. [START_REF] Kwok | The pre-image problem in kernel methods[END_REF]).

For a detailed description of the specific DR methods used in this paper to showcase the proposed methodology, namely principal component analysis (PCA) and kernel PCA, the reader is referred to Section 4.

Surrogate Modelling

In the context of UQ, the physical or computational model of a system can be seen as a black-box that performs the mapping:

Y = M(X), (4)
where X is a random vector that parametrises the variability of the input parameters (e.g. through a joint probability density function) and Y is the corresponding random vector of model responses. One of the main applications of UQ is to propagate the uncertainties from X to Y through the model M. Direct methods based on Monte-Carlo simulation may require that the computational model is run several thousands of times for different realisations x of the input random vector X. However, most models that are used in applied sciences and engineering (e.g. high-resolution finite element models) can have high computational costs per model run. As a consequence, they cannot be used directly. To alleviate the associated computational burden, surrogate models have become a staple tool in all types of uncertainty quantification applications.

A surrogate model M is a computationally inexpensive approximation of the true model of the form:

M(X) = M(X; θ) + , ( 5
)
where θ is a set of parameters that characterise the surrogate model and refers to an error term. The parameters θ are inferred (typically through some form of optimisation process) from a limited set of runs of the original model X = x (1) , . . . , x (N ) , called the experimental design. As an example, θ denotes the set of coefficients in the case of a truncated polynomial chaos expansion, or the set of parameters of both the trend and the covariance kernel in case of Gaussian process modelling. Throughout the rest of the paper, the output of the model M is considered scalar, i.e. y = M(x) ∈ R.

Arguably the most well-known accuracy measure for most surrogates is the relative generalisation error ε gen that reads:

ε gen = E Y -M(X; θ) 2 /Var [Y ] . (6)
This error measure (or, more precisely, one of its estimators) is also the ideal objective function for the optimisation process involved in the calibration of the surrogate parameters θ. In practical situations, however, it is not possible to calculate ε gen analytically. An estimator ε gen of this error can be computed by comparing the true and surrogate model responses evaluated at a sufficiently large validation set X v = x (1) , . . . , x (Nv) of size N v :

ε gen = Nv i=1 M(x (i) ) -M(x (i) ) 2 Nv i=1 M(x (i) ) -µ y 2 , (7)
where

µ y = 1 N Nv i=1 M(x (i)
) is the sample mean of the validation set responses and M(x (i) ) is used in place of M(x (i) ; θ) to simplify the notation.

In data-driven applications, or when the computational model is expensive to evaluate, only a single set S def = {X , Y} is available. The entire set is therefore used for calculating the surrogate parameters. Estimating the generalisation error by means of Eq. ( 7) on the same set, however, corresponds to computing the so-called empirical error, which is prone to underestimate drastically the true generalisation error, due to the overfitting phenomenon.

In such cases, a fair approximation of ε gen can be obtained by means of cross-validation (CV) techniques (see e.g. [START_REF] Hastie | The elements of statistical learning: Data mining, inference and prediction[END_REF]). In k-fold CV, S is randomly partitioned into k mutually exclusive and collectively exhaustive sets S i of approximately equal size:

S i ∩ S j = ∅ , ∀(i, j) ∈ {1, . . . , k} 2 and k i=1 S i = S. (8)
The k-fold cross-validation error ε CV reads:

ε CV = k i=1 x∈Si M(x) -M S\Si (x) 2 x∈S (M(x) -µ y ) 2 , (9)
where M S S\Si denotes the surrogate model that is calculated using S excluding S i . The bias of the generalisation error estimator is expected to be minimal in the extreme case of leave-one-out (LOO) cross-validation [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF], which corresponds to N -fold cross validation. The LOO error ε LOO is calculated as in Eq. ( 9) after substituting the set S i by the singleton x (i) (i.e. k = N ):

ε LOO = N i=1 M(x (i) ) -M \i (x (i) ) 2 N i=1 M(x (i) ) -µ y 2 , (10)
where the term M \i (x (i) ), denotes the surrogate built from the set S\ x (i) , evaluated at

x (i) . The calculation of ε LOO can be computationally expensive, because it requires the evaluation of N surrogates, but it does not require any additional run of the full computational model. For Gaussian process modelling and polynomial chaos expansions, computational shortcuts are available to alleviate such costs (e.g. [START_REF] Dubrule | Cross validation of Kriging in a unique neighborhood[END_REF]; [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]), in the sense that ε LOO in Eq. ( 10) is evaluated from a single surrogate model M calculated from the full data set S.

As a final step in the surrogate modelling procedure, the set of parameters θ of the surrogate model are optimised w.r.t. to one of the generalisation error measures in Eq. ( 9)

or Eq. ( 10) directly, based on the available samples in S, i.e.:

θ = arg min D θ gen (θ; S), (11)
where θ denotes the optimal set of parameters, D θ the feasible domain of parameters and gen refers to the chosen estimator of gen . An important aspect of this optimisation step for many types of recent surrogates is that the number of parameters θ scales with the number of input variables. Therefore, surrogates tend to suffer from the curse of dimensionality in two distinct ways: higher dimensional optimisation and underdetermination. Higher dimensional optimisation is linked to a complex objective-function topology, and is therefore prone to convergence to low-performing local minima. In general it requires global optimisation algorithms, such as genetic algorithms, covariance matrix adaptation, or differential evolution [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF][START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (cma-es)[END_REF][START_REF] Yang | Differential evolution for high-dimensional function optimization[END_REF]. Underdetermination leads the solutions to the minimisation problem to be non-unique due to the lack of constraining data. In other words, surrogate models with more parameters require in general a larger experimental design or sparse minimisation techniques to avoid overfitting.

3 The proposed DRSM approach

Introduction

Consider now the experimental design S = {X , Y} introduced above, and assume that it is the only available information about the problem under investigation. Moreover, the dimensionality of the input space is high, i.e.

x (i) ∈ R M , i = 1, . . . , N where M is large, say O 10 2-4 .
The goal is to calculate a surrogate model that serves as an approximation of the real model solely based on the available samples. This is a key ingredient for subsequent analyses in the context of uncertainty quantification.

To distinguish between various computational schemes, we denote from now on by M|X , Y a surrogate model whose parameters θ are calculated from the experimental design X and associated model response Y. Due to the high input dimensionality, a surrogate M|X , Y may lead to poor generalisation performance or it may not even be computationally tractable.

To reduce the dimensionality, the class of DR methods was introduced in Section 2.1. A DR transformation, expressed by Z = g(X ; w), can provide a compressed experimental design,

i.e. z (i) ∈ R m , i = 1, . . . , N with m M . The surrogate M|Z, Y becomes tractable if m
is sufficiently small. The potential of M|Z, Y to achieve satisfactory generalisation performance depends on (i) the learning capacity of the surrogate itself and (ii) the assumption that the input-output map x → y can be sufficiently well approximated by a smaller set of features via the transformation g(•). This discussion focuses on the latter and assumes that the learning capacity of the surrogate is adequate. In case of unstructured inputs, the importance of each input variable may vary depending on the output of interest. In case of structured inputs, there is typically high correlation between the input components. Hence, in both families of problems a low-dimensional representation may often approximate well the input-output map.

Traditional DR approaches are focused on the discovery of the input manifold and not the input-output manifold. Performing an input compression without taking into account the associated output values may lead to a highly complex input-output map that is difficult to surrogate. In the DRSM (dimensionality reduction for surrogate modelling) approach proposed in this paper, we capitalise on this claim to try and find an optimal input compression scheme w.r.t. the generalisation performance of M|Z, Y.

A nested optimisation problem

The goal of DRSM is to optimise the parameters w of the compression scheme so that the auxiliary variables z = g(x; w) are suitable to achieve an overall accurate surrogate. The general formulation of this problem reads:

w, θ = arg min w∈Dw, θ∈D θ M(•), M (g(•; w), θ) , (12)
where denotes the objective function (a.k.a. loss function) that quantifies the generalisation performance of the surrogate. In practice, if a validation set is available, corresponds to a generalisation error estimator like the one in Eq. ( 7). In the absence of a validation set, then either the LOO estimator in Eq. ( 10) or its k-fold CV counterpart in Eq. ( 9) are used instead. In the following, it is assumed that a validation set is not available and the generalisation error is estimated by the LOO error, hence is substituted by the ε LOO expression in Eq. ( 10).

The proposed approach for solving Eq. ( 12), is related to the concept of block-coordinate descent [START_REF] Bertsekas | Nonlinear programming[END_REF]. During optimisation, the parameters w and θ are updated in an alternating fashion. One of the main reasons for this choice is that the optimisation steps of both DR and SM techniques are often tuned ad-hoc to optimise their performance. Examples include sparse linear regression for polynomial chaos expansions [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF], or quadratic programming for support vector machines for regression [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. A single joint optimisation, albeit potentially yielding accurate results, would require the definition of complex constraints on the different sets of parameters w and θ. Therefore, the problem in Eq. ( 12) is expressed as a nested-optimisation problem. The outer loop optimisation reads:

w = arg min w∈Dw ε LOO (w; θ(w), X , Y), (13)
where ε LOO denotes the LOO error (Eq. ( 10)) of the surrogate M(z; w, X , Y) evaluated at {X , Y} and θ(w) denotes the optimal parameters of M for that particular w value. The term θ(w) is calculated by solving the inner loop optimisation problem:

θ = arg min θ∈D θ ε LOO (θ; w, X , Y). ( 14
)
The nested optimisation approach to DRSM comes with costs and benefits. On the one hand, each objective function evaluation of the outer-loop optimisation becomes increasingly costly w.r.t. the number of samples in the experimental design and the complexity of the surrogate model. On the other hand, the search space in each optimisation step can be significantly smaller, compared to the joint approach, due to the reduced number of optimisation variables. Moreover, this nested optimisation approach enables DRSM to be entirely non-intrusive. Off-the-shelf well-known surrogate modelling methods can be used to solve Eq. ( 14).

Proxy surrogate models for the inner optimisation

Albeit non-intrusive and having a relatively low dimension, the inner optimisation in Eq. ( 14) is in general the driving cost of DRSM. Indeed, calculating the parameters of a single highresolution modern surrogate may require anywhere between a few seconds and several minutes. To reduce the related computational cost, it is often possible to solve proxy surrogate problems, i.e. using simplified surrogates that, while not being as accurate as their full counterparts, are easier to parametrise. A simple example would be to prematurely stop the optimisation in the inner loop in Eq. ( 14), or to use isotropic kernels for kernel-based surrogates such as Kriging or support vector machines instead of their more accurate, but costly to train, anisotropic counterparts. Once the outer loop optimisation completes on the proxy surrogate, thus identifying the quasi-optimal DR parameters w, a single high-accuracy surrogate is then computed on the compressed experimental design {Z = g(X ; w), Y}. Further discussion on this topic can be found in Sections 4.3 and 4.3.1.

Selected compression and surrogate modelling techniques used in this paper

Due to the non-intrusiveness in the design of the DRSM method proposed in Section 3, no specific dimensionality reduction or surrogate modelling technique has been introduced yet. In the following section, two well-known dimensionality reduction (namely principal component analysis and kernel-principal component analysis) and two surrogate modelling techniques (Kriging and polynomial chaos expansions) are introduced to showcase the DRSM methodology on several example applications in Section 5. Only the main concept and notation is reminded so that the paper is self-consistent.

Principal component analysis

Principal Component Analysis (PCA) is a dimensionality reduction technique that aims at calculating a linear basis of X with reduced dimensionality that preserves the sample variance [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF]. Given a sample of the input random vector X = x (1) , . . . , x (N ) , the PCA algorithm is based on the eigen-decomposition of the sample covariance matrix C:

C = 1 N X X , (15)
of the form:

Cv (i) = λ (i) v (i) , i = 1, . . . , M (16)
where X denotes the centred (zero mean) experimental design, λ (i) denotes each eigenvalue of C and v (i) the corresponding eigenvector. The dimensionality reduction transformation reads:

Z = X V (17)
where V is the M ×m collection of the m eigenvectors of C with maximal eigenvalues. Those eigenvectors are called the principal components because they correspond to the reduced basis of X with maximal variance. Based on the general DR perspective that was presented in Section 2.1, PCA is a linear transformation of the form Z = g(X ; w), where the only parameter to be selected is the dimension m of the reduced space, i.e. w = m.

Kernel principal component analysis

Kernel PCA (KPCA) is the reformulation of PCA in a high-dimensional space that is constructed using a kernel function [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. A kernel function applied on two elements x (i) , x (j) ∈ D x has the following form:

κ x (i) , x (j) = Φ x (i) • Φ x (j) (18)
where Φ(•) is a function that performs the mapping Φ : D x → H and H is known as the feature space. Based on Eq. ( 18), the so-called kernel trick is applied, which refers to the observation that, if the access to H only takes place through inner products, then there is no need to explicitly define Φ(•). The result of the inner product can be directly calculated using

κ(•, •).
Kernel PCA is a non-linear extension of PCA where the kernel trick is used to perform PCA in H. The principal components in H are obtained from the eigen-decomposition of the sample covariance matrix C H , analogously to the PCA case in Eq. ( 15).

However, in KPCA the eigen-decomposition problem:

C H v (i) = λ i v (i) , i = 1, . . . , N (19)
is intractable, since C H cannot in general be computed (H might even be infinitely dimensional). This problem is by-passed by observing that each eigenvector belongs to the span of the samples Φ x (1) , . . . , Φ x (N ) , therefore scalar coefficients α

(i)
k exist, such that each eigenvector v (i) can be expressed as the following linear combination [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]:

v (i) = N k=1 α (i) k Φ x (k) , i = 1, . . . , N. (20)
Based on Eq. ( 20) it can be shown that the eigen-decomposition problem in Eq. ( 19) can be cast as:

Kα (i) = λ (i) α (i) , i = 1, . . . , N (21)
where K is the kernel matrix with elements:

K ij = κ x (i) , x (j) . ( 22
)
As for the case of PCA, Z is calculated by projecting X on the m principal axes v (i) , i = 1, . . . , m corresponding to the m largest eigenvalues. [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] showed that Z can be directly computed based only on the values of the eigenvector expansion coefficients α

(i)
k and the kernel matrix K. The k-th component of the i-th sample of Z, denoted by z (i) k is given by;

z (i) k = Φ x (i) T v (k) = N j=1 α (j) k κ x (i) , x (j) (23)
The key ingredient of KPCA is arguably the kernel function κ. In this paper two kernels are considered, namely the polynomial kernel:

κ(x, x ; w) = w 1 x T x + w 2 w3 , w 1 > 0, w 2 ≥ 0, w 3 ∈ N, (24)
and the Gaussian kernel:

κ(x, x ; w) = exp - 1 2 M k=1 1 w 2 k (x k -x k ) 2 , w k > 0 , k = 1, . . . , M. (25)
A special case of the Gaussian kernel is the isotropic Gaussian kernel (also known as radial basis function) that simply assumes the same parameter value w k for all components of x.

Note that KPCA using a polynomial kernel with parameters w 1 = 1, w 2 = 0 and w 3 = 1 is identical to PCA, since Φ(x) = x. A discussion on the equivalence between PCA and KPCA with linear kernel (w 3 = 1) for arbitrary values of w 1 , w 2 can be found in Appendix A. From

Eq. ( 23) it follows that Z can be expressed as Z = g(X ; w) where w encompasses both the kernel parameters and the reduced space dimension m.

In the context of unsupervised learning, two methods to infer the values of w from X are considered. The distance preservation method aims at optimising w in such a way that the Euclidean distances between the samples are preserved between the original and the feature space [START_REF] Weinberger | Learning a kernel matrix for nonlinear dimensionality reduction[END_REF]. This is expressed by the following objective function:

J dist (w; X ) = N i,j=1 (d ij -δ ij ) 2 (26)
where

d ij = x (i) -x (j) (27)
and

δ ij = Φ(x (i) , w) -Φ(x (j) , w) . ( 28
)
By expanding the norm expression in Eq. ( 28) it is straightforward to show that:

δ ij = K ii + K jj -2K ij , (29)
hence the value of δ ij is readily available from the kernel matrix K.

The reconstruction error -based method aims at optimising w in such a way that the socalled pre-image, x = g -1 (z, w ), of z = g(x, w) approximates x as close as possible [START_REF] Alam | Hyperparameter selection in kernel principal component analysis[END_REF]. This is expressed by the following objective function:

J recon (w; X ) = 1 N N i=1 x (i) -x(i) 2 (30)
In contrast to PCA, calculating x is non-trivial, an issue that is known as the pre-image problem (see e.g. [START_REF] Kwok | The pre-image problem in kernel methods[END_REF]). The approach for dealing with this problem is the one adopted by the popular python package scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF], which is based on [START_REF] Weston | Learning to find pre-images[END_REF]. After performing the KPCA transform X → Z, the (non-unique) pre-image of a new point z is computed by kernel-ridge regression using a new kernel function κ pre :

x = β T l(z), (31)
where:

(z) = κ pre (z, z (j) ), j = 1, . . . , N , (32)
and β are the kernel-ridge regression coefficients. They are calculated as follows:

β = (L + rI N ) -1 X L ij = κ pre z (i) , z (j) , i, j = 1, . . . , N (33)
where r is a regularisation parameter and I N is the N -dimensional identity matrix. In [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] and in this paper, we use for simplicity the same kernel for the pre-image problem as for KPCA, i.e.

κ pre (•, •) is chosen equal to κ (•, •).
Note that, in the unsupervised learning literature, the reduced space dimension, m, is typically not part of w, i.e. only the kernel parameters are considered when minimising the objective function in Eq. ( 26) or Eq. ( 30).

Kriging

Kriging, a.k.a. Gaussian process modelling, is a surrogate modelling technique which assumes that the true model response is a realisation of a Gaussian process described by the following equation [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]:

M(x) = β f (x) + σ 2 Z(x) (34)
where β f (x) is the mean value of the Gaussian process, also called trend, σ 2 is the Gaussian process variance and Z(x) is a zero-mean, unit-variance Gaussian process. This process is fully characterised by the auto-correlation function between two sample points R(x, x ; θ).

The hyperparameters θ associated with the correlation function R(•; θ) are typically unknown and need to be estimated from the available observations. Various correlation functions can be found in the literature [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF], including the linear, exponential, Gaussian (a.k.a. squared exponential ) and Matérn functions. In this paper the separable Matérn correlation family is chosen:

R (|x -x | ; l, ν) = M i=1 1 2 ν-1 Γ(ν) √ 2ν |x i -x i | l i ν κ ν √ 2ν |x i -x i | l i , (35)
where x, x are two samples in the input space D x , l = {l i > 0, i = 1, . . . , M } are the scale parameters (also called correlation lengths), ν ≥ 1/2 is the shape parameter, Γ(•) is the Euler Gamma function and κ ν (•) is the modified Bessel function of the second kind (a.k.a.

Bessel function of the third kind). The values ν = 3/2 and ν = 5/2 of the shape parameter are commonly used in the literature. The isotropic variant of the Matérn correlation family assumes a fixed correlation length value l in Eq. ( 35) over all M input variables.

Regarding the trend part β f (x) in Eq. ( 34), the general formulation of universal Kriging is adopted, which assumes that the trend is composed of a linear combination of P preselected functions {f i (x), i = 1, . . . , P }, i.e.:

β f (x) = P i=1 β i f i (x), ( 36
)
where β i is the trend coefficient of each function.

The Gaussian assumption states that the vector formed by the true model responses, y and the prediction, Y (x), at a new point x, has a joint Gaussian distribution defined by: 

 Y (x) y   ∼ N N +1     f (x)β F β   , σ 2   1 r (x) r(x) R     ( 37
)
where F is the information matrix of generic terms:

F ij = f j (x (i) ) , i = 1, . . . , N, j = 1, . . . , P, (38)
r(x) is the vector of cross-correlations between the prediction point x and each one of the observations whose terms read:

r i (x) = R(x, x (i) ; θ), i = 1, . . . , N. (39)
R is the correlation matrix given by:

R ij = R(x (i) , x (j) ; θ), i, j = 1, . . . , N. (40)
The mean and variance of the Gaussian random variate Y (x) (a.k.a. mean and variance of the Kriging predictor) can be calculated based on the best linear unbiased predictor (BLUP)

from [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]:

µ Y (x) = f (x) β + r(x) R -1 (y -F β) , (41)
σ 2 Y (x) = σ 2 1 -r (x)R -1 r(x) + u (x)(F R -1 F ) -1 u(x) (42)
where:

β = F R -1 F -1 F R -1 y (43)
is the generalised least-squares estimate of the underlying regression problem and

u(x) = F R -1 r(x) -f (x). ( 44
)
The mean response in Eq. ( 41) is considered as the output of a Kriging surrogate, i.e.

M(x) = µ Y (x).

It is important to note that the Kriging model interpolates the data, i.e.:

µ Y (x) = M(x), σ 2 Y (x) = 0, ∀ x ∈ X (45)
The equations that were derived for the best linear unbiased Kriging predictor assumed that the covariance function σ 2 R(•; θ) is known. In practice however, the family and other properties of the correlation function need to be selected a priori. The hyperparameters θ, the regression coefficients β and the variance σ 2 need to be estimated based on the available experimental design.

The optimal estimates of the correlation parameters θ are determined by minimising the generalisation error of the Kriging surrogate, based on the leave-one-out cross-validation error [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF][START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model misspecifications[END_REF]:

θ CV = arg min D θ K i=1 M(x (i) ) -µ Y ,(-i) (x (i) ) 2 , ( 46
)
where µ Y ,(-i) (x (i) ) corresponds to the mean value of a Kriging predictor that was built from the samples X \ x (i) , y (i) , evaluated at x (i) . The computational cost for calculating the terms µ Y ,(-i) (x (i) ) can be significantly reduced as shown in [START_REF] Dubrule | Cross validation of Kriging in a unique neighborhood[END_REF]. First, the following matrix inversion is performed:

B =   σ 2 R F F T 0   -1 . (47)
Then µ Y ,(-i) is calculated as follows:

µ Y ,(-i) = - N j=1,j =i B ij B ii y (j) . (48)
In this work we use cross-validation for estimating the correlation parameters instead of the maximum likelihood method [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]). This is motivated by the comparative study in [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model misspecifications[END_REF] between maximum likelihood (ML) and CV estimation methods.

The CV method is expected to perform better in cases that the correlation family of the Kriging surrogate is not identical to the one of the true model. This is typically the case in practice and in the application examples in Section 5.

Determining the optimal parameters θ CV in Eq. ( 46) leads to a complex multi-dimensional optimisation problem. Common optimisation algorithms employed to solve Eq. ( 46) can be cast into two categories: local and global. Local methods are usually gradient-based, such as the BFGS algorithm [START_REF] Bazaraa | Nonlinear programming: theory and algorithms[END_REF], and search locally in the vicinity of the starting point. This makes them prone to get stuck at local minima, although they can be computationally efficient due to the use of gradients. Global methods such as genetic algorithms [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF] do not rely on local information such as the gradient. They seek the global minimum by various adaptive resampling strategies within a bounded domain. This often leads to considerably more objective function evaluations compared to local methods.

As mentioned in Section 3.3, to alleviate the computational costs in the inner loop optimisation in Eq. ( 14), an inexpensive-to-calibrate Kriging surrogate is built. To this end, the isotropic version of the Matérn correlation family is used, combined with low computational budget optimisation of the correlation parameters. For calculating the final, high-accuracy, Kriging surrogate an optimisation with high-computational budget is performed instead, combined with the use of an anisotropic correlation family. The introduction of anisotropy is expected to improve the generalisation performance the metamodel, as shown for instance in the study by [START_REF] Moustapha | Comparative study of Kriging and support vector regression for structural engineering applications[END_REF].

Polynomial chaos expansions

Polynomial chaos expansions represent a different class of surrogate models that has seen widespread use in the context of uncertainty quantification due to their flexibility and efficiency. Consider that X ∈ R M is a random vector with independent components described by the joint PDF f X and that the model output Y in Eq. ( 4) has finite variance. Then the polynomial chaos expansion of M(X) is given by:

Y = M(X) = α∈N M θ α Ψ α (X) (49)
where the Ψ α (X) are multivariate polynomials orthonormal with respect to f X , α ∈ N M is a multi-index that identifies the components of the multivariate polynomials Ψ α and the θ α ∈ R are the corresponding coefficients.

In practice, the series in Eq. ( 49) is truncated to a finite sum, by introducing the truncated polynomial chaos expansion:

M(X) ≈ M(X) = α∈A θ α Ψ α (X) ≡ θ Ψ(x) (50)
where A ⊂ N M is the set of selected multi-indices of multivariate polynomials. A typical truncation scheme consists in selecting multivariate polynomials up to a total degree p, i.e.

A = α ∈ N M : α 1 ≤ p , with α 1 = M i=1 α i .
The corresponding number of terms in the truncated series rapidly increases with M , giving rise to the "curse of dimensionality".

Other truncation strategies effective in higher dimension are discussed, e.g., in [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]; [START_REF] Jakeman | Enhancing 1 -minimization estimates of polynomial chaos expansions using basis selection[END_REF].

The polynomial basis Ψ α (X) in Eq. ( 50) is traditionally built starting from a set of univariate orthonormal polynomials φ (i) k (x i ) which satisfy:

φ (i) j (x i ), φ (i) k (x i ) def = D X i φ (i) j (x i )φ (i) k (x i )f Xi (x i )dx i = δ jk ( 51
)
where i identifies the input variable w.r.t. which they are orthogonal, as well as the corresponding polynomial family, j and k the corresponding polynomial degree, f Xi (x i ) is the i th -input marginal distribution and δ jk is the Kronecker symbol. Note that this definition of inner product can be interpreted as the expectation value of the product of the multiplicands. The multivariate polynomials Ψ α (X) are then assembled as the tensor product of their univariate counterparts:

Ψ α (x) def = M i=1 φ (i) αi (x i ) (52)
For standard distributions, such as uniform, Gaussian, gamma, beta, the associated families of orthogonal polynomials are well-known [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. Orthogonal polynomials can be constructed numerically w.r.t. any distribution (including non-parametric ones like those obtained by kernel density smoothing) by means of Gram-Schmidt orthonormalisation (a.k.a. Stieltjes procedure for polynomials [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]).

The expansion coefficients θ = θ α , α ∈ A ⊂ N M in Eq. ( 50) are calculated by minimising the expectation of least-squares residual [START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF]:

θ = arg min E θ T Ψ(X) -M(X) 2 . ( 53
)
In the context of DRSM, the set of input parameters w for a PCE surrogate consists in w = {p, θ}, i.e. the maximal degree of the truncated expansion and the associated coefficients.

Due to the quadratic programming nature of the minimisation in Eq. ( 53) and the linearity of PCE (see Eq. ( 50)), we adopt the adaptive sparse-linear regression based on least angle regression first introduced by [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF].

As for the case of Kriging, the LOO error (see Eq. ( 10)) is analytically available from the expansion coefficients [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF]:

ε LOO = N i=1 M(x (i) ) -M P C (x (i) ) 1 -h i 2 N i=1 M(x (i) ) -µ Y 2 , ( 54
)
where h i is the i th component of the vector given by:

h = diag A(A T A) -1 A T , ( 55
)
and A is the experimental matrix with entries A ij = Ψ j x (i) ) .

To calculate the proxy PCE surrogates used during the DRSM optimisation phase (see Section 3.3), the input variables in z, are assumed uniformly distributed and independent.

The PCE coefficients are computed by solving Eq. ( 53) using the ordinary least squares method [START_REF] Berveiller | Stochastic finite elements: a non intrusive approach by regression[END_REF]. To calculate the PCE coefficients of the final, highaccuracy, surrogate M(g(x; w)), the distributions of the input variables are fitted using kernel-smoothing, while retaining the independence assumption, motivated by the results in [START_REF] Torre | Data-driven polynomial chaos expansion for machine learning regression[END_REF]. In addition, a sparse solution is obtained by solving the optimisation problem in Eq. ( 53) using least angle regression [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] instead of ordinary least squares.

Applications

The performance of DRSM is evaluated on the following applications: (i) an artificial analytic function with 20 unstructured inputs and approximately known intrinsic dimension, (ii) a realistic electrical engineering model with 80 unstructured inputs and unknown intrinsic dimension and, (iii) a heat diffusion model with 16, 000 structured inputs and unknown intrinsic dimension.

For each example, DRSM is applied using KPCA for compression together with Kriging or polynomial chaos expansions for surrogate modelling. The surrogate performance is then compared, in terms of generalisation error, to the sequential application of unsupervised dimensionality reduction followed by surrogate modelling. To improve readability, various details regarding the implementation of the optimisation algorithms and the surrogate models calibration are omitted from the main text and given in Appendix B instead.

All the surrogate modelling techniques were deployed with the Matlab-based uncertainty quantification software UQLab (Marelli andSudret, 2014, 2018;[START_REF] Lataniotis | The Gaussian process modelling module in UQLab[END_REF].

Sobol' function

The Sobol' function (also known as g-function) is a commonly used benchmark function in the context of uncertainty quantification. It reads:

Y = M i=1 |4X i -2| + c i 1 + c i , ( 56
)
where X = {X 1 , . . . , X M } are independent random variables uniformly distributed in the interval [0, 1] and c = {c 1 , . . . , c M } T are non-negative constants. In this application, we chose M = 20 and the constants c given by Konakli and Sudret (2016a); [START_REF] Kersaudy | A new surrogate modeling technique combining Kriging and polynomial chaos expansions -Application to uncertainty analysis in computational dosimetry[END_REF]: 2, 5, 10, 20, 50, 100, 500, 500, . . . , 500} T .

c = {1,
(57)

It is straightforward to see that the effect of each input variable X i to the output Y is inversely proportional to the value of c i . In other words, a small (resp. large) value of c i results in a high (resp. low) contribution of X i to the value of Y i . For the given values of the constants c, one would expect that, roughly, the first 4 to 6 variables can provide a compressed representation of X with minimal information loss regarding the input-output relationship.

To showcase the performance of DRSM, an experimental design X , consisting of 800 samples, is generated by Latin Hypercube sampling of the input distribution [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. Based on the samples in X and the corresponding model responses Y, several combinations of KPCA, Kriging and PCE are tested within the DRSM framework. An additional set of 10 5 validation samples {X v , Y v } is generated for evaluating the performance of the final surrogates.

The first analysis consists in comparing the generalisation performance as a function of the compressed input dimension m for Kriging and PCE models combined with KPCA with different kernels. Because of the availability of a validation set, the performance of the LOO error estimator in Eq. ( 10) is also assessed by comparing it with the true validation error in Eq. ( 6). Figures 1a and1b show the LOO error estimator of the final surrogate model when using Kriging and PCE, respectively. In each panel the different curves correspond to different KPCA kernels, namely polynomial kernel (Eq. ( 24)) and isotropic (resp. anisotropic) Gaussian (Eq. ( 25)). Figures 1c and1d show the corresponding validation error on the validation set for the same scenarios. At a first glance, it is clear that the top and bottom figures are remarkably similar, both in their trends and in absolute value. Therefore, it is concluded that on this example LOO is a good measure of the generalisation error gen . This is an important observation, because in the general case a validation set is not available, while LOO can always be calculated. Moreover, the intrinsic dimension identified by all the best DR-SM combinations is equal to m = 6, which is a reasonable estimate based on the values of the constants c i in Eq. ( 57).

The DRSM algorithm identifies the anisotropic Gaussian kernel as the best KPCA kernel to be used in conjunction with both Kriging and PCE. However, the performance of PCE is significantly better in terms of generalisation error. The optimal parameters for each case (Kriging and PCE) are highlighted by a black dot in Figure 1, and their numerical values are reported in Table 1. Subsequently, the performance of DRSM is compared against an unsupervised approach, in which dimensionality reduction is carried out first, before applying surrogate modelling.

To facilitate a meaningful comparison between the various methods, the reduced dimension and the optimal KPCA kernel as determined by the first analysis (see Table 1) is used. The results are summarised in Figure 2, while the corresponding list of tested configurations for both DRSM and the sequential DR-SM is given in Table 2.

The experimental design consists of 800 samples. The performance of each method is evaluated in terms of the generalisation error of the final surrogate M(z) evaluated on a validation set {X v , Y v = M(X v )} with 10 5 samples. To evaluate the robustness of the results, this process is repeated 10 times, each corresponding to a different set X , drawn Table 2: Different setups considered for evaluating the final surrogate model performance after using each of them for dimensionality reduction.

Dim. reduction Parameter tuning objective Abbreviation

Kernel PCA ε LOO of Kriging (KG) or PCE surrogate (Eq. ( 13 The DRSM approach consistently shows superior performance compared to the unsupervised approaches. This performance improvement becomes more apparent in the case of PCE surrogate modelling, where the average validation error over the 10 repetitions is reduced by almost two orders of magnitude compared to the other methods.

Due to the analytical nature of the model under consideration, we further evaluate the DRSM-based input compression by means of how the most important input variables are mapped to the reduced space. We adopt the total Sobol' sensitivity indices as a rigorous measure of the importance of each input variable. Sobol' sensitivity analysis is a form of global sensitivity analysis based on decomposing the variance of the model output into contributions that can be directly attributed to inputs or sets of inputs [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. The total Sobol' sensitivity index of an input variable X i , denoted by S T ot i ∈ [0, 1], quantifies the total effect of X i on the variance of Y . In this particular example, the total Sobol' indices can be analytically derived (Saltelli et al., 2000). Their values are shown for reference in Figure 3a.

It is clear from Eq. ( 56) and Eq. ( 57) that all 20 input variables contribute to the output variability, i.e. the intrinsic dimension of the problem is 20. However, the contribution of each input component quickly diminishes with larger values of c i (see Figure 3a in which the values of the 20 total Sobol' indices are plotted, in logarithmic scale, as horizontal bars).

Compressing the inputs in this problem is expected to lead to a mapping where those first few input components have the largest contribution.

In Figure 3 the features in the reduced space Z are compared against the original inputs X. The rationale behind this heuristic analysis is simple: if the features obtained by DRSM are correctly identified, they should depend mostly on the same variables identified as important in the Sobol' analysis in Figure 3a. A simple measure of dependence between the reduced space components {z i , i = 1, . . . , m} and the initial input space components correlates strongly with a specific x i , (ii) the z i 's correlate with the m "most important"

{x i , i = 1, . . . ,
x i 's, and, (iii) the larger m value leads to the discovery of a new input z i that correlates with the next "most important" component of x. As in the previous section, the goal of the first analysis is to determine the generalisation performance of the DRSM surrogate as a function of the reduced space dimension m when KPCA is combined with either Kriging or PCE. In addition, the accuracy of the LOO error in Eq. ( 10) is compared to the validation error in Eq. ( 7). or PCE), evaluated on {X , Y}, whereas Figures 5c and5d show the validation error of the surrogate, evaluated on {X v , Y v }. In each panel, each curve corresponds to a different KPCA kernel, namely anisotropic or isotropic Gaussian, and polynomial. Finally, the optimal configuration for each SM method is illustrated by a black dot. Similarly to the Sobol' function, the use of an anisotropic kernel in KPCA results in significantly reduced generalisation error.

Electrical resistor network

Indeed this is expected from a physical standpoint. The effect of the resistors on the voltage V will decay with distance (in terms of the number of preceding resistors) from V , which implies anisotropy in terms of the effect of each input variable to the output. As in the previous application example, the LOO error in Figures 5a and5b provides a reliable proxy of the generalisation error in Figures 5c and5d and the same optimal parameters are identified w.r.t. the two error measures. The optimal DRSM configuration for each surrogate model is given in Table 3. Next, the performance of DRSM is compared to unsupervised approaches considering the setups in Table 2. The results of this comparative study are given in Figure 6 using box plots. They are obtained by the repeated random selection of 500 samples from the available 1, 000, leading to 10 separate surrogate models for each case. The performance of each method is determined by means of the ε gen of the final surrogate M(z) evaluated on the validation set {X v , Y v = M(X v )}, that corresponds to the remaining 500 samples of each split. Hence, each box-plot provides summary statistics of the validation error over the different splits. Each of the setups is tested both for Kriging (Figure 6a) and PCE surrogates (Figure 6b). In this application example the DRSM-based surrogates outperform the others by several orders of magnitude in both cases (Kriging, PCE). This highlights the difference between the unsupervised and supervised compression: compressing the input using only the information in X appears inefficient when followed by surrogate modelling.

2D heat diffusion

This last application consists in a 2-dimensional stationary heat diffusion problem. The problem is defined in a square domain, D = [-0.5, 0.5] × [-0.5, 0.5], where the temperature field T (v), v ∈ D is the solution of the elliptic partial differential equation:

-∇ • (d(v)∇T (v)) = 500 I A (v), (58)
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. ( 58),

A corresponds to a square domain (see Figure 7) and I A is the indicator function equal to 1 if v ∈ A and 0 otherwise. The diffusion coefficient d(v) is a lognormal random field defined by: To solve Eq. ( 58), the Gaussian random field g(v) is first discretised using the expansion optimal linear estimation (EOLE) method [START_REF] Li | Optimal discretization of random fields[END_REF]. Consider a grid in D with nodes {v 1 , . . . , v n }. By retaining the first p terms in the EOLE series, g(v) is approximated by:

d(v) = exp (a d + b d g(v)) , (59)
g(v) = p i=1 ξ i √ l (i) φ (i) C vv (v), (60)
where {ξ 1 , . . . , ξ p } are independent standard normal random variables, C vv is a vector i) , i = 1, . . . , n} are the eigenvalues and eigenvectors of the correlation matrix C vv with elements C (i,j) vv = R(v i , v j ) for i, j = 1, . . . , n. In the following analysis the Gaussian random field realisations are computed using p = 30 terms in the EOLE series in Eq. ( 60), which allows to represent 93.69% of the variance of the original field.

with elements C (k) vv = R(v, v k ) for k = 1, . . . , n and { l (i) , φ (
The underlying deterministic problem is solved with an in-house finite-element analysis code developed in Matlab. The mesh shown in Figure 7a consists of 16, 000 triangular T3 elements. Figure 7b shows a realisation of the diffusion coefficient random field which corresponds to the input of the model. The corresponding model output, shown in Figure 7c, is the mean temperature in the highlighted square region B. Each realisation of the diffusion coefficient random field is discretised over the mesh in Figure 7a. In the following analysis, the system is treated as a black-box, with the discretised heat diffusion coefficient as a high-dimensional input (M = 16, 000) and the average temperature in square B as the scalar model output. A single set of 500 experimental design samples and model responses is available. This example mimics a realistic scenario in which various maps of spatially varying parameters measured on a regular grid, are input to a computational model that analyses some performance of the system.

dimensionality.

A similar convergence behaviour is observed between Kriging-and PCE-based DRSM.

The corresponding optimal parameter values are highlighted in Figure 8 and their numerical values are reported in Table 4. The linear polynomial kernel performs best in both cases and leads to the same reduced space dimension m = 20. This significantly low dimension can be interpreted by Eq. ( 60). The heat diffusion coefficient, although 16, 000-dimensional, is a non-linear combination of p independent standard normal random variables. Moreover, the LOO and validation error curves show similar behaviour both in terms of their trend and their absolute value. Hence, the LOO error served as a reliable proxy of the validation error, as was observed in the previous application examples too.

In the subsequent analysis we compare the performance of the DRSM approach against other sequential approaches listed in Table 2. To test each setup, we repeat the calculation posed a novel approach for effectively combining dimensionality reduction with surrogate modelling, called DRSM. DRSM consists of three steps: (i) the DR and SM parameters are calculated by solving a nested optimisation problem, where only low-accuracy surrogates are considered to reduce the associated computational cost, (ii) the optimal configuration parameters, including the dimension of the reduced space, are empirically estimated based on the surrogate model performance, and, (iii) a final high-accuracy surrogate is calculated using the optimal values of all the aforementioned parameters.

The performance of DRSM was compared on three different benchmark problems of varying complexity against the classical approach of tuning the dimensionality reduction and surrogate modelling parameters sequentially. DRSM consistently showed superior performance compared to the others in all the benchmark applications.

The novelty of the proposed methodology lies in its non-intrusive way of combining dimensionality reduction and surrogate modelling. This allows for the combination of various techniques without the need of tweaking the dedicated optimisation algorithms on which each of them capitalises. A practical implication of the non-intrusiveness of DRSM is that off-the-shelf surrogate modelling methods (or even software) with sophisticated calibration algorithms can be directly used within this framework.

The focus was given to data-driven scenarios where only a limited set of observations and model responses is available. We demonstrated that the leave-one-out cross-validation error of the surrogate models can serve as a reliable proxy for estimating the generalisation error in order to tune the DR parameters, but also to assess the overall accuracy of the resulting surrogate.

It is noteworthy to mention that in application-driven scenarios where the goal is to obtain a surrogate with optimal performance (regardless of its type) for that specific problem, the proposed approach could be extended in a way that the surrogate type itself is included as one of the parameters that DRSM needs to optimise. However, special care would need to be given to the error metric used during the DRSM optimisation in this case, because the LOO error estimations by different surrogates may have widely varied levels of bias (see e.g. [START_REF] Tibshirani | A bias correction for the minimum error rate in cross-validation[END_REF]).

In future extensions of this work, focus will be given to capitalising on available HPC resources to optimise for different combinations of surrogate models and dimensionality reduction methods. In addition, the cost of training surrogate models increases with the number of available experimental design samples. Therefore, research efforts will also be directed towards dealing with large experimental designs, possibly within a big data framework.

where V H is the collection of the m eigenvectors of C H = cov [Φ(X )] with maximal eigenvalues. Notice that in case of a = 1 and b = 0, from Eqs. ( 61), (64) follows that z = q.

The covariance matrix C H can be expressed as:

C H =      0 . . . 0 . . . a C 0      . ( 65
)
Hence, excluding the eigenvector that corresponds to the zero eigenvalue, it is straightforward to show that

V H =   0 . . . 0 V   . (66)
Based on Eqs. ( 63) and (66), Eq. ( 64) can be written as follows:

q = √ b √ a x   0 . . . 0 V   (67) = √ a x V (68)
= √ a z (from Eq. ( 61))

Therefore, the dimensionality reduction using kernel PCA with a linear kernel provides a scaled version of standard PCA and the constant b has no effect.

B Implementation details

This section provides an extensive list of the configuration parameter values that were used to produce the results in Section 5. Table 5 (resp. Table 6)lists the configuration parameters of Kriging (resp. polynomial chaos expansions) surrogate models. For each surrogate method a distinction is made, in terms of the parameters used, between the proxy (i.e. low computational cost) surrogate and the high-accuracy one. The proxy surrogates were used for solving the nested optimisation problem of DRSM in Eqs. ( 13), ( 14). The same configuration was used to calculate the high-accuracy surrogates regardless of the input compression method (DRSM or disjoint PCA/KPCA).

The parameters of the DRSM-based optimisation are listed in Table 7. Note that the exact same optimisation algorithm and parameters were used for optimising w w.r.t. the KPCA reconstruction and point-wise distance error in the box-plots used to compare the various approaches. The optimisation constraints differ from the ones reported in Table 7 when a polynomial kernel is used in KPCA, as in Eq. ( 24), for improved numerical stability of the solver. On top of the bound constraints reported in the table, that still apply for w 1

Figure 1 :

 1 Figure 1: Error estimates of the DRSM surrogate as a function of the reduced space dimension. Kernel PCA is used with isotropic (resp.anisotropic) Gaussian as well as polynomial kernels.

  random using the Latin Hypercube sampling method. On the left (resp. right) panel, a Kriging (resp. PCE) surrogate is calculated using one of the methods in Table 2. Each box 20 plot in Figure 2 provides summary statistics of the generalisation error that was achieved by each configuration over the 10 repetitions. The central mark indicates the median, and the bottom and top edges of the box indicate the 25 th and 75 th percentiles, respectively. The whiskers extend to the most extreme data points up to 1.5 times the inter-quartile range above or below the box edges. Any sample beyond that range is considered an outlier and plotted as a single point.

Figure 2 :

 2 Figure 2: Sobol' function: estimates of the generalisation error.

Figure 3 :

 3 Figure 3: Sobol' function: visualisation of the sample-based Spearman correlation coefficient (absolute value) between the model inputs X and the reduced space inputs Z.

Figure 4 :

 4 Figure 4: The resistor networks application example

Figure 6 :

 6 Figure 6: Electrical resistor networks: estimates of the generalisation error.

  where g(v) is a Gaussian random field and the parameters a d , b d are such that the mean and standard deviation of d are µ d = 1 and σ d = 0.3 respectively. The random field is characterised by a Gaussian correlation function R(v, v ) = expvv 2 / 2 , with = 0.2. The output of interest is the average temperature in the square domain B within D (see Figure 7).

  Figure 7: 2D heat diffusion problem: illustration of the model input and output.

Figure 8 :

 8 Figure 8: 2D heat diffusion problem: Error estimates of the DRSM surrogate as a function of the reduced space dimension. Kernel PCA is used with isotropic Gaussian and polynomial kernels.

Table 1 :

 1 Sobol' function: optimal DRSM configurations for Kriging-and PCE-based surrogate models

	SM method	KPCA kernel	m ε LOO	ε gen
	Kriging	Anisotropic Gaussian 6 0.0704 0.0830
	PCE	Anisotropic Gaussian 6 0.0096 0.0083

  M } is provided by the metric |ρ (z i , x i )|, where ρ denotes the Spearman correlation coefficient. Figures 3b -3e represent graphically the quantity |ρ (z i , x i )| for the best surrogate identified in Table 1, namely a PCE coupled with KPCA using an anisotropic Gaussian kernel, evaluated on the validation set {X v , Y v }. Each figure corresponds to a different selection of reduced space dimension m. Figure 3 clearly shows that (i) each z i

  The samples are randomly split into 500 pairs {X , Y} used during the DRSM calibration and 500 pairs {X v , Y v } used for validation. Figures 5a and 5b show the LOO error estimator of the final surrogate model (Kriging

Table 3 :

 3 Resistor networks: optimal DRSM configurations for Kriging and PCE surrogate models

	SM method	KPCA kernel	m	ε LOO	ε gen
	Kriging	Anisotropic Gaussian 24 2.000e-04 2.402e-04
	PCE	Anisotropic Gaussian 32 3.621e-05 3.249e-05

Table 4 :

 4 2D diffusion: optimal DRSM configurations for Kriging-and PCE-based surrogate modelsAs in the previous application examples, the goal of the first analysis is to determine the optimal DRSM configuration in terms of the KPCA kernel and the reduced space dimension, as well as test the effectiveness of the LOO error as a proxy of the validation error. In this analysis, the available samples are randomly split into 300 pairs to be used during the DRSM optimisation and 200 pairs to be used for validation. The results are shown in Figure8.

	SM method KPCA kernel m	w (Eq. (24))	ε LOO	ε gen
			w 1		w 2	w 3
	Kriging	Polynomial	20 131.3681 112.0040 1 0.0205 0.0216
	PCE	Polynomial	20 17.5225	15.1853	1 0.0340 0.0356
	Figures 8a and 8b show the LOO error estimator of the final Kriging (resp. PCE) surrogate,
	evaluated on {X , Y}, whereas Figures 8c and 8d show the validation error of the surrogate
	evaluated on {X

v , Y v }. Each curve corresponds to a specific type of KPCA kernel, namely isotropic Gaussian and polynomial, and a specific surrogate, namely Kriging and PCE. We omitted the anisotropic Gaussian kernel for KPCA which is intractable due to the large input

Acknowledgements

Dr John Jakeman (Sandia National Laboratories) is gratefully acknowledged for having provided the data sets used in the electrical resistor networks application example (Section 5.2).

process 10 times. In each case the 500 available samples are split randomly into 300 samples for calculating the surrogate and 200 samples for validation. The optimal KPCA kernel that was determined by DRSM is used in all methods that involve KPCA. Also, for the sake of comparison, the same reduced space dimension m = 20 is assumed for all methods.

The results of this comparative study are given in Figure 9 using box plots to provide summary statistics of the validation error over the different splits of the samples. In case of Kriging surrogate modelling, DRSM consistently provides superior results compared to the other methods. Notice that KPCA with linear kernel is equivalent to PCA on a scaled version of the experimental design with scaling factor √ w 1 (see Appendix A for more details).

The Kriging surrogates, in contrast to the PCE ones, are affected by this scaling. This also explains the performance improvement compared to the case of PCA-based DR. In case of PCE surrogate modelling, the performance improvement gained by DRSM is marginal compared to PCA and KPCA with distance preservation-based tuning of w.

Overall, DRSM consistently provides more accurate or at least comparable results compared to the other approaches. The main difference with a standard UQ setting in which the thermal conductivity is supposed to be sampled from a random field with known properties, is that the proposed DRSM methodology is purely data-driven, i.e. it would be applied identically in a case when the input maps are given without knowing the underlying random process.

Summary and Conclusions

Surrogate modelling is a key ingredient of modern uncertainty quantification. Due to the detrimental effects of high input dimensionality on most recent surrogate modelling techniques, the input space needs to be compressed to make such problems tractable. We pro-A Relationship between PCA and KPCA with linear kernel Consider the PCA-based dimensionality reduction x ∈ R M → z ∈ R m . As discussed in Section 4.1, z is calculated as follows:

where

is the experimental design.

Next, consider the kernel PCA mapping x ∈ R M → q ∈ R m using the linear kernel function:

It is straightforward to show that the following transformation is equivalent to the linear kernel in Eq. (62):

because κ (x, x ) = Φ(x) Φ(x ). A sample q in the reduced space is calculated as follows (see Section 4.2): Max. iterations: 50 for both GA and BFGS and w 2 , the variable w 3 (degree) is constrained to integer values 1 ≤ w 3 ≤ 4 instead. In addition, the following non-linear constraint is included: