Paper deep dive
Adaptive surrogate modeling for high-dimensional spatio-temporal output
Berkcan Kapusuzoglu, Shunsaku Matsumoto, Yoshitomo Miyagi, Daigo Watanabe, Sankaran Mahadevan
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 8/19/2026, 4:31:17 AM
Summary
This paper presents an adaptive surrogate modeling method designed for high-dimensional spatio-temporal outputs in multi-physics systems. The approach utilizes randomized singular value decomposition (rSVD) for dimension reduction to map high-dimensional outputs to a low-dimensional latent space, followed by surrogate model construction in that space. A novel adaptive sampling technique combining exploration and exploitation is introduced to select training points efficiently, minimizing the need for expensive physics-based model runs. The method is validated using thermo-mechanical analysis of a gas turbine engine blade.
Entities (7)
Relation Signals (6)
Berkcan Kapusuzoglu → authored → Adaptive Surrogate Modeling
confidence 99% · Citation: Kapusuzoglu, B., ... Adaptive surrogate modeling for high-dimensional spatio-temporal output.
Adaptive Surrogate Modeling → publishedin → Structural and Multidisciplinary Optimization
confidence 98% · Citation: Kapusuzoglu, B., ... Structural and Multidisciplinary Optimization 65, 290 (2022).
Adaptive Surrogate Modeling → appliedto → Gas Turbine Engine Blade
confidence 95% · Thermo-mechanical analysis of a gas turbine engine blade is used to analyze the effectiveness of the proposed method.
Adaptive Surrogate Modeling → uses → Randomized Singular Value Decomposition
confidence 95% · For dimension reduction, the method of randomized singular value decomposition (rSVD) is used to identify the important features in the output space
Adaptive Surrogate Modeling → handles → Spatio-Temporal Output
confidence 93% · This paper develops an adaptive surrogate modeling method for problems with very high-dimensional spatio-temporal outputs.
Adaptive Surrogate Modeling → employsstrategy → Exploration and Exploitation
confidence 90% · We present a novel adaptive sampling technique that combines exploration and exploitation to improve the surrogate model accuracy
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper develops an adaptive surrogate modeling method for problems with very high-dimensional spatio-temporal outputs. The analysis of spatio-temporal multi-physics systems is computationally expensive and consists of a large number of inputs and outputs. Surrogate models are often constructed to replace the physics-based model to achieve computational efficiency in analyses such as uncertainty quantification and optimization that require many function calls. In order to address the challenge introduced by the high dimensionality of spatio-temporal output, a dimension reduction method is first employed to map the high-dimensional output to a low-dimensional latent space. This is followed by the construction of the surrogate model in the low-dimensional space. The prediction error in the original space, which includes both the reconstruction error and surrogate model error, is evaluated using different error metrics. Based on the prediction accuracy of the surrogate model, new training points are identified for adaptive improvement of the surrogate model. We present a novel adaptive sampling technique that combines exploration and exploitation to improve the surrogate model accuracy with the fewest possible runs of the expensive physics-based model. Thermo-mechanical analysis of a gas turbine engine blade is used to analyze the effectiveness of the proposed method.
Tags
Links
- Source: https://arxiv.org/abs/2608.17250v1
- Canonical: https://arxiv.org/abs/2608.17250v1
Trouble viewing inline? Open PDF directly →
Full Text
91,574 characters extracted from source content.
Expand or collapse full text
Adaptive Surrogate Modeling for High-Dimensional Spatio-Temporal Output Thanks: Citation: Kapusuzoglu, B., Matsumoto, S., Miyagi, Y., Watanabe, D., & Mahadevan, S. Adaptive surrogate modeling for high-dimensional spatio-temporal output. Structural and Multidisciplinary Optimization 65, 290 (2022). DOI: 10.1007/s00158-022-03402-x Berkcan Kapusuzoglu Thanks: Corresponding author: berkcan.kapusuzoglu@vanderbilt.edu Affiliation: Department of Civil and Environmental Engineering, Vanderbilt University, Nashville, TN 37235, USA Shunsaku Matsumoto Affiliation: Strength Research Department, Research and Innovation Center, Mitsubishi Heavy Industries, Ltd.,Nagasaki, 851-0392, Japan Yoshitomo Miyagi Affiliation: Strength Research Department, Research and Innovation Center, Mitsubishi Heavy Industries, Ltd.,Takasago, 676-8686, Japan Daigo Watanabe Affiliation: Strength Research Department, Research and Innovation Center, Mitsubishi Heavy Industries, Ltd.,Nagasaki, 851-0392, Japan Sankaran Mahadevan Affiliation: Department of Civil and Environmental Engineering, Vanderbilt University, Nashville, TN 37235, USA Abstract This paper develops an adaptive surrogate modeling method for problems with very high-dimensional spatio-temporal outputs. The analysis of spatio-temporal multi-physics systems is computationally expensive and consists of a large number of inputs and outputs. Surrogate models are often constructed to replace the physics-based model to achieve computational efficiency in analyses such as uncertainty quantification and optimization that require many function calls. In order to address the challenge introduced by the high dimensionality of spatio-temporal output, a dimension reduction method is first employed to map the high-dimensional output to a low-dimensional latent space. This is followed by the construction of the surrogate model in the low-dimensional space. The prediction error in the original space, which includes both the reconstruction error and surrogate model error, is evaluated using different error metrics. Based on the prediction accuracy of the surrogate model, new training points are identified for adaptive improvement of the surrogate model. We present a novel adaptive sampling technique that combines exploration and exploitation to improve the surrogate model accuracy with the fewest possible runs of the expensive physics-based model. Thermo-mechanical analysis of a gas turbine engine blade is used to analyze the effectiveness of the proposed method. Keywords Phsyics-informed Neural Network (PINN), Design of experiments (DoE), Adaptive sampling, Machine learning, Cross-validation, Surrogate modeling, Active learning 1 Introduction Engineering analyses such as design optimization, reliability assessment, and system health diagnosis and prognosis often require techniques such as uncertainty quantification, model calibration, and optimization, which require multiple runs of the physics model. Computational models are often used to analyze the response of an engineering system for a variety of input realizations, since conducting experiments to directly measure the true response for many input realizations is often not affordable. For large mechanical systems, expensive finite element-based or computational fluid dynamics-based physics models are commonly employed. When it is not affordable to run the expensive physics-based model many times, an inexpensive surrogate model becomes necessary to carry out repetitive analyses such as design optimization, uncertainty quantification (UQ), probabilistic diagnosis and prognosis, and risk analysis to support the decision making [33, 32]. The construction of an accurate surrogate model requires adequate amount of training data that can be generated by evaluating the physics-based model at multiple settings in the input space. The challenge regarding computational effort becomes even greater in the case of nonlinear behavior under extreme environments, as in gas turbine engine components, where the structural response is governed by a coupled multidisciplinary system of equations (fluid mechanics, heat transfer, and structural mechanics) with high-dimensional system outputs that vary over space and time. The quantity of interest (QoI) is often a multivariate output that is a field quantity (exhibits spatial dependence) and/or a stochastic process (exhibits temporal dependence). Surrogate models are often needed in such cases as mentioned above, but the surrogate model prediction quality is highly dependent on the size and distribution of the training data, which is obtained by running the original physics-based model at different input settings. For nonlinear systems with high-dimensional output, such repeated evaluations pose significant challenge w.r.t. computational resources and time, thus motivating the minimization of the number of training runs of the expensive physics-based model. Several types of surrogate models are used in the literature, e.g., polynomial chaos expansion (PCE) [53], Gaussian process (GP) regression [46], support vector regression (SVR) [6], deep neural networks (DNNs), etc. The first two approaches PCE and GP are computationally demanding in problems with high-dimensional outputs. Co-Kriging [41] has been used for multivariate outputs, but it can be computationally demanding for high-dimensional field outputs [15]. As for SVR and DNN and other machine learning models, the accuracy is dependent on the quality and quantity of the training data, which is often limited if the physics model runs are expensive [31]. Several studies on surrogate modeling techniques have addressed high-dimensional output by mapping the output to the space of principal directions, where the top few components capture most of the variance in the output, using methods such as principal components analysis (PCA) [52, 22], or singular value decomposition (SVD)) [42, 16, 24]. PCA and SVD are challenging for very high-dimensional problems since they require large storage and memory for the large covariance matrix. On the other hand, randomized SVD (rSVD) [18, 17, 10, 49] offers an efficient way to approximate the dominant singular components for high-dimensional outputs, thus allowing for a scalable architecture for modern “big data” applications. Figure 1: Schematic describing the proposed procedure for surrogate model construction with high-dimensional field data using rSVD The quality of any surrogate model is critically influenced by the set of experiments or computer simulations of the original model used for training the surrogate model. Further, if the surrogate model is constructed in a low-dimensional space, in addition to the surrogate model error in the low-dimensional space, there is also reconstruction error in translating the surrogate model prediction to the high-dimensional original space. Thus, errors need to be quantified in a systematic and useful manner to assess the quality of the surrogate model. However, often the available resources prevent obtaining a large amount of training data when the original model is computationally expensive. Therefore, it is important to construct an accurate surrogate model with the fewest possible experiments or computer simulations. Design of experiments (DoE) [48], which refers to the selection of the input values at which to conduct these experiments, becomes crucial for sufficiently accurate model construction with minimum computational or experimental cost. One of the common DoE approaches is the space-filling approach that uses a distance metric to spread the samples evenly in the design space without considering any information about the physics of the problem. In general, DoE methods can be classified as being non-adaptive or adaptive. The adaptive DoE techniques consider the physics of the system (as indicated by the previously sampled training points) while sampling new training points. Examples of non-adaptive DoE methods are the full factorial [2], Latin hypercube [39], orthogonal array [43], minimax and maximin-distance designs [27], where the focus is on coverage of the input space and not the physics of the input-output relationship. As a result, the non-adaptive DoE methods may under/oversample and waste computational resources. Whereas, in sequential adaptive sampling, the information obtained from the previous samples is used to populate new samples in regions of high interest, e.g., where the underlying function is highly nonlinear or exhibits abrupt changes and the surrogate model accuracy is poor. Moreover, with adaptive sampling it is possible to stop the computationally expensive sampling process as soon as the surrogate model accuracy reaches an acceptable level. Two components are essential in identifying additional sampling points to improve the surrogate model: (1) quantification of surrogate model error, and (2) decision criterion or learning function to identify the additional sampling points. Several approaches have been proposed in the literature to address both components. For example, cross-validation (CV) methods use estimates such as Leave-One-Out (LOO) CV error [25, 36] and new training points are located in the region of the maximum value of the CV error. Adaptive approaches based on cross-validation variance (CVV) select new samples based on the predicted maximum value of the CVV among candidate future samples [35, 47]. Hombal and Mahadevan [21] proposed adaptively selecting training points that focus on minimizing the bias in the prediction. The LOO error-based Accumulative Error (ACE) approach [36] uses a weighted combination of LOO errors. Another method is based on mean squared error (MSE) [29, 25] which chooses the next training point based on the maximum value of the estimated mean squared error in the response predicted by a GP surrogate model. The Cross-Validation Voronoi (CVVor) is based on the combination of a cross-validation exploitation with a distance-based exploration, in which a Voronoi tessellation is employed to divide the whole input space into a set of Voronoi cells [54]. The maximin scaled distance (MSD) method [25] uses a modification of maximin distance-based sampling that assigns weights to the important variables by using the available information. The Expected Improvement (EI) method uses geometry-based exploration and exploitation based on the variance of the prediction. The main goal of EI is to predict the global minimum value of the response accurately [28]. Another geometry-based exploitation method is Local Linear Approximation (LOLA)-Voronoi, which is a discontinuous adaptive sampling approach based on an exploitation feature with gradient estimation and exploration given by the volume of Voronoi tessellation cells [7]. Adaptive methods using query-by-committee-based exploration such as Mixed Adaptive Sampling Algorithm (MASA) [11] are also studied in the literature, where a new sample point is found by combining a local exploitation contribution based on Query by Committee (QBC) fluctuation and a global exploration based on distance. Adaptive sequential sampling design (also referred to as active learning [5, 55]) has also been investigated for reliability analysis, i.e., for the estimation of probability of failure; the focus of such studies has been accurate surrogate modeling of only the limit state, often formulated as g()=0g()=0, which is the boundary between regions of success and failure, instead of modeling the function g()=0g()=0 over the entire domain. The efficient global reliability analysis (EGRA) [1] method extends the expected improvement (EI) [28] mentioned above to define a new learning function Expected Feasibility Function (EFF), in the context of Gaussian process (GP) surrogate modeling of the limit state g()=0g()=0. Other active learning methods for reliability analysis are AK-MCS [12] and AK-MCS + U [44] that adaptively select new training points for the GP model improvement. Based on the AK-MCS + U method, Hu and Mahadevan [23] proposed an enhanced surrogate model-based reliability analysis method based on global sensitivity analysis to further improve the computational efficiency. All of the above-mentioned methods focus on approximating only the limit state g()=0g()=0, whereas our concern in this paper is with a general prediction model, i.e., predicting the output accurately over the entire input space. Overall, the limitations of the above-mentioned adaptive learning methods can be summarized as follows: (a) looking at a single limit state instead of the entire input space, (b) considering only single scalar output, and (c) not easily applicable to time-dependent multi-physics dynamic systems with high-dimensional spatio-temporal outputs. In order to address these limitations an adaptive surrogate modeling strategy is developed in this paper to predict the high-dimensional spatial and temporal output quantities of interest (QoIs). The proposed approach as shown in Fig. 1 addresses two challenges: high-dimensionality of the output and adaptive selection of training runs for the surrogate model. For dimension reduction, the method of randomized singular value decomposition (rSVD) [18] is used to identify the important features in the output space to give a lower dimensional representation of the original QoIs. These important features are then used to construct the surrogate model in the low dimensional space. Subsequently, the trained model is used to predict the QoIs in the original space. Error analysis is performed both in the lower dimensional latent space and the high-dimensional original space. This helps to quantify the contributions of surrogate model error and reconstruction error separately. This lays the foundation for adaptive improvement of the surrogate model. A novel approach that combines the ideas of exploration and exploitation for adaptive training point selection was developed in a manner that is applicable to time-dependent multi-physics dynamic problems with high-dimensional spatio-temporal outputs. The important contributions of this paper can be summarized as (1) Dimension reduction for very high-dimensional spatio-temporal outputs; (2) Surrogate model error quantification in two spaces; and (3) New sequential sampling technique for adaptive surrogate model improvement for multivariate spatio-temporal outputs. The rest of this paper is organized as follows. A brief introduction to dimension reduction techniques is given in Section 2. The proposed method for adaptive surrogate model construction with high-dimensional spatio-temporal output is presented in Section 3, which includes two aspects: error quantification and adaptive sampling strategy. Section 4 compares the proposed adaptive sampling technique with existing adaptive sampling algorithms, using several benchmark problems. Section 5 demonstrates the application of the proposed method for gas turbine blade thermo-mechanical analysis (with spatio-temporal, multivariate output). Section 6 provides concluding remarks and identifies future research needs. 2 Background: Output Dimension Reduction The model outputs considered in this paper vary over space and time, and are very high-dimensional. When there are many high-dimensional output QoIs at numerous spatial and temporal locations, it is not trivial to build surrogate models for every output at every location and time instant. Instead of building surrogate models directly in the original space, first the dimensionality of the high-dimensional response field can be reduced by mapping it to an uncorrelated latent space, and then a surrogate model can be constructed only for the important features in this latent space. This section discusses two techniques for dimension reduction: singular value decomposition (SVD) and randomized SVD (rSVD). Singular value decomposition (SVD) is a generalized eigen-decomposition technique for describing a large amount of high-dimensional data by mapping to a low-dimensional latent space [3]. SVD is applicable to non-square matrices and can be used to handle the spatial correlation of the response. Performing a basic SVD on a very large matrix is not only computationally expensive but also memory-intensive due to the need to store and invert very large matrices. Therefore, the randomized SVD (rSVD) algorithm, which requires less memory and avoids the high computational cost while not sacrificing accuracy, is pursued in this study to obtain a low-rank approximation of the large response matrix [18]. The randomized SVD (rSVD) method [18] maps the original data matrix to a small random subspace. The main computational operations in the rSVD-based approach used here are matrix-matrix multiplications, and QR decomposition and SVD on small matrices. More specifically, the rSVD involves five main components: (1) constructing a random matrix; (2) projection of the original data matrix to the random matrix; (3) QR decomposition of the resulting matrix; (4) multiplying the original data matrix with the resulting orthogonal basis; and (5) basic SVD on the resulting matrix. In order to obtain a low-rank, say rank r(r<<N)r\ (r<<N), approximation of a matrix ∈M×NX ^M× N, ≈rrr′X _r _rV _r (1) the target dimension, i.e., the first r pairs of singular value and singular vectors, can be obtained with the following steps. First, a random projection matrix ∈ℝN×rP ^N× r is used to get an orthonormal basis for the data matrix X, =,Z=XP, (2) where ∈ℝM×rZ ^M× r is a much smaller matrix (r<<Nr<<N) than X and approximate the column space of X with high probability due to the randomness of P. Then, a QR decomposition, which is normally employed to compute the SVD, of Z can be obtained: =QR,Z=QR^QR, (3) where ∈ℝM×rQ ^M× r and QR∈ℝr×rR^QR ^r× r. Then, the data matrix X is projected onto the subspace Q and SVD is performed on the projection ∈ℝr×NY ^r× N to obtain r _r and rTV_r^T: =T, =Q^TX, =YrrT. =U_Y _rV_r^T. (4) The matrices r _r and rTV_r^T in Eq. (2) are the same for X since Q approximates the column space of X [18]. Thus, the r left singular vectors Ur∈ℝM×rU_r ^M× r of X is r=Y.U_r=QU_Y. (5) A lower dimensional representation (dimension r) in place of the original data (dimension N) can be taken as ^LD=rr, X^LD=U_r _r, (6) Randomized SVD (rSVD) is also a linear mapping like the basic SVD. Since Z is of a much smaller size compared to X and the basis Q is low rank, the QR decomposition in Step (3) and the basic SVD in Step (5) are not time consuming. This makes rSVD a powerful tool when dealing with a very large data matrix. 3 Proposed Methodology The proposed methodology aims to efficiently construct a surrogate model for high-dimensional spatio-temporal response prediction. In this section, we first discuss dimension reduction to obtain a lower dimensional representation, features, where we adapt the rSVD method for high-dimensional spatio-temporal output using a two-step mapping strategy. Next, we investigate multiple options for surrogate modeling and select the best surrogate model for a multivariate output. This is achieved by performing a cross-validation analysis for each multivariate surrogate model constructed in a low-dimensional space for the important features. Following this, we develop the adaptive sampling strategy for identifying additional training points to improve the surrogate model, by combining exploitation and exploration. Thus the proposed overall methodology consists of four steps: 1. Dimension reduction (two-step mapping) 2. Surrogate model construction 3. Surrogate model error quantification (through cross-validation) 4. Adaptive training point selection The following subsections describe these steps in detail. 3.1 Dimension Reduction Given nsn_s data points (i.e., total number of nodes) over the spatial domain for ntn_t time domain realizations and for ndn_d realizations of design domain (i.e., training points), a data matrix for the k-th output can be written as =[(ti,θ1),(ti,θ2),…,(ti,θnd)]T _k=[D_k(t_i, _1),D_k(t_i, _2),…,D_k(t_i, _n_d)]^T =[Dk(ti,s1,θ1)Dk(ti,s1,θ2)…Dk(ti,s1,θnd)Dk(ti,s2,θ1)Dk(ti,s2,θ2)…Dk(ti,s2,θnd)⋱Dk(ti,sns,θ1)Dk(ti,sns,θ2)…Dk(ti,sns,θnd)]T = bmatrixD_k(t_i,s_1, _1)&D_k(t_i,s_1, _2)&…&D_k(t_i,s_1, _n_d)\\ D_k(t_i,s_2, _1)&D_k(t_i,s_2, _2)&…&D_k(t_i,s_2, _n_d)\\ & & & \\ D_k(t_i,s_n_s, _1)&D_k(t_i,s_n_s, _2)&…&D_k(t_i,s_n_s, _n_d) bmatrix^T (7) where (ti,θj)=[Dk(ti,s1,θj),Dk(ti,s2,θj),…,Dk(ti,sns,θj)]D_k(t_i, _j)=[D_k(t_i,s_1, _j),D_k(t_i,s_2, _j), …,D_k(t_i,s_n_s, _j)] is the i-th temporal location of the j-th realization for the k-th field response, tit_i represents the i-th temporal location and θ stands for different realizations in the design domain. The response S of the original physics-based model can be collected at the training points as follows: =[,,…,]TS=[D_1,D_2,…,D_n_k]^T (8) where nkn_k is the total number of physics model outputs. This large amount of high-dimensional response can be mapped to a low-rank approximation by using rSVD as =TS=UMV^T, where U is a (nt×nd)×ns(n_t× n_d)× n_s matrix, V is a ns×nsn_s× n_s orthogonal matrix and M is a ns×nsn_s× n_s rectangular diagonal matrix with non-negative real numbers λ=[λ1,λ2,…,λrλ=[ _1, _2,…, _r], r=min((nt×nd),ns)r=min((n_t× n_d),n_s) in the diagonal. The diagonal elements λ of M are called singular values and are arranged in descending order. The number of important features q(q≤r)q(q≤ r) to be used to represent the abstract features of the original data matrix S are determined based on the magnitudes of λ. By using the first q largest singular values, S can be reconstructed as ~ S: (ti,θj)T (t_i, _j)^T ≈~(ti,θj)=∑p=1qξp(ti,θj)p, ≈ S(t_i, _j)= _p=1^q _p(t_i, _j)V_p, (9) ∀i=1,2,…,nt;j=1,2,…,ns, ∀ i=1,2,…,n_t;\ j=1,2,…,n_s, where (ti,θj)TS(t_i, _j)^T is the (i×j)(i× j)-th row of S, ξp(ti,θj) _p(t_i, _j) is the element of ξ=ξ=UM at (i×j)(i× j)-th row and p-th column, and pV_p is the p-th important feature vector used to approximate S. The rSVD method described in Section 2 is extended in this work for dimension reduction. The following additional steps are proposed to adapt this method for a high-dimensional spatio-temporal output, in a manner that facilitates surrogate model construction. A two-step dimension reduction approach is employed to identify the important features (principal components) in the low-dimensional output space to give a lower dimensional representation of S. The dimension of the output data matrix ∈ℝ(nt×nd)×nsS ^(n_t× n_d)× n_s has a large number of features and it has more features (columns) than observations (rows). Thus, in the first step of the two-step approach the original features of the output data matrix can be reduced using rSVD to a smaller subset of features that are most relevant to the prediction problem. The result is matrix (1)∈ℝ(nt×nd)×n1I^(1) ^(n_t× n_d)× n_1 with a lower rank n1n_1 that is said to approximate the original spatio-temporal output S. Let us denote the truncated singular value matrix, eigenvector matrix, and orthogonal matrix obtained using n1n_1 important features as (1)∈ℝ(nt×nd)×n1U^(1) ^(n_t× n_d)× n_1, (1)∈ℝn1×n1 ^(1) ^n_1× n_1, (1)∈ℝn1×n1V^(1) ^n_1× n_1. The compressed form of (1)I^(1) is (1)≈(1)(1)((1))TI^(1) ^(1) ^(1)(V^(1))^T. The resulting matrix (1)I^(1) still has time-dependency. The rows (1)I^(1) include the temporal information and the realizations of design domain. The columns of (1)I^(1) represent the time series features. It was not feasible to perform a single rSVD and reduce dimensions for all quantities of interest both in time and space at the same time. We preferred to remove this temporal correlation in the second step of the two-step approach for two reasons. First, surrogate model construction was simpler for time-independent features. Second, the accuracy of the surrogate model increased as we removed the temporal variation with the proposed two-step dimension reduction approach. When rSVD is performed twice on the transformed version of matrix (1)I^(1), i.e., (1)′∈ℝnd×(n1×nt)I^(1) ^n_d×(n_1× n_t), where the rows represent the realizations of design domain, and the columns represent the time series features for each design, we obtain (2)I^(2). The columns of this low-rank matrix (2)I^(2) would represent the important features (i.e., the number of features n2n_2 that are not time series to build the surrogate model) and the rows represent the FEM simulations performed with different realizations of the design domain. Thus, we obtain (2)≈(2)(2)((2))TI^(2) ^(2) ^(2)(V^(2))^T, where (2)U^(2) ∈ ℝnd×n2R^n_d× n_2, (2) ^(2) ∈ ℝn2×n2R^n_2× n_2, (2)V^(2) ∈ ℝn2×(n1×nt)R^n_2×(n_1× n_t). The accuracy of the above two-step dimension reduction needs to be evaluated. To do this, the low-rank approximated matrix (2)I^(2) in the second feature space is mapped to the original space (i.e., first it is mapped to the first feature space and then mapped to the original space). The reconstruction accuracy is quantified using root mean square error (RMSE). The RMSE value decreases with the number of important features used to approximate the original matrix S. The number of important features is chosen based on the reconstruction accuracy and the percentage of variance explained by the top few singular vectors and singular value pairs. The proposed two-step dimension reduction approach is used in the numerical example in Section 5. Next, the surrogate model construction described in Section 3.2 is performed in the second low-dimensional space using the important features (2)=(2)(2) ξ^(2)=U^(2) ^(2). 3.2 Surrogate Model Construction Any of the available surrogate modeling techniques mentioned earlier (e.g., PCE, GP, SVR, DNN etc.) can be used to construct the surrogate model. In addition to these techniques, other prominent machine learning algorithms based on boosting (which is one of ensemble learning algorithms) are also available for improving the performance of a simple machine learning model. This paper explores several of them such as Light Gradient Boosting Machine (lightGBM) [34], Extreme Gradient Boosting (XGBoost) [4], Categorical Boosting (CatBoost) [9], and random forest (RF) [19]. Gradient Boosting Decision Tree (GBDT) is a relatively new decision tree-based ensemble learner. LightGBM, XGBoost, and CatBoost are different variations of gradient boosting methods. The main difference between these techniques is how they build the decision tree. A tree-based regression technique known as extremely randomized tree regressor or simply Extra-Trees regressor [13] is found to give the best performance in the numerical example (Section 5), based on comparing the k-fold cross validation error results for various techniques. Tree-based ensemble methods randomly construct more than one decision tree to achieve an increase in the generalization performance. For example, random forest (RF) regression [19] trains the trees with bootstrap samples for each candidate split based on randomly selected subset of the features. The main idea behind Extra-Trees is to randomly create a number of different trees with randomly chosen features. The randomization helps to achieve a greater reduction in the variance of the model prediction, compared with other ensemble methods like RF [19]. The main differences between the Extra-Trees algorithm and RF regression are: (a) Extra-Trees randomly splits nodes using a random subset of the features selected at every node, rather than the best split used in RF; and (b) RF applies the bagging procedure to iteratively generate sub-training sets (bootstrap samples) with replacements, while Extra-Trees uses all the training samples to construct each tree with varying numbers of parameters. In this work, the surrogate model is not only built for the design stage, but also for future health management of the component. Depending on how the component is used different quantities could be critical at different times. Thus, the surrogate model is envisioned to handle such changing requirements. 3.3 Surrogate Model Error Quantification The surrogate model prediction error is used to inform about areas with the highest error. The Leave-One-Out Cross-Validation (LOOCV) method is a special case of the k-fold cross-validation (CV), with k=Nk=N. For each observation i∈[1,N]i∈[1,N], a separate surrogate model ℳ−iM_-i is trained on N−1N-1 observations consisting of the reduced set −i=\(i,i)D_-i=D (x^i,y^i). The surrogate model accuracy is then evaluated on the test point ix^i. The spatio-temporal physical quantities are predicted for the test data using the surrogate model trained in the latent space, where high-dimensional spatio-temporal outputs are projected to. The surrogate model outputs in the feature space are then mapped to the original space (reconstruction) to evaluate the surrogate model accuracy in the original space. Different error metrics are suitable in different situations. For example, an engineer may be more interested in the relative error of creep damage response in the regions where strong creep behavior is expected since a small amount of change can lead to detrimental effects. Whereas the engineer is less interested in the relative error of the predicted response in regions where weak creep behavior (i.e., physical quantities are close to zero) is dominant. Thus, for this purpose, two error metrics are calculated using LOOCV, namely mean absolute error (MAE) and relative mean absolute error (rMAE). The MAE metric is defined as follows: MAEj=mean(AEj)=1ntime∑intime|Si,j−S^i,j|MAE_j=mean(AE$_j$)= 1n_time _i^n_time S_i,j- S_i,j (10) where the subscript j refers to the jjth FEM node (e.g., MAEj is the mean absolute error of jjth node, where the mean is taken across time for each spatial location). Using MAE will allow the engineers to choose a different threshold that has a physical meaning for each quantity of interest. The main challenge of the MAE metric is the identification of thresholds. Moreover, the MAE metric is scale-dependent. In order to address such limitations, the relative mean absolute error (rMAE) metric that is scale independent and less sensitive to outliers is investigated. The rMAE metric is defined as follows: rMAE=MAEmean|−¯|rMAE= MAEmean - S (11) where ¯ S is the mean of the time series output at each spatial location. An advantage of this method is its interpretability. For example, rMAE measures the possible improvement of the proposed model relative to the benchmark model (i.e., the denominator of Eq. (11)). When rMAE is less than 1, the proposed model is better than the benchmark model, and when the value of this error metric is greater than 1, then the proposed model is worse than the benchmark model. The only circumstance under which rMAE would be infinite or undefined is when all historical observations are equal (i.e., static problem). For each node in the FEM model, the error metrics are evaluated based on the LOOCV predictions. Note that the surrogate model predicts features in the low-dimensional latent space, which are then mapped to the original space to evaluate the error metrics for the spatio-temporal physical outputs. Then, PNMAE, the percentage of nodes having MAE values greater than some threshold value, and PNrMAE, the percentage of nodes having rMAE values greater than one, are calculated for the current training points using LOOCV. However, there is no way to assess these metrics for points not in the current design, i.e., unobserved points (or candidate new points, cx^c). Therefore, additional surrogate models are built to learn the relationship between the training points and the corresponding LOOCV prediction error metrics in the original space. In this study, Gaussian process (GP) surrogate models are used (allowing the use of Expected improvement (EI) [28]) to predict ePNMAEe_PNMAE and ePNrMAEe_PNrMAE for unobserved points. These GP models are developed using the initial training points. As new points are added to the design D, the training data used to train these GP models are updated. 3.4 Proposed Adaptive Sample Selection Method First, the approach generates an initial design. The initial design is chosen using a space-filling criterion with a given initial number of points. In this work, we use maximin Latin hypercube sampling (LHS) [40, 26], though any other space-filling metric can be used, to generate initial training points that cover the input space uniformly. Based on this initial design, the LOOCV error eLOOCVe_ LOOCV is calculated for each point in the design space by building a surrogate model each time leaving one point out. A secondary GP surrogate model is constructed for eLOOCVe_ LOOCV to predict eLOOCVe_ LOOCV for candidate points for additional samples. A straightforward approach is to simply use the LOOCV error for adaptive sampling, i.e., exploit regions with high prediction error to select new training points [51]. However, a purely exploitation-based approach can result in clustered samples in the input space. In order to avoid this problem, an exploration strategy based on a space-filling criterion can be added to spread the new samples over the entire input domain while exploiting the regions of interest. Thus two strategies, namely exploitation and exploration, are combined here for surrogate model improvement through adaptive selection of additional training points. Exploration aims to discover regions of the input space not covered in previous training points in order to obtain knowledge of the response over the entire design space; thus exploration does not use the simulation outputs for the previous training points. Whereas exploitation uses the output information gained from previous training points to identify high-interest subregions to generate new samples in the vicinity of these regions. These subregions can be associated with large prediction error, significant non-linear behavior, discontinuity, etc. New samples can be added in the region of interest depending on the aim of the analysis to build a surrogate model that predicts the response accurately in the exploited high-interest regions. Instead of considering these strategies individually, it can be beneficial to consider them together in hybrid adaptive learning to leverage both of their strengths, such that exploration adds points from previously unexplored regions (geometry-based) and exploitation adds points in regions of interest pertaining to surrogate model accuracy (physics-based). 3.4.1 Space-Filling Criterion A space-filling criterion is used to avoid the clustering of samples that could occur in a purely exploitation-based approach and ensure uniform coverage of the design space. The space-filling criterion used in this work is based on Euclidean distance, specifically the maximin distance, wherein the minimum non-zero distance dmind_min of a candidate point from all other points in the current set of training points is computed. A candidate point (cx^c) with the maximum of these minimum distances is selected as the new training point from ∗x^*. Thus, the exploration strategy is given by: c ^c =argmax∗∈(dmin) = _x^* (d_min ) dmin d_min =mini∈∗∈‖i−∗‖ = _ subarraycx^i \\ x^* subarray||x^i-x^*|| (12) 3.4.2 Expected Improvement In order to overcome the overexploitation problem, which could cause new samples to be clustered in regions with large mean LOOCV prediction error, a well-known approach is to use the Expected improvement (EI) [28] function, which helps to trade-off between exploitation and exploration. It can be expressed as EI()=(∗−μ())Φ(u)+σ()ϕ(u),if σ()>0.0,if σ()=0.EI(x)= cases (y^*-μ(x) ) (u)+σ(x)φ(u),&if $σ( x)>0$.\\ 0,&if $σ( x)=0$. cases (13) where u=(∗−μ())/σ()u= (y^*-μ(x) )/σ(x), ∗y^* being the current best solution chosen from among the true function values at the training points and ϕ()˙CLOSEφ( ) and Φ()˙CLOSE ( ) represent the PDF and CDF of the standard normal distribution, respectively. EI is a non-negative, parameter-free function and is zero at the training points. 3.4.3 Proposed Learning Function Combining the EI function defined in Section 3.4.2 and the space-filling criterion described in Section 3.4.1, a point with the largest expected improvement EIPNMAEEI_ PNMAE (i.e., expected improvement based on PNMAE) and EIPNrMAEEI_ PNrMAE (i.e., expected improvement based on PNrMAE) and the maximum dmind_min is selected as the new sample point. Thus, a learning function that facilitates the active learning process is proposed here as: ℒ()=(β1EIPNMAE+β2EIPNrMAE)α×(dmin())γL(x)=( _1EI_ PNMAE+ _2EI_ PNrMAE)^α×(d_min(x))^γ (14) where β1,β2∈[0,1] _1, _2∈[0,1] control the relative contributions of EIPNMAEEI_ PNMAE and EIPNrMAEEI_ PNrMAE respectively, and α,γ∈[0,1]α,γ∈[0,1] are used to control the trade-off between exploitation and exploration respectively. (Other mathematical formulations of the learning function using these metrics are also possible). In this formulation of the learning function, EI provides both local exploitation and local exploration based on the mean and variance of the prediction error, whereas the space-filling criterion only provides geometry-based global exploration. The surrogate model variance quantifies the closeness between the new training points and all existing training points. Reduction of such variance may only be a local exploration depending on the error metrics used. The error metrics used in the engineering application are dependent on the actual values of the QoIs as discussed in Section 5.4. Thus, here we used the space-filling criterion that focuses on the sample distances in the input space and guarantees global exploration. Note that ℒL provides a trade-off between the space-filling criterion, which avoids clustering, and the largest estimation of prediction error. The two extreme cases, α=0α=0 and α=1α=1, respectively denote that exploitation is not considered and exploitation is considered. Exploration of the input space is controlled by γ. When γ=0γ=0, the minimum Euclidean-distance is not considered (i.e., no exploration), whereas the larger the γ value is, the higher the weight of the exploration term. The parameters of the learning function can be estimated by using an algorithm similar to the Expectation-Maximization (EM) algorithm [8]. We have 2 sets of parameters. First, we fix the first set β1,β2 _1, _2 and compute the second set α,γα,γ in one step. Then, we compute the first set based on the computed values of the second set in the previous step. These two steps are repeated until the change in the parameter estimates between consecutive iterations does not exceed a desired threshold. The new sample point cx^c that results in the maximum learning function value ℒmaxL_max is identified as follows: ℒmax _max =max∗∈[ℒ(∗)] = _x^* [L(x^*) ] c ^c =argℒmax. = _max. (15) Note that the evaluations of candidate points using the proposed method are computationally inexpensive since these evaluations are based on the current surrogate model, not the original expensive physics model. 3.4.4 Stopping Criteria The computational resource limit (i.e., number of training points, N), and the LOOCV error can be used as the stopping criteria. The latter criterion is used to evaluate the surrogate model performance on the test set. The surrogate model accuracy is evaluated as new samples are added to the design based on the proposed adaptive sampling approach. If the LOOCV error metrics ePNMAEe_ PNMAE and ePNrMAEe_ PNrMAE achieve an acceptable value, then stop adding new points. In this work, both the computational resource limit and surrogate model accuracy are used as the stopping criteria. The detailed sequence of steps of the proposed adaptive learning approach is outlined in Algorithm 1. Algorithm 1 Identifying new training points for adaptive improvement of the surrogate model Input: D Output: cx^c 1: Generate an initial design using maximin LHS 2: procedure Output Dimension Reduction 3: Construct the data matrix, S 4: Perform a two-level rSVD on the matrix S to obtain important features 5: end procedure 6: 7: procedure Leave-One-Out Cross-Validation (LOOCV) 8: Split the design set D into k disjunct sets iD_i, i=1,…,ki=1,…,k, with k=Nk=N (LOO) 9: for i=1toNi=1\ to\ N do 10: Train an Extremely Randomized Trees machine learning model ℳ−iM_-i with N−1N-1 features consisting of the reduced set −i=\(i,i)D_-i=D (x^i,y^i). 11: Calculate ePNMAEe_ PNMAE and ePNrMAEe_ PNrMAE for the current training points in the reduced set −iD_-i: 12: ePNMAE=e_ PNMAE= % of nodes having MAE ≥ 2.5e-4 13: MAEj=mean(AEj)=1ntime∑intime|Si,j−S^i,j|MAE_j=mean(AE$_j$)= 1n_time _i^n_time S_i,j- S_i,j 14: ePNrMAE=e_ PNrMAE= % of nodes having rMAE ≥ 1 15: rMAE=MAEmean|−¯|rMAE= MAEmean - S 16: end for 17: end procedure 18: 19: procedure Adaptive Learning 20: while Average PNrMAE ≤ 15% or N≥66N≥ 66 do 21: Construct GP models for ePNMAEe_ PNMAE and ePNrMAEe_ PNrMAE (LOOCV errors) to make predictions e^PNMAE e_ PNMAE and e^PNrMAE e_ PNrMAE for candidate points 22: Define the learning function: 23: ℒ=[(β1EIPNMAE+β2EIPNrMAE)α×(dmin(∗))γ]L= [( _1EI_ PNMAE+ _2EI_ PNrMAE)^α×(d_min(x^*))^γ ] 24: Sample a large set of random points from the design space and evaluate ℒL for these points 25: Identify the largest learning function value; ℒmax=max∗∈ℒ(∗)L_max= _x^* \ L(x^*) 26: Choose the new training point cx^c that results in ℒmaxL_max, i.e., c=argℒmaxx^c= _max. 27: end while 28: end procedure 3.5 Summary of Methodology The proposed methodology for adaptive surrogate modeling with high-dimensional spatio-temporal output as shown in Fig. 1 consists of the following steps: 1. Construct the output data matrix S based on the available simulations where each row represents a time instant corresponding to one FEM run, and the columns represent multiple spatial locations of multiple output quantities for that FEM run; 2. Perform the two-step dimension reduction approach explained in Section 3.1 to find the lower dimensional representation (2)I^(2); 3. Build a single surrogate model for all features vs. process inputs in the feature space; 4. For the test set points, project the surrogate model prediction in the feature space back to the original space to obtain the predictions in the original space, S; 5. Use cross-validation to evaluate the surrogate model accuracy in the original space; 6. If the accuracy of the surrogate model is not acceptable, then use the proposed adaptive sequential sampling based on the learning function given in Section 3.4.3 to propose additional training points for the surrogate model; 7. Repeat steps 1 to 6 until the model accuracy reaches an acceptable value or until the computational resource limit reached. 4 Evaluation of Proposed Approach using Benchmark Problems In this section, the proposed adaptive sampling technique is compared to several existing techniques described in Section 1 (i.e., EI, CVV, LOLA, MSD, MASA, CVVor). To provide a fair comparison, the performance of each sampling technique in accurately capturing seven benchmark test functions of different complexity is analyzed (see Appendix A). Since the R2R^2 values for all cases were close to each other, normalized root mean square error (NRMSE = RMSE/(ymax−ymin)RMSE/(y_max-y_min), where RMSE = ∑i=1n(yi−y^i)2/n _i=1^n(y_i- y_i)^2/n) is used to be able to distinguish them from each other. NRMSE is first computed for each fold of LOOCV, followed by the computation of the average LOOCV error, by taking the average across all the LOOCV folds. This average LOOCV error is used to assess the performance of all the adaptive sampling techniques. Adaptive sampling techniques with higher exploitation component are strongly dependent on the size of the initial sample. A few empirical formulas have been proposed in the literature for choosing the initial sample size [28, 38]. In this work, the number of samples included in the initial design (m) is chosen based on the rule m=10nm=10n [28], where n is the dimension of the input space. For each two-dimensional test function in Appendix A, an initial training set consisting of 20 points is generated using maximin LHS and used to start all the adaptive sampling techniques. The positions of the 20 initial sample points are shown by black dots in Figs. 3 and 4. In order to have a fair comparison between the proposed method and the existing adaptive sampling techniques, a Gaussian process (GP) regression model is built with the exact same properties in each fold of LOOCV-based techniques. Then the number of training points is increased up to 45 by adding a new point in each step for each adaptive sampling technique, following each technique’s procedure. The NRMSE results for different sampling techniques for 45 total samples are shown in Fig. 2 for different benchmark functions. For illustration purposes, the adaptively selected samples are highlighted by red dots (larger points indicate a sample very close to a neighboring point) and the corresponding surrogate model accuracy in terms of NRMSE are shown as contours in Figs. 3 and 4 for the Branin and Schwefel functions in Appendix A respectively, for each sampling technique. Figure 2: Evolution of the NRMSE from the initial surrogate model (20 samples) to the final surrogate model (45 samples) using different adaptive sampling techniques, for the seven benchmark functions in Appendix A Figure 3: Adaptively selected training point locations for the Branin’s function and the final surrogate model accuracy in terms of NRMSE are plotted in contours (Black dots: initial 20 samples, Red dots: additional 25 samples) Figure 4: Adaptively selected training point locations for the Schwefel function and the final surrogate model accuracy in terms of NRMSE are plotted in contours (Black dots: initial 20 samples, Red dots: additional 25 samples) The Expected Improvement Cross-Validation (EICV) and Space-Filling Expected Improvement Cross-Validation (SFEICV) techniques are special cases of the proposed method; these two methods are compared against other adaptive improvement methods in this section. In EICV, γ=0γ=0, therefore the weighted minimum Euclidean distance is not considered (i.e., no additional exploration other than the exploration already inherent in EI). Whereas in SFEICV both θ and γ equal 1, thus the weight of exploration is higher. (Note that two separate GP models are constructed for EICV and SFEICV, where the first GP model is to predict the original output QoI, and the second GP model is to predict the LOOCV error as described in Section 3.4). An issue with some techniques with the exploitation component is that they may focus on insignificant features with respect to the problem of interest. When this is the case, the technique cannot capture the proper behavior of the entire function and its performance becomes highly dependent on the initial dataset. In most complex engineering problems, the initial dataset is limited, thus it is important to select an adaptive technique such as the proposed EICV that does not focus on insignificant features of the problem and exploits the regions of interest. The gradient estimation approach used in Local Linear Approximation (LOLA)-Voronoi gets increasingly complex as the problems become high-dimensional and complex. The performance of LOLA is also not as good as EICV (see Figs. 2(a)-(g)), which is a more robust and computationally efficient technique for high-dimensional problems. It should be noted that all adaptively selected samples (i.e., red dots in Figs. 3 and 4) using the LOLA technique are in regions where the initial surrogate model error was already low. This behavior of the LOLA technique can be observed in Fig. 3(e). The computational complexity of Mixed Adaptive Sampling Algorithm (MASA) depends on the number of committee members (e.g., in the context of GP models, committee members could be autocorrelation functions). In these test problems, three GP models with different autocorrelation functions (Matérn, squared exponential kernel with a different correlation length for each coordinate, and squared exponential kernel with same correlation length in all coordinates) are considered as committee surrogates (e.g., GP models based on different autocorrelation functions). MASA shows a higher emphasis on exploration (see Figs. 3(g) and 4(g)) and is not as dependent on the initial sample size as some other methods such as EI and MSD. The Maximin scaled distance (MSD) method performs poorly because the method is not capable of exploiting as much as the other methods (as seen in Figs. 2, 3 and 4). SFEICV is almost as good as MASA and MSD at exploring the design space (see Figs. 3(b) and 4(b)), and it can also exploit the regions where there is high prediction error due to local non-linearity (note the additional samples selected on the bottom left corner in Fig. 3(b)). It should be noted that Cross-Validation Voronoi (CVVor) and EICV yield the most complete performance and best accuracy across all the investigated test problems (see Fig. 2). They achieve the smallest NRMSE values, close to zero for all test functions, with the updated surrogate model (i.e., updated with 25 samples). An important difference between these adaptive techniques is the computational effort. The computational effort required by methods based on the LOOCV error such as CVVor or EICV is relatively higher, especially as the dimensions increase. However, they yield accurate results and require less development and user-knowledge. Depending on the complexity and dimension of the problem, the SFEICV technique may be useful in the first few adaptively selected samples, but as more samples are collected, the EICV method becomes more promising. This can be seen in Figs. 2(a)(d)(e)(f) where SFEICV yields similar NRMSE values to EICV after 10 adaptively selected training points, i.e., 30 total training points, whereas, EICV results in better improvement in the surrogate model accuracy after 30 samples. Note that two separate GP models are constructed in EICV and SFEICV (which are special cases of the general method proposed in Section 3.4.3), thus they require more computational effort than the other adaptive sampling techniques that are not based on the LOOCV error. However, they do not require the tessellation of the input space into Voronoi cells, which also requires high effort, as in CVVor and LOLA. In addition, the accuracy of EICV and SFEICV is better or comparable to the other methods for different benchmark functions. 5 Engineering Application A simplified gas turbine engine blade model (see Fig. 5) is studied in this section to demonstrate the proposed surrogate modeling approach for an engineering application with high-dimensional spatio-temporal output. Six output quantities from the physics model (creep equivalent strain (CEEQ), creep damage (Dc), von Mises stress (Mises), and displacements in x, y, and z directions) are of interest, as shown in Table 1; these are available at a large number of spatial and temporal points. The input space is four-dimensional and consists of turbine blade coating thickness (TBC thickness), turbine output rate (T1T rate), and transient and constant creep rates. The ranges of the inputs are: TBC thickness [2/3*nominal,4/3*nominal], T1T rate [80*nominal,100*nominal], transient and constant creep rate [-3σ,3σ]. These model inputs are considered as uncertain. The boundary conditions are assumed to be known with certainty. The total number of nodes is 29,374 for each output quantity of interest (QoI) and the total number of time steps is 54. Thus, the problem has over 176,244 spatial dimensions (i.e., 29,374 × 6) in the output space for a single FEM simulation. Considering also the temporal dimension, the problem dimension is over 9,000,000 (i.e., 176,244 × 54) just for a single FEM simulation. Table 1: Physics model input and output parameters Input Values TBC thickness 23× 23× nominal to 43× 43× nominal T1T rate 80 × nominal to 100 × nominal Transient creep rate -3σ to 3σ Constant creep rate -3σ to 3σ Output Nomenclature Creep equivalent strain CEEQ Creep damage Dc Von Mises stress Mises X displacement X Y displacement Y Z displacement Z Figure 5: FEM model of engine rotor blade 5.1 Generation of Initial Training Data The initial surrogate model is built with 46 training points (from 46 FEM runs). Prior to this study, nine initial input settings had been generated using Taguchi L9 orthogonal design [50] for exploratory analysis. To this set of 9 points, we added 37 points, generated using maximin Latin hypercube sampling (LHS). These 46 sets of inputs and the corresponding FEM outputs are then used to build the initial surrogate model. In order to obtain the maximin LHS for the 4-dimensional input space, 100 LHS designs are generated. Based on the spread of the training points in the design space, the design with the highest spatial coverage is selected as the optimum. To quantify how spread out the points are in each design set, the sum of pair-wise Minkowski distances when p = 1 within the set is selected as a metric. For two row vectors, ∈ℝnx ^n and ∈ℝny ^n in the matrix, the Minkowski distance is defined as: d(,)=(∑i=1n|i−i|p)(1/p).d(x,y)= ( _i=1^n _i-y_i ^p )^(1/p). (16) The higher the d(,)d(x,y) is, the higher the occupancy rate (coverage) in the design space. The sum of the pair-wise Minkowski distances was calculated for each LHS design and the optimum design is selected to give the 37 additional training points for the initial surrogate model. 5.2 Dimension Reduction of the Output Space Based on the correlation analysis, CEEQ and Dc are found to be perfectly correlated with each other at all time steps [30]. Similarly, Y and Z displacements are negatively correlated, and the correlation becomes more significant at higher time steps. The features obtained in single-physics problem like structural dynamics have physical meaning, but in multi-physics problems where different output quantities are combined together, the physical meaning is not clear. In this work, dimension reduction is performed using rSVD on the entire multi-response output data in order to take this correlation between the output QoIs into account and improve the computational efficiency. Thus the obtained features represent all the QoIs. Then a single surrogate model, Extra-Trees regressor, is built to predict multivariate time series for all 6 QoIs. To be able to jointly predict the QoIs at each time step for untested configurations, the output data matrix S is constructed as follows: =[123456](ntime×nsim)×(nodes×nQoI)S=[D_1D_2D_3D_4D_5D_6]_(n_time× n_sim)×(n_nodes× n_QoI) (17) where i∈D_i∈ ℝntime×nodesR^n_time× n_nodes, i=(1,2,…,nQoI)\ i=(1,2,...,n_QoI) is the iith QoI data matrix. The matrix ∈S∈ ℝ(nsim×ntime)×(nodes×nQoI)R^(n_sim× n_time)×(n_nodes× n_QoI) is constructed such that the number of rows is equal to the total number simulations (nsimn_sim) times the total number of time steps of each simulation (ntimen_time), and the number of columns is equal to number of nodes (nodesn_nodes) times the total number of QoIs (nQoIn_QoI). For example, for the final surrogate model we have output data for 66 simulations and the output data matrix (54×66)×(29374×6)S_(54× 66)×(29374× 6) consists of approximately 629 million data points. We have used RMSE to evaluate the reconstruction accuracy and the percentage of variance explained by the top important features for further analysis. The percentage of variance explained by the top 50 and 20 features in the first and second feature spaces, respectively, is 99% and the RMSE value for reconstruction error is approximately 0.0006, which is small compared to the average magnitude of the predicted quantities (0.01). We have not used the MAE metric for calculating the reconstruction error. MAE will allow the engineers to choose a different threshold that has a physical meaning for each individual quantity of interest. However, the values of MAE thresholds for reconstruction error in individual quantities are not known to the engineer. Whereas, when we use RMSE, we can quantify one error measure for all the quantities together and assess the quality of reconstruction. 5.3 Surrogate Model Construction and Cross-Validation After the important features are identified, several types of surrogate models are investigated as discussed in Section 3.2 and based on the cross-validation results given in Table. 2, the Extra-Trees regressor is observed to give the best results with an average RMSE value of 0.10513 (15% better than the next best model). We have used 7-fold in Table 2 to partition the data evenly (k=7 is a divisor of the sample size). We tried both with 5-fold and 10-fold. The cross-validation results were almost identical. The LOOCV error is used to evaluate the surrogate model accuracy. The inputs to the surrogate model are TBC thickness, T1T rate, transient and constant creep rates and the outputs of the surrogate model are the top 20 important features. The proposed approach is able to distinguish between errors stemming from bad approximation of the surrogate model and errors stemming from insufficient quality of the project space. We have not distinguished them in our numerical example, but have included both of them together in our analysis, especially since the reconstruction error was negligible compared to the surrogate model error. Table 2: 7-fold cross-validation results in terms of RMSE Model Fold # Extra-Trees LightGBM XGBoost CatBoost DNN GP 1 0.1722 0.2029 0.1823 0.1780 0.2102 0.2205 2 0.1885 0.2146 0.2113 0.1843 0.2243 0.2315 3 0.1329 0.1829 0.1732 0.1362 0.1856 0.1795 4 0.1090 0.1257 0.1219 0.1155 0.1351 0.1583 5 0.1123 0.1286 0.1199 0.1136 0.1209 0.1348 6 0.1039 0.1273 0.1166 0.1065 0.1385 0.1421 7 0.1426 0.1835 0.1462 0.1503 0.1857 0.1812 Average 0.1373 0.1665 0.1531 0.1406 0.1714 0.1782 For illustration purposes, the predicted feature values versus the actual feature values are shown in Fig. 6 for FEM run 23. Note that the pairs of observations and predictions are close to the 45-degree line, showing good agreement. The R2R^2 and RMSE values of the predictions shown in Fig. 6 are 0.97 (quantifies how close the observed vs. predicted data are to the 45-degree line) and 0.45 (small error compared to the actual feature values), respectively. Figure 7 shows actual Feature 1 vs. the input parameters TBC thickness and T1T rate. Although the general trend shows a decrease of Feature 1 values with an increase in TBC thickness, there are still small Feature 1 values across all values of TBC thickness. Feature 1 values increase as T1T rate increases, indicating a possible quadratic relationship. From the physics of the problem, it is known that T1T values greater than 90 result in strong creep cases (i.e., large CEEQ and Dc values) and Feature 1 values corresponding to these cases are significantly larger than Feature 1 values corresponding to the weak creep cases (i.e., T1T << 90). Moreover, the correlations between the model inputs and the top 5 important features are computed (not shown in the paper due to space limitation), and the most significant correlation is found to be between T1T rate and Feature 1, with a value of 0.85. Figure 6: Predicted vs. observed values of important features for FEM run 23 Figure 7: Actual Feature 1 values vs. TBC thickness and T1T rate 5.4 Adaptive Surrogate Modeling Next, the training points corresponding to additional FEM runs (47-66) are adaptively selected based on the proposed methodology of adaptive training point selection. These runs are completed in 4 batches (see Fig. 8 and Table 3). From the surrogate model improvement point of view, it is more desirable to select one training point at a time. However, this is not possible in industrial settings due to the fact that running a single FEM model at each step is inefficient use of human and computational resources. Table 3: Optimal hyperparameter values of each batch Batch id β1 _1 β2 _2 α γ First batch: FEM 47-56 0.30 0.03 1 1 Second batch: FEM 57-60 0.56 0.23 1 1 Third batch: FEM 61-62 1 1 1 1 Fourth batch: FEM 63-66 0.53 0.04 1 1 A more efficient strategy in this application was to run between 2 and 10 FEM models at the same time. Thus, multiple new training points were proposed in each batch where the number of points for each batch is chosen based on error quantification, space filling, as well as expert opinion regarding inputs of interest. The simplest strategy would be to rank the candidate points and select the desired number of best points. However, this is not ideal since it does not consider the information overlapping of the new points in that batch, which can lead to clustered batch points [37]. (This is because the space-filling component of the learning function in Section 3.4.3 only gives the score for each single candidate point in terms of its distance from the existing points; whereas when we consider multiple points in a batch, we also have to consider the distances among them). In our case, the selection of the new training points in each batch considers how informative and diverse they are [37], from the perspectives of both geometry and physical behavior. This way not only are the new points sampled in the regions of interest but they are also far away from each other. Efficient global optimization (EGO) techniques that select the next points in batches extended the expected improvement (EI) such that EI is maximized when multiple points are added to the data set [20, 45]. However, the optimization of the multiple point learning functions is computationally challenging [14]. Here we considered a simpler and computationally cheaper approach (i.e., multiple points were selected based on their learning function values out of all candidate points) for selecting the next points in batches. For example, in the second batch, only four (57-60) samples were selected that satisfied the criteria of informativeness and diversity based on the proposed learning function values out of all candidate points. Similarly, in the third and fourth batches, only two points each were selected for the same reasons. The first batch of new training points are selected by performing the proposed methodology described in Algorithm 1 on the previous 46 FEM runs, and the second batch of new training points are obtained by doing the same analysis on the previous 56 FEM runs (i.e., 46 initial points plus 10 points from the first batch of additional training points), and so on. For most of the adaptively sampled points, i.e., FEM 47-60, and 63-66, the parameters of the learning function are estimated by using the two-step algorithm described in Section 3.4.3. The expert opinion is only used for some of the training points to guide the learning. For example, when the relative error was more important for the engineers, we have adjusted beta in an ad-hoc fashion. Most of these samples are from regions where the absolute error has high values (i.e., physical quantities take medium to high values, which is of interest to the analyst in this case). However, for the problem of interest it is also equally important to be able to accurately predict physical quantities that are close to zero. (For example, in the gas turbine engine, some input values might lead to low creep; such cases could be quickly screened and ignored in further detailed analysis). Thus, for a few of the adaptively selected additional samples such as FEM 61 and 62 the parameters of the learning function shown in Eq. (14) are chosen as unity to give equal weight to all metrics in order to facilitate choosing more diverse training points, including from regions where the relative error is higher, within each batch (e.g., FEM 64 and 65). Figure 8: Initial FEM training data and adaptively selected additional samples The surrogate model prediction accuracy is further analyzed in the original space by mapping the surrogate model predictions in the features space back to the original space. As discussed earlier, the percentage of nodes having MAE >> threshold (PNMAE) and the percentage of nodes having rMAE >> unity (PNrMAE) are defined for error analysis. The threshold value of 0.00025 is chosen based on expert opinion. The accuracy of the final surrogate model is shown in Fig. 9 using the two error metrics PNMAE and PNrMAE in a two-dimensional space, i.e., TBC thickness vs T1T rate. The error metrics have different values in different regions. For example, the absolute error metric, PNMAE, has smaller values in the bottom half region, where the physical quantities are close to zero, and has larger values in the top left region, where the physical quantities take medium to high values. Whereas, the relative error of the predicted creep response, PNrMAE, has greater values in regions where weak creep behavior is dominant. This is caused by dividing with a very small number. Since the physical quantities are close to zero in weak creep regions (i.e., the bottom half region T1T <90<90) the denominator of Eq. (11) would be close to zero. Figure 9: Final surrogate model prediction accuracy in the original space for leave-one-out FEM runs The final surrogate model predictions after 66 training points at an important node (id 5057: hot spot, based on analyst’s knowledge) are compared against the FEM results of run 36, where high creep damage is observed, as shown in Fig. 10. It is observed that the predictions are in good agreement with the FEM results at all time steps (similar trends in some cases if not actual values), which reflects the effectiveness of the proposed methodology. Figure 10: Final surrogate model prediction compared against FEM run 36 at node 5057 The LOOCV results are shown in Table 4 using three different error metrics. In order to have a fair comparison, the improvement in the surrogate model accuracy is assessed on the same data set (i.e., first 46 FEM runs). The columns of Table 4 represent the average LOOCV error values on the first 46 FEM runs. In general, the accuracy of the surrogate model is observed to improve with each batch of adaptively selected samples. The average RMSE value, which is obtained by averaging the RMSE values obtained for each LOOCV fold, decreases with each new batch of adaptively selected training points. For example, it decreases from 0.00254 to 0.00222 with the addition of 10 new adaptively selected training points (samples 47-56) based on the LOOCV on the first 46 FEM runs. Similarly, PNMAE and PNrMAE values decrease as we include adaptively selected training points during the surrogate model construction process. Table 4: Average LOOCV results of the first 46 FEM runs in terms of PNMAE, PNrMAE, and RMSE Training data Avg. PNMAE Avg. PNrMAE Avg. RMSE FEM 1-46 27.9 19.3 0.00254 FEM 1-56 26.5 17.6 0.00222 FEM 1-60 26.3 17.3 0.00219 FEM 1-62 24.9 16.1 0.00209 FEM 1-66 24.6 15.8 0.00193 The number of training points N and the surrogate model accuracy are used as stopping criteria. The sequential adaptive sampling procedure is stopped after 66 training points, N=66N=66. The adaptively improved surrogate model does not successfully reach the desired level of accuracy of average PNrMAE ≤ 15% (chosen based on expert opinion). However, the average PNrMAE value is 15.8% after 66 training points (for the first initial 46 FEM runs), which suggests that with a few more FEM runs the desired level of accuracy could be reached. In addition to the average PNrMAE value being distinctly close to the desired level of accuracy, the PNMAE and RMSE values are also acceptable for the analyst after 66 training points. The largest contributor to the overhead cost is the adaptive sampling procedure. The second largest contributor is the validation of the surrogate model for each QoI using the appropriate error metric. All calculations except the FEM simulations are performed on a single desktop computer with an Intel 8-core CPU, 3.00 GHz base frequency, and 8 GB memory. 6 Conclusion This paper developed an approach for adaptive surrogate model construction in engineering problems with high-dimensional spatio-temporal output. The important features are obtained by performing a two-step dimension reduction using rSVD to map the original high-dimensional spatio-temporal output to an uncorrelated space. Cross-validation error is used to identify the most accurate surrogate model type. Once the most accurate surrogate model type is identified for the problem of interest, the prediction error in the original space is evaluated using LOOCV with different error metrics (i.e., PNMAE and PNrMAE). The subsequent adaptive sampling technique combines exploration and exploitation for adaptive improvement of surrogate model accuracy with the fewest possible runs of the expensive physics-based original model. The proposed adaptive sampling technique is first compared with some of the existing adaptive schemes (Section 4) using several test problems. The proposed method yields very good results across all investigated test problems. The effectiveness of the proposed methodology is next demonstrated for a turbine blade example (Section 5), using a time-dependent multi-physics dynamic system model with high-dimensional spatio-temporal outputs. It is observed that the adaptively improved surrogate model reaches an acceptable level of accuracy using the proposed strategy. Future work can explore other large-scale engineering applications in order to further evaluate the performance of the proposed method, gain more insight about the method’s behavior for different problems, and further improve the methodology. The proposed adaptive sampling technique can be extended in future to include parallelization of sequential sampling for global optimization. Future work can also isolate the reconstruction error from the overall prediction error in the original space, which includes both the reconstruction error and surrogate model error. Acknowledgments This study was funded by Mitsubishi Heavy Industries in Nagasaki and Takasago, Japan. The support is gratefully acknowledged. In addition, we acknowledge our paper in the AIAA SCITECH 2022 Forum [30]; note that in that paper, we only presented preliminary results but did not discuss the methodology. Declarations The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Conflict of Interest The authors declare that they have no conflict of interest. Replication of Results Unfortunately, we cannot share the in-house code used in this work because of confidential proprietary nature as identified by the sponsor. Appendix A Benchmark Test Functions Two-dimensional benchmark functions used in this paper to evaluate the performance of the proposed method and compare its performance to existing adaptive surrogate model improvement methods are defined in Table 5. Table 5: Test problems Branin functiony=[x2−5(x12π)2+5x1π−6]+10(1−18π)cos(x1)+10x1∈[−5,10],x2∈[0,15] aligned &Branin function\\ &y= [x_2-5 ( x_12π )^2+ 5x_1π-6 ]+10 (1- 18π )cos(x_1)+10\\ &x_1∈[-5,10],\ x_2∈[0,15] aligned Goldstein & Price functiony=[1+(x1+x2+1)2(19−14x1+3x12−14x2+6x1x2+3x22)]×[30+(2x1−3x2)2(18−32x1+12x12+48x2−36x1x2+27x22)]x1,x2∈[−2,2] aligned &Goldstein \& Price function\\ &y= [1+ (x_1+x_2+1 )^2 (19-14x_1+3x^2_1-14x_2+6x_1x_2+3x^2_2 ) ]\\ & × [30+ (2x_1-3x_2 )^2 (18-32x_1+12x^2_1+48x_2-36x_1x_2+27x^2_2 ) ]\\ &x_1,\ x_2∈[-2,2] aligned Sasena’s functiony=2+0.01(x2−x12)2+(1−x1)2+2(2−x2)+7sin(0.5x1)sin(0.7x1x2)x1,x2∈[0,5] aligned &Sasena's function\\ &y=2+0.01(x_2-x_1^2)^2+(1-x_1)^2+\\ & 2(2-x_2)+7sin(0.5x_1)\ sin(0.7x_1x_2)\\ &x_1,\ x_2∈[0,5] aligned Alpine functiony=sinx1sinx2(x1x2)x1,x2∈[0,10] aligned &Alpine function\\ &y=sinx_1\ sinx_2 (x_1x_2)\\ &x_1,\ x_2∈[0,10] aligned Modified form of function in Meckesheimer et al. [meckesheimer2001computationally]y=e(x1−x2)2+e(10−x1)2−x1x2x1,x2∈[0,10] aligned &Modified form of function in Meckesheimer et al.~ [cite][\@@bibrefmeckesheimer2001computationally]\\ &y=e^(x_1-x_2)^2+e^(10-x_1)^2-x_1x_2\\ &x_1,\ x_2∈[0,10] aligned Modified form of function in Jin et al. [25]y=cos(10x12)+3.1|x1−0.7|+2x12+sin(1|x1−0.7|+0.31)+2x22x1,x2∈[0,1] aligned &Modified form of function in Jin et al.~ [cite][\@@bibrefjin2002sequential]\\ &y=cos(10x_1^2)+3.1 x_1-0.7 +2x_1^2+sin ( 1 x_1-0.7 +0.31 )+2x_2^2\\ &x_1,\ x_2∈[0,1] aligned Schwefel functiony=418.9829n−∑(xisin|xi|CLOSExi∈[−500,500],i=1,…,n,n=2 aligned &Schwefel function\\ &y=418.9829n-Σ(x_i\ sin x_i \\ &x_i∈[-500,500],\ i=1,...,n,n=2 aligned References [1] B. J. Bichon, M. S. Eldred, L. P. Swiler, S. Mahadevan, and J. M. McFarland (2008) Efficient global reliability analysis for nonlinear implicit performance functions. AIAA journal 46 (10), p. 2459–2468. External Links: Document Cited by: §1. [2] G. E. Box and J. S. Hunter (1961) The 2 k—p fractional factorial designs. Technometrics 3 (3), p. 311–351. External Links: Document Cited by: §1. [3] A. Chatterjee (2000) An introduction to the proper orthogonal decomposition. Current Science 78 (7), p. 808–817. External Links: ISSN 00113891, Document Cited by: §2. [4] T. Chen and C. Guestrin (2016) XGBoost: a scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, New York, NY, USA, p. 785–794. External Links: ISBN 9781450342322, Document Cited by: §3.2. [5] D. A. Cohn, Z. Ghahramani, and M. I. Jordan (1996) Active learning with statistical models. Journal of Artificial Intelligence Research 4, p. 129–145. External Links: Document Cited by: §1. [6] C. Cortes and V. Vapnik (1995) Support-vector networks. Machine learning 20 (3), p. 273–297. External Links: Document Cited by: §1. [7] K. Crombecq, D. Gorissen, D. Deschrijver, and T. Dhaene (2011) A novel hybrid sequential design strategy for global surrogate modeling of computer experiments. SIAM Journal on Scientific Computing 33 (4), p. 1948–1974. External Links: Document Cited by: §1. [8] A. P. Dempster, N. M. Laird, and D. B. Rubin (1977) Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society: Series B (Methodological) 39 (1), p. 1–22. External Links: Document Cited by: §3.4.3. [9] A. V. Dorogush, V. Ershov, and A. Gulin (2018) CatBoost: gradient boosting with categorical features support. External Links: 1810.11363 Cited by: §3.2. [10] P. Drineas, R. Kannan, and M. W. Mahoney (2006) Fast monte carlo algorithms for matrices i: computing a low-rank approximation to a matrix. SIAM Journal on Computing 36 (1), p. 158–183. External Links: Document Cited by: §1. [11] J. Eason and S. Cremaschi (2014) Adaptive sequential sampling for surrogate model generation with artificial neural networks. Computers & Chemical Engineering 68, p. 220–232. External Links: Document Cited by: §1. [12] B. Echard, N. Gayton, and M. Lemaire (2011) . Structural Safety 33 (2), p. 145–154. External Links: ISSN 0167-4730, Document Cited by: §1. [13] P. Geurts, D. Ernst, and L. Wehenkel (2006) Extremely randomized trees. Machine learning 63 (1), p. 3–42. External Links: Document Cited by: §3.2. [14] D. Ginsbourger, R. Le Riche, and L. Carraro (2010) Kriging is well-suited to parallelize optimization. In Computational Intelligence in Expensive Optimization Problems, Y. Tenne and C. Goh (Eds.), p. 131–162. External Links: Document Cited by: §5.4. [15] C. Gogu and J. Passieux (2013) Efficient surrogate construction by combining response surface methodology and reduced order modeling. Structural and Multidisciplinary Optimization 47 (6), p. 821–837. External Links: Document Cited by: §1. [16] Y. Guo, S. Mahadevan, S. Matsumoto, S. Taba, and D. Watanabe (2021) Surrogate modeling with high-dimensional input and output. In AIAA Scitech 2021 Forum, p. 0182. External Links: Document Cited by: §1. [17] N. Halko, P. G. Martinsson, and J. A. Tropp (2011) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review 53 (2), p. 217–288. External Links: Document Cited by: §1. [18] N. Halko, P. Martinsson, and J. A. Tropp (2011) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM review 53 (2), p. 217–288. External Links: Document Cited by: §1, §1, §2, §2, §2. [19] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman (2009) The elements of statistical learning: data mining, inference, and prediction. Springer Series in Statistics, Springer New York, New York City, USA. External Links: ISBN 9780387216065, Document Cited by: §3.2, §3.2. [20] N. Henkenjohann and J. Kunert (2007) An efficient sequential optimization approach based on the multivariate expected improvement criterion. Quality Engineering 19 (4), p. 267–280. External Links: Document Cited by: §5.4. [21] V. Hombal and S. Mahadevan (2011) Bias minimization in gaussian process surrogate modeling for uncertainty quantification. Visualization of Mechanical Processes: An International Online Journal 1 (4). External Links: Document Cited by: §1. [22] V. Hombal and S. Mahadevan (2013) Surrogate modeling of 3d crack growth. International Journal of Fatigue 47, p. 90–99. External Links: Document Cited by: §1. [23] Z. Hu and S. Mahadevan (2016) Global sensitivity analysis-enhanced surrogate (gsas) modeling for reliability analysis. Structural and Multidisciplinary Optimization 53 (3), p. 501–521. External Links: Document Cited by: §1. [24] Z. Hu and S. Mahadevan (2017) A surrogate modeling approach for reliability analysis of a multidisciplinary system with spatio-temporal output. Structural and Multidisciplinary Optimization 56 (3), p. 553–569. External Links: Document Cited by: §1. [25] R. Jin, W. Chen, and A. Sudjianto (2002) On sequential sampling for global metamodeling in engineering design. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. 36223, p. 539–548. External Links: Document Cited by: §1. [26] R. Jin, W. Chen, and A. Sudjianto (2003) An efficient algorithm for constructing optimal design of computer experiments. In International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Vol. 37009, p. 545–554. External Links: Document Cited by: §3.4. [27] M. E. Johnson, L. M. Moore, and D. Ylvisaker (1990) Minimax and maximin distance designs. Journal of Statistical Planning and Inference 26 (2), p. 131–148. External Links: Document Cited by: §1. [28] D. R. Jones, M. Schonlau, and W. J. Welch (1998) Efficient global optimization of expensive black-box functions. Journal of Global Optimization 13 (4), p. 455–492. External Links: Document Cited by: §1, §1, §3.3, §3.4.2, §4. [29] D. R. Jones (2001) A taxonomy of global optimization methods based on response surfaces. Journal of Global Optimization 21 (4), p. 345–383. External Links: Document Cited by: §1. [30] B. Kapusuzoglu, Y. Guo, S. Mahadevan, S. Matsumoto, M. Yoshitomo, S. Taba, and D. Watanabe (2022) Dimension reduction for efficient surrogate modeling in high-dimensional applications. In AIAA SCITECH 2022 Forum, p. 1440. External Links: Document Cited by: §5.2, Acknowledgments. [31] B. Kapusuzoglu and S. Mahadevan (2021) Information fusion and machine learning for sensitivity analysis using physics knowledge and experimental data. Reliability Engineering & System Safety 214, p. 107712. External Links: ISSN 0951-8320, Document Cited by: §1. [32] B. Kapusuzoglu, P. Nath, M. Sato, S. Mahadevan, and P. Witherell (2022) Multi-Objective Optimization Under Uncertainty of Part Quality in Fused Filament Fabrication. ASCE-ASME J Risk and Uncert in Engrg Sys Part B Mech Engrg 8 (1). Note: 011112 External Links: ISSN 2332-9017, Document, https://asmedigitalcollection.asme.org/risk/article-pdf/8/1/011112/6820272/risk_008_01_011112.pdf Cited by: §1. [33] B. Kapusuzoglu, M. Sato, S. Mahadevan, and P. Witherell (2020) Process optimization under uncertainty for improving the bond quality of polymer filaments in fused filament fabrication. Journal of Manufacturing Science and Engineering 143, p. 1–46. External Links: ISSN 1087-1357, Document Cited by: §1. [34] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T. Liu (2017) LightGBM: a highly efficient gradient boosting decision tree. In Advances in Neural Information Processing Systems, I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Eds.), Vol. 30, Long Beach, California, USA, p. 3149–3157. Cited by: §3.2. [35] J. P. Kleijnen and W. C. Van Beers (2004) Application-driven sequential designs for simulation experiments: kriging metamodelling. Journal of the Operational Research Society 55 (8), p. 876–883. External Links: Document Cited by: §1. [36] G. Li, V. Aute, and S. Azarm (2010) An accumulative error based adaptive design of experiments for offline metamodeling. Structural and Multidisciplinary Optimization 40 (1-6), p. 137. External Links: Document Cited by: §1. [37] H. Liu, Y. Ong, and J. Cai (2018) A survey of adaptive sampling for global metamodeling in support of simulation-based complex engineering design. Structural and Multidisciplinary Optimization 57 (1), p. 393–416. External Links: Document Cited by: §5.4. [38] H. Liu, S. Xu, Y. Ma, X. Chen, and X. Wang (2016) An adaptive bayesian sequential sampling approach for global metamodeling. Journal of Mechanical Design 138 (1), p. 011404. External Links: Document Cited by: §4. [39] M. D. McKay, R. J. Beckman, and W. J. Conover (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), p. 239–245. External Links: Document Cited by: §1. [40] M. D. Morris and T. J. Mitchell (1995) Exploratory designs for computational experiments. Journal of Statistical Planning and Inference 43 (3), p. 381–402. External Links: Document Cited by: §3.4. [41] D. E. Myers (1982) Matrix formulation of co-kriging. Journal of the International Association for Mathematical Geology 14 (3), p. 249–257. External Links: Document Cited by: §1. [42] P. Nath, Z. Hu, and S. Mahadevan (2017) Sensor placement for calibration of spatially varying model parameters. Journal of Computational Physics 343, p. 150–169. External Links: Document Cited by: §1. [43] A. B. Owen (1992) ORTHOGONAL arrays for computer experiments, integration and visualization. Statistica Sinica 2 (2), p. 439–452. External Links: ISSN 10170405, 19968507, Document Cited by: §1. [44] Z. Peijuan, W. C. Ming, Z. Zhouhong, and W. Liqi (2017) A new active learning method based on the learning function u of the ak-mcs reliability analysis method. Engineering Structures 148, p. 185–194. External Links: ISSN 0141-0296, Document Cited by: §1. [45] W. Ponweiser, T. Wagner, and M. Vincze (2008) Clustered multiple generalized expected improvement: a novel infill sampling criterion for surrogate models. In 2008 IEEE congress on evolutionary computation (IEEE world congress on computational intelligence), p. 3515–3522. External Links: Document Cited by: §5.4. [46] C. E. Rasmussen (2004) Gaussian processes in machine learning. In Advanced Lectures on Machine Learning: ML Summer Schools 2003, Canberra, Australia, February 2 - 14, 2003, Tübingen, Germany, August 4 - 16, 2003, Revised Lectures, O. Bousquet, U. von Luxburg, and G. Rätsch (Eds.), p. 63–71. External Links: ISBN 978-3-540-28650-9, Document Cited by: §1. [47] D. A. Romero, C. H. Amon, and S. Finger (2006) On adaptive sampling for single and multi-response bayesian surrogate models. In ASME 2006 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, p. 393–404. External Links: Document Cited by: §1. [48] J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn (1989) Design and analysis of computer experiments. Statistical science, p. 409–423. External Links: Document Cited by: §1. [49] T. Sarlos (2006) Improved approximation algorithms for large matrices via random projections. In 2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), Vol. , p. 143–152. External Links: Document Cited by: §1. [50] G. Taguchi and M. S. Phadke (1989) Quality engineering through design optimization. In Quality Control, Robust Design, and the Taguchi Method, p. 77–96. External Links: ISBN 978-1-4684-1472-1, Document Cited by: §5.1. [51] F. A. Viana, C. Gogu, and T. Goel (2021) Surrogate modeling: tricks that endured the test of time and some recent developments. Structural and Multidisciplinary Optimization 64 (5), p. 2881–2908. External Links: Document Cited by: §3.4. [52] S. Wold, K. Esbensen, and P. Geladi (1987) Principal component analysis. Chemometrics and Intelligent Laboratory Systems 2 (1), p. 37–52. Note: Proceedings of the Multivariate Statistical Workshop for Geologists and Geochemists External Links: ISSN 0169-7439, Document Cited by: §1. [53] D. Xiu and G. E. Karniadakis (2002) The wiener–askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing 24 (2), p. 619–644. External Links: Document Cited by: §1. [54] S. Xu, H. Liu, X. Wang, and X. Jiang (2014) A robust error-pursuing sequential sampling approach for global metamodeling based on voronoi diagram and cross validation. Journal of Mechanical Design 136 (7), p. 071009. External Links: Document Cited by: §1. [55] X. Yang, Y. Liu, Y. Gao, Y. Zhang, and Z. Gao (2015) An active learning kriging model for hybrid reliability analysis with both random and interval variables. Structural and Multidisciplinary Optimization 51 (5), p. 1003–1016. External Links: Document Cited by: §1.