Advances in Bayesian Probabilistic Modeling for Industrial Applications
Sayan Ghosh Piyush Pandita Steven Atkinson Waad Subber Yiming Zhang Natarajan Chennimalai Kumar Suryarghya Chakrabarti Liping Wang
stat.ML cs.LG


Bag of biterms modeling for short texts[This paper is an extension version of the PAKDD2016 Long Presentation paper “Mai, K., Mai, S., Nguyen, A., Van Linh, N.,& Than, K. (2016). Enabling Hierarchical Dirichlet Processes to Work Better for Short Texts at Large Scale. In Advances in Knowledge Discovery and Data Mining (pp. 431-442). Springer International Publishing.]

    Anh Phan Tuan Bach Tran Thien Nguyen Huu Linh Ngo Van Khoat Than

    Received: date / Accepted: date
==================================================================================================================================================================================================================================================================================================================================================================================





Industrial  applications frequently pose a notorious challenge for state-of-the-art methods in the contexts of optimization, designing experiments and modeling unknown physical response. This problem is aggravated by limited availability of clean data, uncertainty in available physics-based models and additional logistic and computational expense associated with experiments. In such a scenario, Bayesian methods have played an impactful role in alleviating the aforementioned obstacles by quantifying uncertainty of different types under limited resources. These methods, usually deployed as a framework, allows decision makers to make informed choices under uncertainty while being able to incorporate information on the the fly, usually in the form of data, from multiple sources while being consistent with the physical intuition about the problem. This is a major advantage that Bayesian methods bring to fruition especially in the industrial context. This paper is a compendium of the Bayesian modeling methodology that is being consistently  developed at GE Research. 
The methodology, called GE's Bayesian Hybrid Modeling (GEBHM), is a probabilistic modeling method, based on the Kennedy and O'Hagan framework, that has been continuously scaled-up and industrialized over several years.
In this work, we explain the various advancements in GEBHM's methods and demonstrate their impact on several challenging industrial problems. 













§ INTRODUCTION


Reliability of industrial applications necessitates physical experiments of complex engineering systems. 
In conjunction with their cost, measurements  generated from these complex physical tests can be limited, or even sometimes infeasible, in practice due to lack of accessibility and repeatability. 
Advances in numerical algorithms and the availability of high-performance computing platforms, on the other hand, facilitate solution of complex engineering systems with less cost, in general, than the physical experiments. 
For credible numerical simulations, computers models should be verified, validated and calibrated. 
The numerical simulations can only  be  predictive simulations once uncertainties due to model and input parameters are quantified <cit.>.
In abstract sense, both the physical experiments and predictive computer models can be viewed as data generation mechanisms. 
As the complexity of the industrial applications increases, a unified framework linking experimental measurements and numerical simulations  is becoming more critical than ever. 

To this end, the Bayesian paradigm provides a rich probabilistic approach to combine information from  physical experiments and computer simulations in a way that permits engineers to continuously update their prior belief on the process under consideration <cit.>. In Bayesian framework, a probabilistic meta-model can be constructed utilizing data from physical experiments as well as low and high fidelity  computer simulations. Optimally, design of computer experiments (DoE) is performed to guide the construction of such meta-models. The DoE algorithms seek to balance between the cost and accuracy of the computer simulations <cit.>. Industrial-scale computer simulations often have large input dimensions. The input space may include control variables that are set to control the desired process and physical parameters reflecting the media of operation <cit.>. Therefore, sensitivity analysis can estimate the relative influence of each input parameter on the variability of the prediction. The most influential parameters are then calibrated under uncertainty for accurate representation of important aspects of the physical process <cit.>. Once calibrated, uncertainties due to modeling errors and data are propagated  to predict the quantities of interest with confidence bounds.


For the Bayesian framework, building surrogate models utilizing both high and low-fidelity  simulations as well as experimental data are often based on Gaussian Processes (GP) which has become an indispensable tool in the Bayesian paradigm <cit.>.  In this direction, GP surrogates have been applied for various engineering problems <cit.>, such as optimization of composite laminates<cit.>, structural optimization<cit.>, combustion synthesis of composite materials <cit.>, design applications <cit.> and reliability analysis <cit.>, and additive manufacturing <cit.>, to name a few.


The Kennedy and O'Hagan approach <cit.> can be viewed as one of the well established unified frameworks for Bayesian calibration, meta-modeling utilizing both simulation model and observed data. A computer  implementation (named gpmsa) to this approach was released by Los Alamos National Lab  <cit.>. Deploying gpmsa to challenging industrial problems is not a straightforward task, as many technical challenges exist, such as the curse of dimensionality, model-identifiability, transient problems with multiple outputs, lack of data and extrapolation, multiple sets of experimental data, etc. To overcome these challenges, we at GE Research have built the GE Bayesian Hybrid Modeling (GEBHM) framework (shown in Fig. <ref>) to carry out Bayesian analysis for industrial problems. The core components of GEBHM are based on gpmsa, but with several significant enhancements on top that are required for robust employment to many industrial problems.

Presented in Fig. <ref>, the distinguished features of GEBHM are: (1)  using full-Bayesian approach (2) an explicitly formulated model discrepancy (3) building the GP surrogates during the calibration and updating process (4) calibration of  dynamic models with multiple outputs (5) forward uncertainty quantification and propagation (6) model validation (7) real time visualization of the calibration process and outcomes (8) scaling to high dimensional problems (i.e., hundreds of parameters) (8) parallelization  on high-performance computing platforms. These capabilities of GEBHM have been successfully demonstrated on several challenging engineering applications, such as multi-stage engine Blade Row Models, transient model calibration , engine system-level thermal model calibration and validation, combined cycle plant (gas turbine and steam turbine), control model calibration, combustor global sensitivity, predictive emissions modeling and optimization, airfoil cooling system uncertainty quantification and model validation, engine cycle deck performance data matching, Airframe Digital Twin, and more.


 
In this work, we present the  GEBHM  framework and demonstrate its capability in tackling many industrial applications. In particular, Section <ref> briefly reviews the theoretical background of the framework and an application on structural dynamic problem. In Section <ref> we show GEBHM applicability to transient problem, robust design, sensitivity analysis and portable GP, multi-source legacy data modeling, and parallel MCMC. Section <ref> we discuss some of our future enhancements to GEBHM. 



 
 




§ GE'S BAYESIAN HYBRID MODELING (GEBHM)


 


Before delving into the advancement, we will first  discuss the the Kennedy and O'Hagan theoretical framework.

In this section, we discuss the Kennedy and O'Hagan theoretical framework framework for data assimilation. In this framework, let y denote the outputs from either a physical process or a high fidelity simulator, and x represents a design vector of p dimension. The training  output y_i, at training  input point x_i can be expressed as:

    y_i±ϵ_ij(x_i)=η(x_i,θ^*) + δ(x_i),    i=1,  …  ,n,

where  n is the number of training data points, η is the simulation model, θ^* are the true values of calibrated parameters, δ is the discrepancy model between the calibrated  model and experimentally observed data, and ϵ is noise in the observed data.
Generally, the data from the simulation model are available on a set of m design and calibration parameter combination η(x_j,θ_j) for j=1,  …  ,m, which may not be  at the same setting of experimentally observed data.  
As proposed by Kennedy and O'Hagan <cit.>, and as described by Higdon et al. <cit.>, the simulation output and the model discrepancy are modeled as GP models. 
The simulator model (η) is approximated as a GP model using known simulator output at m design locations  with a zero mean and covariance matrix given by,

    Σ_ij^η  = 1/λ_η_zexp( β_η  | x^sim_i - x^sim_j |^2 ) + I 1/λ_η_s    i,j = 1,…, m,

where x^sim is the combined vector of design variables and calibration parameter (x^sim_j = (x_j,θ_j)), the parameters λ_η_z and λ_η_s characterize the marginal data variance captured by the model and by the residuals, respectively, and β_η characterizes the strength of dependence of the outputs on the design variables.
The squared exponential function ensures that the GP model is smooth and is infinitely differentiable. 
Similarly, the experimental observation is modeled as zero mean and covariance matrix given by,

    Σ_ij^y  = 1/λ_y_zexp( β_y  | x^obs_i - x^obs_j |^2 ) + I 1/λ_y_s    i,j = 1,…,  n,

where x^obs is the combined vector of observed design variables of experiments and calibration parameters (x^obs_i = (x_i,θ_i^*)) and θ_i^* is the prior values of calibration parameter. The parameters λ_y_z, λ_η_s, and β_y are for observation GP model.
To represent the correlation between simulator output and experimental observations, the cross covariance matrix is given as


    Σ_ij^η y  = 1/λ_y_zexp( β_y  | x^obs_i - x^sim_j |^2 ) + I 1/λ_y_s
          i = 1,… n, j = 1,…, m.

The discrepancy parameter δ is modeled as 

    Σ_ij^δ  = 1/λ_δ_zexp( β_δ  | x_i - x_j |^2 ) 
           i,j = 1,…, n

where, λ_δ_z and β_δ are the parameters associated with discrepancy model GP. 

Instead of operating on the raw data from the two sources namely, a) experimental, and b) simulation, we scale the observations to a standard normal using the mean and standard deviation of the simulation data. This is done in order to remain consistent with the zero mean function of the GP based models for the simulation and the discrepancy terms. Another important feature of the data-processing phase is the transformation of the experimental normalized observed data into two parts: a) underlying function and, b) discrepancy term. The transformed simulation data is denoted by y_s.
The underlying function part from the transformed experimental data is denoted by y_o and the discrepancy part by y_d.
This is done using a linear basis transformation, the details of which are more involved and outside the scope of the discussion here.
More interested readers are referred to Section 3.2  of <cit.> where the aforementioned data preprocessing steps are exhaustively discussed.
Intuitively, the (y_o^T, y_s^T) and  y_d^T, both of which are assumed to be samples from multivariate zero-mean Gaussian distributions, serve as the training data for the simulation and the discrepancy GP models, respectively.
 
The likelihood of combined data z = (y_o^T, y_s^T, y_d^T) is then given as follows:

    𝕃(D|λ_η_z, β_η, λ_η_s,λ_δ_z,β_δ) = 1/|Σ|^1/2exp( - 1/2 D^T Σ ^-1 D )

where D = (y_o^T, y_s^T, y_d^T) and Σ = [    Σ^y_d        0        0;        0    Σ^y_o Σ^y_oy_s;        0 Σ^y_sy_o    Σ^y_s ]
For multiple outputs y, additional processing such as Singular Value Decomposition (SVD), support kernels for computing δ, etc. are required for which readers can refer to  <cit.>.
The posterior distribution of all the hyper-parameters is given by 

    π (λ_η_z, β_η, λ_η_s,λ_δ_z,β_δ | D) = 
       𝕃(D|λ_η_z, β_η, λ_η_s,λ_δ_z,β_δ) π(λ_η_z) π(β_η) π(λ_η_s) π(λ_δ_z) π(β_δ)

where π(.) on the right-hand side of the equation is the prior distribution of the parameters, where all the hyperparameters are assumed to be independent. The target posterior distribution is evaluated using Metropolis-Hastings within the Gibbs algorithm with univariate proposal distributions for the MCMC posterior updates <cit.>.
The covariance matrices are updated with current realizations of the hyperparameters at every MCMC step.

It is important to emphasize  a fundamental principle of the mathematical formulation in the calibration framework, which is based on the so-called Kennedy O'Hagan formulation <cit.>. This formulation models the simulation model, η(x, θ), as a GP regression model using all the available data i.e. simulation data and observation data. 
This is why the likelihood formulation, in Eq. <ref>, has a cross-covariance between the observed data and simulation data.  This approach is different from other model calibration methods that directly calibrate the known physical parameters of the simulation model using the observed data.



For a test input point matrix X_*, the predicted mean and standard deviation for the gth output is given as,

    μ(X_*) = Σ(X_*, X)Σ(X, X)^-1y   
       Σ(X_*, X_*) = Σ(X_*, X_*) - Σ(X_*, X) Σ(X, X) ^-1Σ(X, X_*)
 
where X is the observed data matrix and y is the vector the observed output used for training the model.



  §.§.§ Application on Structural Dynamics Problem


This section presents a new demonstration for the GEBHM calibration capability. In particular, we consider a physical experiment to model the underplatform dampers used widely in turbine blades and investigate the damping effectiveness of different damper geometries over a wide range of contact preload and dynamic excitation levels. Fig. <ref> illustrates the experiment setup, which consists of two rectangular beams placed adjacent to each other. The beams have a platform at roughly one third of the beam length from the fixed end. A vibration damper is placed between them such that it contacts the underside of the beam platforms. The damper is held in place by two threaded rods which are attached to surrounding stationary frame. The damper can be loaded up against the platforms by pretensioning the rods and this load is measured by loadcells between the pretensioning nuts and the stationary frame. Dynamic excitation is provided to one of blades using an electrodynamic shaker.  Analytical simulation of the experiment is performed by generating a Craig-Bampton type reduced order model from the finite element model of the test setup shown in Fig. <ref> and performing a nonlinear vibration analysis using a multi-harmonic balance hybrid frequency time domain solver similar to that described by Poudou et al. <cit.> This helps improve the solution times by several orders of magnitude over the traditional approach of using a full FEA model with a transient time stepping analysis.

 
The problem has two design variables: 1) the static damper preload and 2) the dynamic excitation. Furthermore, there are  two unknown calibration parameters that are 1) the friction coefficient (μ) and 2) backround structural damping which is applied as purely stiffness proportional damping with the proportionality constant β.  The design variables preload and excitation are denoted by x_1 and x_2, respectively.

The friction coefficient μ and the background damping β are the parameters being inferred. Both the test and simulation provide response as a function of excitation frequency. The frequency responses are parameterized by two parameters. The maximum response and the damping expressed as Q-factor. The former captures the height of the response peaks while the latter captures the width of the frequency response function peaks. We make use of the test data obtained for 24 designs  and simulation data at 30 designs. calibration_vibration. shows the posterior densities of the calibration parameters μ and β obtained using GEBHM calibration framework. Whereby, GEBHM estimates the posterior distributions by fusing information from the two sources and using the fully-Bayesian MCMC discussed in gebhm_calibration.





The nonlinear vibration analysis is run at the test conditions using the mean values of the calibrated parameters μ and β. Fig. <ref> shows the comparison of the calibrated model results and the test data at two different preloads (50 lbf and 100 lbf). For the 50 lbf case, analytical responses obtained with the calibrated parameters match test data well in terms of overall frequency response for three out of four dynamic excitation levels (0.05, 0.1, 0.5 lbf). For the  0.2 lbf excitation the analytical response overpredicts the test data by 30%. The frequency at which the maximum response occurs is predicted accurately for all excitation levels. For the the 100 lbf case, good matches on the maximum response are obtained for three out of four forcing levels. For the 0.5 lbf forcing level the model overpredicts the max response by 37%. In all other cases, maximum response prediction errors are less than 5%. It is possible to draw samples from the distribution of μ and β and run the forced response analysis on these samples to generate confidence bounds on model predictions if the computational cost does not become prohibitive. This is illustrated in  Fig. <ref>(b) where it can be seen that most of the test data falls within 95% model prediction bounds. Overall, with the calibrated parameters good match between data and model was observed over a wide range of preload and dynamic excitation levels.

 




§ ADVANCES AND APPLICATIONS




 §.§ Transient Response


In this section, we incorporate the ability to handle multiple transient responses in the GEBHM framework. We only present a brief introduction to the capability here and refer the interested reader to <cit.> for more information on the theory and industrial applications. Similar to the original GEBHM framework, we utilized Higdon et al's <cit.> implementation of the transient calibration capability. Instead of dealing with the outputs directly in the original GEBHM implementation, we initially perform a singular value decomposition (SVD) of the outputs to identify the principal components to reduce the overall dimensionality of dealing with large transients variables, as follows:

    η(X,θ,t)_n_t× m = U_n_t× n_tΣ_n_t× mV^T_m× m,


where n_t are the number of  time-points in the transient outputs, m are the number of simulations, and θ are the calibration parameters of the simulation model η, U and V are the left and right singular vectors, and Σ is the diagonal matrix of singular values. Then we retain only certain number of principal components (say n_p_u) to approximate the transient simulation model as,
 
    η(X,θ,t)_n_t× m≈ U_n_t× n_p_uΣ_n_p_u× n_p_uV^T_n_p_u× m,

which can then be further linearized as,

    η(X,θ,t) ≈  K_simw,

where K_sim is UΣ and w is K^-1_simη. Next we build a GP model for the principal components w term. The simulation output, the experimental output can also be linearized as,

    y(X,t) ≈  K_obsu,

where K_obs is computed by interpolating K_sim from the simulation time grid onto the observed time grid. This ensures generality of the overall method such that the simulation and observed data are measured at different points in time. The components of the observed data  is modeled as a GP with covariance matrix as shown in Eq. <ref>.

In order to model the discrepancy between the simulation model and observations for the transient outputs, we utilize a specific number of Gaussian basis for the linear model as,


    δ(X,θ,t) = Δ(t) v(X),









 
Discrepancy components v are modeled again as a GP similar to that shown in Eq. <ref>. Essentially, we have utilized linearization  to represent the transient problem using principal components for the simulation model, and Gaussian basis components for the discrepancy model in order to develop the GP models. Once these covariance functions are constructed, the hyperparameters can be inferred in a similar manner as described in the earlier sections. We have enhanced Los Alamos's original implementation to handle multiple transient outputs and have demonstrated good success in applying it to industrial applications, as shown in <cit.>.

 



  §.§.§ Application: Transient Response


We provide a brief description of one of the problems addressed in <cit.>, which is a transient heat transfer problem on a 1D rod. Let us consider a rod of length L, as shown in Fig. <ref> with a heat flux on the left end and an initial temperature condition. Our aim here is to calibrate the heat flux and initial temperature based on the temperatures monitored on five different locations along the length rod. Fig. <ref>(a) shows the calibrated posterior distributions of heat flux and the initial temperatures - we can clearly see that the GEBHM method can identify the strong correlation between the flux and initial temperature clearly. 
 
Figure <ref>(b) shows the transient temperature signals with the discrepancy correction.  We have summarized the multiple transient calibration capability that is part of the overall GEBHM framework here, and interested readers are referred to our publication here: <cit.>.


 §.§ Robust Gaussian Process

In many industrial problem, the input x_* cannot be observed directly and the uncertainty around it is modeled with a Gaussian distribution as x_* ∼𝒩(u,S), with mean u and covariance S.
In such a scenario is it important to quantify the effect of input uncertainty on the output, which is generally done by uncertainty propagation using a meta-model such as a GP.
The resulting predictive distribution of the the output, p(y_*), when the input is a random variable x_* given by p(x_*|u,S), is obtained by marginalizing over the input distribution:

    p(y_*|u,S,𝔻)= ∫ p(y_* | x_*,𝔻) p(x_*|u,S)dx_*   , 
       
        p(y_*| x_*,𝔻)= 𝒩(y̅_*,σ_y_*)   ,

where 𝔻 is the training data, y̅_* and σ_y_* are the predicted mean and standard deviation of output. 

Traditionally, marginalization is estimated by Monte-Carlo sampling approach. The sampling approach using GP is fast, however, in a scenario like uncertainty-based design optimization, the double-loop process (optimization loop and sampling loop) can make the overall process expensive when the system has a large number of random input variables.
Robust GP using GEBHM overcomes this challenge by collapsing the double-loop uncertainty based design optimization problem into a single-loop problem. This is accomplished by exactly calculating the uncertainty measures E_x[y_*] and V_x[y_*]. The robust GP can directly, with one function evaluation, estimate the first and second moments of the output with the consideration  the effect of input uncertainty around. An example of this distinction is shown in Fig. <ref>.  Note y̅_* and σ_y_* are probabilistic GP model predictions and cover model uncertainty only, not the added uncertainty due to uncertain inputs.  


Assuming that x_* is a normally distributed random variable, GEBHM can be used to obtain exact analytical expressions for the mean, m(u,S), and variance, v(u,S), of the marginalized predictive distribution of output. 
This can enable the formation of single-level uncertainty-based optimization problems, where the objective is some function of m(u,S) and/or v(u,S).

The method to solve m(u,S) and/or v(u,S) implemented in GEBHM is based on work developed by Quinonero-Candela, et al. <cit.>.
Using the law of iterated expectations, and the law of total variance:

    m(u,S)    = E_x[E_y[y|x]]=E_x[y̅]  , 
       
    
    v(u,S)    = E_x[V_y[y|x]] + V_x[E_y[y|x]]   , 
       = E_x[σ_y^2]+V_x[y̅]   ,

where E_x[·], V_x[·] denote the expectation and variance wrt. x, respectively. When using Gaussian kernels in GPs, the expressions for y̅ are Gaussian  functions of x and the expressions for σ_y^2  are products of Gaussian  functions of x.  Therefore, the integrand involved in determining m(u,S) and v(u,S) are products of Gaussian  functions, which allows for an analytical calculation.
The analytical form of the mean is given as: 

    m(u,S)    = Ψ^⊤l  , 
    Ψ   = (Σ^η)^-1y .

Elements of vector l=[l_1,…,l_m]^⊤ are given by,

    l_j =    1/λ_η_z2(𝐁)^-1S+I^-1/2
       ×exp(-(u-x_j)^⊤(2S+𝐁)^-1(u-x_j))+ I1/λ_η_s  ,

where 𝐁 = diag[β_y1^-1,…,β_yp^-1] and I is a m by m identity matrix. The variance is given by,

    v(u,S) =    Σ_**^η-(((Σ^η)^-1-ΨΨ^⊤)L) -(ll^⊤ΨΨ^⊤)   ,

where, Σ_**^η is the covariance matrix constructed using <ref> for new design point x_*,  Σ_**^η = Σ^η (x_*,x_*). (·) is the trace operator and L is a m by m matrix with elements given by the following:

    L_i,j   =(1/λ_η_z)^24(𝐁)^-1S+I^-1/2
       ×exp(-[(u-x_d)^⊤(2S+𝐁/2)^-1(u-x_d).. 
       +..(x_i-x_j)^⊤(2𝐁)^-1(x_i-x_j) ])+ I1/λ_η_s,

where we define x_d = 1/2(x_i+x_j). 




  §.§.§ Application: Robust Gaussian Process

An application presented here is a customized subset of the proposed NASA UQ Challenge <cit.>.
The optimization problem is to minimize the instability with a constraint on the magnitude of the instability uncertainty.
The instability is a function of four Gaussian random variables (x = [d_6,d_13,p_4,p_5]). The design variables are the respective mean values of (d_6,d_13,p_4,p_5). For optimization, a Genetic Algorithm (GA) was used with a population of 100, and set to run for 10 generations.  At each optimization step, the single-loop approach calculates E_x[] and V_x[] directly from a robust GP with a single function evaluation. In comparison, with the double-loop approach, uncertainty metrics are calculated using 1000 Monte-Carlo samples from a standard GP. The time-history of GA's iterations for the double-loop and single-loop approaches, are provided in Figure <ref>.  The robust optimum was found between generations 3 and 4 for both approaches.  The single-loop approach takes 336 function calls, while the double-loop approach reaches the minimum instability value after 347,000 function calls. The  double-loop approach costs about 1000 times more than the single-loop approach  to find the robust optimum. For more details on theory and the application on robust GP, please refer to <cit.>.














 


















 §.§ Bayesian Sobol Sensitivity Analysis


Design and analysis of complex engineering systems often presses for a 
thorough know-how of the relative order of importance of the 
various inputs. 
In this work, the model of an engineering system could be a physical experiment 
and/or a computer simulation which can be evaluated at a specific set of input conditions.
Such models are known as expensive black-box functions.
Additional challenges usually associated with such models are a high-dimensional input space
and a lack of gradient of information.
The inputs to the model can include some design variable(s) and some unknown parametric 
constant(s) of the model.
Knowledge of the relative sensitivity of the inputs is critical for designers
while making decisions under uncertainty, as it helps in expediting design exploration, 
identifying a set of critical variables and making simulation decisions under limited-data.









The variance-based methods, commonly known as the Bayesian Sobol sensitivity analysis (BSSA),
are more flexible to undertake the task of obtaining sensitivity information
for high dimensional problems with limited data under uncertainty.
The BSSA investigates how the uncertainty in the black-box model of a system can be divided among the different
uncertain inputs and their interactions.  An addendum to the above facet of the Sobol sensitivities is the ability 
to infer sensitivity of the underlying physical response to higher-order interactions.
All expressions for the main-order and higher-order effects are
computed using analytical expressions. 













We assume that the inputs are uncorrelated and have the following distributions.

    x_i∼𝒰(0,1) ∀ i ∈ d.

where d is the dimensionality of the input space and 𝒰(0,1) is a standard uniform distribution. Since we do not have full access to the true response, we resort to using a probabilistic 
representation of y() as described in koh_y. Denoting the posterior predictive mean of the predictive simulator by η(; ), where   represents the  design variable(s) and  denotes the calibration parameter(s), we make use of the high-dimensional model representation (HDMR) <cit.> to represent y() as follows:


    y()    =η(; )
       =z_0 + ∑_p=1^dz_p(x_p) + ∑_p<qz_p, q(x_p, x_q)
        + z_p,q,⋯,d(x_p, x_q,⋯,x_d),


where,

    z_0 = 𝔼(Y),  
       z_p = 𝔼(Y|x_p) - 𝔼(Y), 
       z_p, q = 𝔼(Y|x_p, q) - z_p - z_q -  𝔼(Y), 
       z_p, q, r = 𝔼(Y|x_p, q, r) - z_p, q - z_q, r - z_r, p - z_p - z_q - z_r
    -𝔼(Y).


To avoid confusion between design variables and calibration parameters we concatenate  and  to represent the set of variables as ={, }.
This is consistent with the formulation since the BSSA follows the same analytic formulation for both sets of variables.
The  terms denoted by z in hdmr are known as effect functions and are orthogonal under the assumption that the distribution over the input space is uniform and that the true response is square-integrable <cit.>.
This  property simplifies the derivation and computation of the BSSA indices for the scenario with uncorrelated inputs.

    𝔼(Y|_p) = ∫_𝒳_-pη()_-p|p(_-p|_p),



    S_p = var{𝔼(Y|_p)}/var{Y},
 


    S_pq = var{𝔼(Y|_pq)}/var{Y}  - S_p - S_q.
 

Analytical tractability of the Gaussian form of η(), derived in <cit.>, allows one to further represent the BSSA indices as follows:


    S_p = 𝔼^*var{𝔼(Y|_p)}/𝔼^*var{Y}
 
and two-way interactions as:

    S_pq = 𝔼^*var{𝔼(Y|_pq)}/𝔼^*var{Y}  - S_p - S_q
 
where, 𝔼^* is the expectation with respect to the probabilistic surrogate model.

Similarly, we can define the BSSA for higher-order interactions as:

    S_pqr = 𝔼^*var{𝔼(Y|_pqr)}/𝔼^*var{Y}  - S_pq - S_qr
    - S_rp - S_q - S_p - S_r





  §.§.§ Application: Bayesian Sobol Sensitivity Analysis


The structural dynamics problem discussed in the previous section is used as a test-bed to demonstrate the BSSA capability of GEBHM. The BSSA obtained with respect to the design variables and the calibration parameters are shown in bssa_vibration.







 §.§ Bayesian Sobol Sensitivity Analysis with Correlated Inputs


The work discussed in the previous sections assumes that inputs are uncorrelated.
In that case, the HDMR has orthogonal effect functions, and
one may compute each Sobol index independently.
Presence of correlated inputs mandates sampling from joint distributions of the inputs which increases the computational cost tremendously.
At the same time,  ignoring correlations can have a detrimental effect on the inferred BSSA indices as has been highlighted in <cit.>.

Saltelli et. al. <cit.> and other authors
have described this effect to be carried over due to correlation
and stress the need to be careful when making conclusions
with Sobol indices when inputs are correlated.
The BSSA approach developed in Xu et. al. <cit.> splits the 
contribution of an individual input
to the uncertainty of the model output into two components: correlated and uncorrelated.
This approach was carried forward by 
Li et al. <cit.> with a
relaxed HDMR concept to compute variance contribution
due to correlation and problem structure for non-linear models.
Chastaing et al. <cit.> develop a hybrid method for decomposing the variance by utilizing the  hierarchical
orthogonality of component functions and the
Gram–Schmidt method.
A hybrid method using polynomial chaos expansions and copula method was proposed in <cit.>.

We proceed with the framework proposed by Li and Rabitz <cit.> which enforces orthogonality only on the lower-order component functions, a technique commonly known as hierarchical orthogonality.
In other words, z_pq(x_p, x_q) is only required to
be orthogonal to z_p(x_p) and z_q(x_q) and not to any other component functions. 
Following the above approach, the output variance can be decomposed as follows:

    var{Y} = ∑_=1^2^d-1var{z_} + cov{z_, ∑_=1^2^d-1 z_}


Here,  and  are multi-indices that denote the different interaction terms.
The total index S_ is only a sum of the
structural and correlative contribution.


    S_^a   = var{z_}/var{Y},
    
    S_^b   = cov{z_, ∑_=1^2^d-1 z_}/var{Y},
 
where S_ = S_^a + S_^b. Note that, the structural and correlative component functions cannot
be estimated independently of each other in the presence of input correlation and more advanced methods would be needed to infer the BSSA indices for additive models.

The use of GEBHM's BSSA with uncorrelated and correlated inputs has been successfully demonstrated on several industrial applications<cit.>.
The usefulness of this framework is best highlighted in problems with reasonably high-dimensionality and limited observed or simulation data under uncertainty, while being fully integrated with GEBHM's model calibration framework.
We demonstrate the approach on industrial applications in the following sections.



  §.§.§ Application: Bayesian Sobol Sensitivity Analysis with Correlated Inputs




To highlight the impact of the correlated BSSA capability of GEBHM we apply the framework on a problem of engine vibration.
The problem has one output and six design variables.
The correlation exists between inputs denoted by x_1 and x_4. 
The impact of including the correlative term on the overall ordering of the BSSA is highlighted in bssa_vibration_corr (a) and (b).
The  structural BSSA is unable to incorporate the effect of the correlation between the two.
However, the final BSSA, with the correlative part included, reorders some of the interaction effects involving x_1 and x_4.

 



 §.§ Portable Gaussian Process


It is often the case that a small number of effects as identified in Sec. <ref> explain the majority of the GEBHM metamodel's predictive trends.
In addition, we have encountered a number of applications where the m^2 storage footprint and m computational cost associated with computing predictions with the full GEBHM make it infeasible to apply.
In this section, we turn our attention to a method for approximating the posterior mean m of the GEBHM model by leveraging the observation that it may be approximated by the combination of a relatively small number of effect functions;
each of these effect functions is approximated through a generalized linear model.
Doing so eliminates the undesirable nonparameteric characteristics of the GEBHM while incurring a negligible penalty in predictive performance.
Due to the lightweight nature of this model, we refer to it as a  portable Gaussian process  or  portable Bayesian Hybrid Model  (pBHM)<cit.>.
The remainder of this section describes the methodology by which one constructs a pBHM.

The approach takes as input a trained GP model from GEBHM and the set Ω_x ⊂^p over which one expects to make predictions, equipped with a density p().
We consider a scalar posterior mean η(·) of the GP model, recognizing that one may create a multi-output model by repeating the process outlined in this section for the other N output dimensions.
The output of the approach is a parametric approximation of η(·) that may be used within Ω_x.

Carrying out the sensitivity analysis of Sec. <ref>, we compute the Sobol indices of the HDMR of Eq. (<ref>) and identify those whose magnitude are greatest; let the set of Sobol indices be represented by = { S_p }, where S_p ⊂{1, …, p} is the set of dimensions associated with an index.
We retain the effect functions corresponding to the largest Sobol indices, picking the number of terms such that the sum of the associated indices exceeds some threshold (typically 99%), ensuring that this truncated series captures the majority of the variation in η over p().
Note that if the GP kernel function is of the form in Eqn. <ref> and p() factorizes into uniform or Gaussian marginals in each of the input dimensions, then the Sobol indices may be computed analytically.

As shown in Sec. <ref>, if k(·, ·) is parametric with an exponentiated quadratic form [e.g. as in Eq. (<ref>)] and p() is either uniform or Gaussian, then the effect functions may be computed analytically.
We now consider the task of fitting a parametric model to a given effect function z_S_p.
The following task is conducted independently for all S_p contained in  and may therefore be trivially parallelized.
We sample m̃ points from some experimental design over the input space in the dimensions contained in S_p.
This may come from a space-filling design (e.g. Latin Hypercube) or simply by sampling from the marginal p(_S_p).
We collect the inputs 
∈^m̃× |S_p| 
and corresponding outputs under the effect function 
∈^m̃ 
as a training dataset 
= {^(i), ỹ^(i)}_i=1^m̃.
We seek to learn a generalized linear model 

    z_S_p≈_S_p_S_p(_S_p),

where _S_p contains a set of basis functions and α_S_p the corresponding coefficients.
We have used polynomial basis functions as well as splines in previous work <cit.>.
We determine _S_p by minimizing a mean squared error loss on  with an additional L_1 penalty to induce sparsity in _S_p;
the magnitude of the regularization hyperparameter is tuned by cross-validation. 
Note again that it is comparably cheap to generate data from the effect function to be approximated, so cross-validation is not subject to troubles arising from data scarcity.
























  §.§.§ Application: Portable Gaussian Process

Here, we demonstrate the application of portable Gaussian process or pBHM to model crack growth for an airframe component.
We focus on constructing a surrogate model for a stress intensity factor (SIF), a key quantity driving crack growth, as a function of load and geometry. 
Because the surrogate will be used in a many-query setting, quick predictions are highly desired.

A GP model was trained using GEBHM on m=360 data generated via finite element simulations, parameterized by p=7 inputs.
A sensitivity analysis using GEBHM reveals that third- and higher-order interactions are negligible.
We consider several parametric fits for the effect functions, eventually choosing spline fits.
Fig. <ref> shows the fit obtained for a representative main effect function.
 


The resulting spline model is tested on 360 held-out validation data.  
The percent errors of the original GP model from GEBHM model and the spline version are shown in Fig. <ref>.  
The mean absolute error is 1.9% for the GP and 2.9% for the spline model. 
The root-mean-square-error is 5.2% for the GP and 11.9% for the spline model. 
In exchange for a modest decrease in accuracy, the the spline model enjoys a 50x speedup over the GP model.





 §.§ Multi-Source Legacy Data Modeling


For a new industrial design, typically only sparse data is available due to cost and time associated with high-fidelity analysis and experimentations. For such complex systems, sparse data may not be sufficient to build an accurate model of design performance. Typically, a multi-fidelity approach is  sought when one can generate relatively-abundant data from a cheaper lower-fidelity model and integrate it with the sparse high-fidelity data. However, in an industrial setting, high-fidelity data are generally available from legacy designs which are not exactly same as the new design, but they belong to a similar family. In this case, we extend the capability to handle data from multiple legacy sources to improve the predictive capability of the new system. This enhancement can be applied even if low-fidelity data are available from legacy systems.

In multi-source legacy modeling with GEBHM, separate models are built for each legacy data by combining the legacy data with new system data.
For k^th legacy data, the legacy model is built as 𝕄_k(x) = η_k^legacy(x) + δ_k^legacy(x), by treating the data in k^th legacy data source as simulation data and new system data as experimental or high-fidelity data.
It is assumed that the input variables (design, operational, etc.) and the output variables (performance, cost, etc.)  are the same for the new system and all the legacy systems.
In this scenario, when legacy-data has more variables than the new system, i.e. dim(x_legacy) > dim(x_new), then the extra variables in legacy data can be treated as calibration parameters.

If K legacy models are generated, then the multi-source prediction at new location is evaluated as 

    y(x) = ∑_k=1^Kw_k(x) (η_k^legacy(x) + δ_k^legacy(x))

where w_k is the model validity which is a function of the input variables such that ∑_k=1^Kw_k(x)=1 for any value of x.

The model validity associated with a legacy model signifies relative confidence of legacy models at a given design location x. There can be multiple ways to evaluate model validity like expert opinion, Bayes' factor, relative age of the legacy data, etc. 
In GEBHM, we use two statistical metrics to evaluate model validity; likelihood-based model validity and  uncertainty-based model validity.
In the likelihood-based model validity, the validity of a legacy model at a input setting x_new^* where data for the new design is available is proportional to the probability of the legacy model to predict the output or performance of the new design.
For example, the likelihood of the k^th legacy model, 𝕄_k at a location x_new^*, for which output of new design is known to be y_new is given as

    w_k^lik(x_new^*) ∝1/√(2 πσ_k^2(x_new^*))exp(- ( y^* - μ_k(x_new^*))^2/2 σ_k^2(x_new^*))

where μ_k(x_new^*) and σ_k(x_new^*) are the predictive mean and standard deviation of the k^th legacy model 𝕄_k at x_new^*.

In the uncertainty-based model validity, the model validity at x_new^* is inversely proportional to the predictive uncertainty as

    w_k^unc(x_new^*) ∝1/σ_k^2(x_new^*)


The overall model validity of the k^th legacy model at a given design point σ_k^2(x_new^*) is then given as

    w_k(x_new^*) = κ(x_new^*) w_k^lik(x_new^*) w_k^unc(x_new^*).

where κ(x_new^*) is a proportionality factor which is constant across all legacy model at a given x_new^*.
GP models of model validity,  w_k(x),  are then built for each legacy model by using the overall model validity evaluated at design points x_new at which new system data are available. 
For more details, please refer to <cit.>.



  §.§.§ Application: Multi-Source Legacy Data Modeling

In one of the application, we have used multi-source legacy data modeling to predict pump efficiency for a new pump design.
These pumps are required for heavy duty services in refinery, petro-chemical and other industrial services, pumping fluids under wide range of pressures and particularly at high temperatures.
The goal is to build a predictive model of a maximum efficiency of new pump design as a function of nine input parameters such as blade angle, eye diameter, hub diameter, pump speed, number of impeller blades, etc. 
For the new pump design we have 18 performance data and different setting of input parameters. 
We also have data available from five legacy pumps, with number of data varying between 10 to 55.
We used all these data to build a multi-source legacy data model for the new design. 
We carried out 10-fold cross validation to demonstrate the improvement in accuracy of using multi-source legacy data modeling over other models.
Fig. <ref> shows the cross-validation error against standard error of cross-validation for different models. 
The “New-Sys" is the model built on new system data only, “Legacy-k"  (for k ∈ 1, …, 5) are model based on multi-fidelity modeling by combining k^th pump legacy data with new system data, and “MS-Legacy" is the model based on multi-source legacy modeling. 
As seen, using multi-fidelity modeling using data from legacy pump-3 or pump-4 improves the model prediction when compared to model just build on new system data. 
However, using multi-source modeling using all the legacy data, improves the model significantly.





 §.§ Parallelizing GEBHM's MCMC


While the GP does not assume a functional form of the problem, it is defined by a set of parameters, called hyperparameters, that need to be learned from the data.
The hyperparameters define the characteristics of the objective function, such as smoothness, magnitude, periodicity, etc. Accurately estimating these hyperparameters is a key ingredient in developing a reliable and generalizable surrogate model. 
Markov chain Monte Carlo (MCMC) is a ubiquitously used Bayesian method to estimate these hyperparameters. 
GEBHM is very effective on problems of small and medium size, typically less than 1000 training points.
However, traditional GP models do not scale well with respect to dimension and size of the dataset.
For some challenging industry applications, the predictive capability of the GP is required but each second during the training of the GP costs thousands of dollars.  
In this work, we apply a scalable MCMC-based methodology enabling the modeling of large-scale industry problems. 
To this end, we extend and implement in GEBHM an Adaptive Sequential Monte Carlo (ASMC) methodology for training the GP. This implementation saves computational time (especially for large-scale problems) while not sacrificing predictability over the current MCMC implementation. 
For more details on the theoretical workings on the method we refer interested readers to <cit.>.
We demonstrate the effectiveness and accuracy of GEBHM with ASMC on on a two-objective challenging industry application in asmc_gebhm_app.



  §.§.§ Application: Parallelizing GEBHM's MCMC


The problem is a GE combustion test data problem, where we aim to build a GP model to predict two emission quantities as a function of three measured temperatures, three fuel flow parameters, and air flow. Thus, there are nine hyperparameters that need to be estimated for each individual objective. 
The total number of hyperparameters for the GP model is 19. 
The MCMC settings for this problem are: 10000 steps with 1000 steps during the initialization.

 

 

The significance of the results in combustions_1 and combustions_2 is noticeable from the perspective of the number of objectives. Not only does the scaled-up HPC ASMC do well in terms of computational time but it improves the predictive accuracy of the GP model compared to GEBHM's MCMC. 
The number of steps per ASMC particle is one for the workstation and five for the HPC runs. 
The number of sequential distributions has been fixed to 10 for both the computing environments.

This provides a holistic overview of ASMC's performance on different types of problems. Problems of varying input dimensionality, varying number of objectives, and different training data sizes have been treated with the extended ASMC and compared to the current MCMC implementation. The ASMC does equally well or better compared to the MCMC in terms of predictive accuracy, while saving time by half or more in most of the challenging problems.





§ CONCLUSION AND FUTURE DIRECTION


We discussed GE's Bayesian probabilistic modeling framework, GEBHM, and various advanced methods developed for industrial applications. Although we only discussed the core capabilities of GEBHM in this work, these capabilities have also enabled development on other frontiers such as intelligent/adaptive design optimization  <cit.>. The calibration and sensitivity analysis capabilities of GEBHM are demonstrated on a new  challenging industrial problem of structural dynamics as discussed. A summary of the numerous enhancements to the GEBHM framework are provided in Table (<ref>). With the increase in complexity of industrial applications and demand of further improvement in probabilistic modeling methods,  some future directions of enhancement are as follows:







 §.§ Preconditioned Conjugate Gradient Method



The demanding computational cost and memory requirement constitute a significant challenge for GP to train and predict based on large datasets. The complexity stems from the requirement of multiple solutions of a large-scale system of linear equations involving the kernel matrix. Inspired by the concepts of divide-and-conquer, we will develop a scalable GP algorithm based on domain decomposition methods. The proposed algorithm will partition the hard-to-solve problem into smaller easy-to-solve ones which can be tackled concurrently. The solutions of the local problems will be utilized to construct a scalable preconditioner for the preconditioned conjugate gradient method (PCGM) to tackle the global problem. For complex industrial applications, a multi-level and multi-fidelity preconditioner will be essential for cost-effective performance and robustness of the algorithm <cit.>.





 §.§ Scalable Gaussian process


With the rapid development of high-performance computing, Big Data and Internet of Things, it has become necessary to scale up GPs for modeling large dataset. Various frameworks have been proposed to scale up GPs such as sparse Gaussian Process for global approximation and distributed Gaussian Process for local approximation. Most existing schemes focus on speeding up the training while balancing the overall prediction accuracy. However, there are certain cases requiring special treatment. For example, the dataset could be from non-uniform distribution in which the dataset are scattered as clusters among the input space. For modeling large datasets with non-uniform distributions, the scalable schemes might overlook the clusters with relatively less data. This will deteriorate the predicative performance at the overlooked area. It's possible to fix the issue with intelligent selection of a subset for GP development with a reasonable coverage among the whole dataset  <cit.>. The intelligent selection should take into account both input variables and the output values. A small subset is expected to be adequate while assuming the large dataset contains more-than-enough data for model development.    





 §.§ Deep Gaussian processes

Gaussian process models have remarkable data efficiency properties due to the very narrow, yet reasonable, prior beliefs that they place on the function being learned.
A core element of this prior is the GP kernel, which encodes beliefs about the regularity and differentiability of the function being modeled.
However, traditional kernels are often selected out of convenience rather than out of the belief that they accurately reflect one's prior beliefs.
There are many systems that are known to violate familiar assumptions such as stationarity, homoscedasticity, and Gaussian marginal densities.
One way of overcoming these limitations is by learning new representations of the data being modeled to enter as inputs to the Gaussian process simultaneously with training the GP itself.
In this spirit, recent model architectures including deep kernel learning <cit.> and deep Gaussian processes <cit.> are have been shown to improve predictive performance and are becoming increasingly accessible due to the advent of advanced approximate inference techniques <cit.>, software platforms <cit.>, and modern computing hardware.
We see considerable opportunity as the community continues to experiment with combining these techniques with traditional modeling approaches.



 §.§ Modeling multiple datasets

Up to this point, our discussion has focused mainly on modeling a single input-output relationship.
However, there are many general classes of problems inside which small permutations alter the details of the input-output relationship being modeled.
One could imagine that the similar problems previously solved might be used to inform the model for the problem at hand.
We are investigating the use of compositional models in which task labels are embedded as latent variables in multi-dimensional Euclidean space; these latent variables augment the input space, allowing one to simultaneously learn a single model over the augmented input space while learning an embedding of tasks that enables a parsimonious functional representation.
We recognize that while one may assemble very large datasets by combining tasks, the complexity of the input space will grow at the same time.
Thus, there remains a strong case to be made for incorporating uncertainty associated with both the learned embedding as well as the function relating the augmented input space to the outputs.
By following the Bayesian modeling groundwork laid out in the previous sections, we anticipate that one will be able to learn credible models that can be trained with even fewer data than what the BHM framework enables today.



§ ACKNOWLEDGMENTS

The authors acknowledge the contributions from Dr. Ankur Srivastava, Dr. Felipe Viana, Dr. Arun K. Subramaniyan, Dr. Issac Asher, Dr. You Ling, Dr. Kevin M. Ryan, and Dr. Jesper Kristensen for their critical contribution in development and enhancement of GEBHM. 

asmems4




