{'fraction_non_alphanumeric': 0.041987716359575655, 'fraction_numerical': 0.010329424902289224, 'mean_word_length': 4.616494198808404, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fixed Inducing Points Online Bayesian Calibration for Computer Models with an Application to a Scale-Resolving CFD Simulation

This paper proposes a novel fixed inducing points online Bayesian calibration (FIPO-BC) algorithm to efficiently learn the model parameters using a benchmark database. The standard Bayesian calibration (STD-BC) algorithm provides a statistical method to calibrate the parameters of computationally expensive models. However, the STD-BC algorithm scales very badly with the number of data points and lacks online learning capability. The proposed FIPO-BC algorithm greatly improves the computational efficiency and enables the online calibration by executing the calibration on a set of predefined inducing points. To demonstrate the procedure of the FIPO-BC algorithm, two tests are performed, finding the optimal value and exploring the posterior distribution of 1) the parameter in a simple function, and 2) the high-wave number damping factor in a scale-resolving turbulence model (SAS-SST). The results (such as the calibrated model parameter and its posterior distribution) of FIPO-BC with different inducing points are compared to those of STD-BC. It is found that FIPO-BC and STD-BC can provide very similar results, once the predefined set of inducing point in FIPO-BC is sufficiently fine. But, the FIPO-BC algorithm is at least ten times faster than the STD-BC algorithm. Meanwhile, the online feature of the FIPO-BC allows continuous updating of the calibration outputs and potentially reduces the workload on generating the database.


Introduction
Computational modelling is now an integral tool to predict and understand complex physical systems in many branches of science and engineering. Computational modelling enables engineers and scientists to analyse physical phenomena which cannot be easily measured by experiments and reduces the time and financial cost of analysing complex systems. However, a computational model is normally subject to various types of errors, such as discretisation error, parametric uncertainty error, and model error. The discretisation error, parametric uncertainty error and model error all contribute to the total error in a mathematical or computational model. It is only by assessing and reducing each of the components to the total error that one can assess the accuracy of a given mathematical or computational model. The discretisation and parametric uncertainty errors are relatively easy to reduce, depending upon the computational resources and experimental data available. However, the model error is normally linked to the underlying mathematical model not representing the exact physics of the system being analysed. Therefore, model error is particularly challenging to both rigorously analyse and reduce in magnitude. In many applications, such as turbulence modelling, the model parameter is determined based upon approximations or simplified conditions. Therefore, model parameter calibration with respect to (w.r.t.) the benchmark data becomes important to improve the performance of the model.
The determination of the model parameters through calibration, data adjustment or data assimilation enables the model more accurately represent the underlying physics of a system being analysed. By this definition, it is natural to calibrate the model parameters using the best-fit principle. However, the best-fit principle does not consider the uncertainties that the data may be subjected to, such as the observational/numerical discretisation errors, which should be considered in determining the agreement between two databases [1][2][3][4]. Kennedy and O'Hagan [5] proposed a Bayesian approach to calibrating the computational model by representing model bias and computer model outputs as the Gaussian processes. The method automatically integrates the uncertainties into the calibration by treating them as Gaussian noises. For simplicity, we will refer to Kennedy and O'Hagan's approach as standard Bayesian Calibration (STD-BC).
Since the development of STD-BC, the concept was adopted and further developed by many researchers in different areas of science and engineering. For example, Bayarri et al [6][7][8] applied the method to validate models for resistance spot welding, dynamic stress and vehicle collision modelling. In computational fluid dynamic (CFD) simulation, it was applied to tuning the model parameters and inferring the true physical value [9,10], quantifying the inverse uncertainty [11][12][13][14][15][16][17][18] and informing the design of supercomputer simulations [19].
Recently, Wu et al. [20,21] modularised the STD-BC by dividing the database into two different sets; a calibration set and an inverse uncertainty quantification set. This framework has been successfully applied in Liu et al [22] to carry out the inverse uncertainty quantification of the correlation in the two-fluid model based multiphase-CFD simulation.
Wu et al. [23] furtherly discussed the "lack of identifiability" issue of inverse uncertainty quantification using the Bayesian calibration. Karagiannis and Lin [24] proposed a mixture Bayesian Calibration framework, which can represent the real system by mixture the outputs of available computer models with different fidelity level.
Despite the successful applications listed above, the STC-BC is impractical for science and engineering problems involving very large data sets [25]. This is due to the poor scaling of inverting the covariance matrix with respect to (w.r.t.) the size of the data sets involved.
Assuming we have a data set of size N (including numerical and benchmark data), there will be an N-by-N covariance matrix to invert. The computational complexity to invert an N-by-N matrix is N 3 . To solve this problem, Higdon et al. [26,27] followed the definition of Kennedy and O'Hagan but suggested to decompose the model output and model discrepancy using singular value decomposition (SVD). The weights generated in the SVD of the numerical output is treated as a function of the model parameters and spatial coordinate. The concept of the STD-BC is then applied to the weight functions. Zhang and Fu [16,17] improved the efficiency of the original Bayesian UQ framework by using the adaptive HDMR technique to decompose the original high dimensional problems into several lower-dimensional subproblems. These treatments subsequentially improves efficiency, however, it does not invoke the online learning capability. Once more benchmark data becomes available, users have to execute, from scratch, the whole procedure with a larger dataset. The work presented in this paper aims to provide a solution to improve the efficiency of the STD-BC and enable the online capability using methods developed within the research literature on Gaussian processes (GP).
GP is a naïve Bayes supervised machine learning algorithm. By assuming the data follows an adjoint multivariate Gaussian distribution, it offers the potential to determine the mean and variances of the unknown. As is the case with the STD-BC, the poor scalability of inversion of the covariance matrix has limited the usage of the GP. In recent years, many efforts had been made to develop a sparse Gaussian processes for very large data sets. A detailed review can be found in [28]. In general, there are several ways to increase the sparsity of a GP, such as the Nystrӧm GP [29,30], deterministic training conditional (DTC) GP [31], and kernel interpolation for scalable structured GP (KISS-GP) [32]. In this paper, we propose to adapt the stochastic variational inference GP (SVI-GP) method [33] to improve the efficiency of the Bayesian calibration and enable the online learning capability. We name the new method 'fixed inducing points online Bayesian calibration' (FIPO-BC).
The rest of this paper is organised as follows: A brief introduction to the SVI-GP and STD-BC is given in section 2 which is followed by the description of FIPO-BC in section 3. In section 4, we summarise the results of a toy case, while section 5 is contributed to calibrate a key model parameter in a scale resolving turbulence model (SAS-SST model). This paper is closed by a summary and discussion of future work in Section 5.

Background
According to Kennedy and O'Hagan [5], the physical process (•), noisy observation ∈ ℝ, and model outputs ℎ (•,•) at ∈ ℝ has the following form: where is the observation error and ~ (0, −1 ), δ(•) is the numerical model inadequacy and ∈ ℝ is the vector of model parameters. The Bayesian calibration assumes the numerical dataset and observations is a joint GP and then proceeds to find the most plausible via maximise the marginal likelihood.

Stochastic variational inference Gaussian processes (SVI-GP)
Before further description of the Bayesian calibration method, a concise description of Gaussian processes [34,35] and stochastic variational inference Gaussian processes (SVI-GP) [33] are given in this section. This allows us to introduce notation and expressions which will be used in the description of STD-BC and FIPO-BC methods.
The standard GP assumes that any finite subset of { ( )| ∈ ℝ } follows an adjoint Gaussian distribution. Considering the vector of noisy observations = { ∈ ℝ} =1 1 taken at In the standard GP, we have Using Eq.(1 & 2), we define the prior of the real-valued function ( ) as an adjoint multivariant Gaussian distribution with zero means ( ) and covariance matrix ( ). is an 1 × 1 matrix, whose elements are calculated using a predefined covariance function which is typically dependent upon a set of hyperparameters. Hyperparameters can be selected by maximise the marginal likelihood ( | ) = ( | , covariance vector and * * is calculated at * using the covariance function.
The computational complexity of GP is governed by the size of . Assuming, we can find inducing variables that are conditioned on inducing points = { } =1 1 , preferably with 1 ≪ 1 . We then augment the GP prior as Similar to ( | ), , where is a 1 × 1 covariance matrix evaluated between all the inducing points. Let (a 1 × 1 matrix) be the covariance matrix between and . We will then have, , where ̃= − −1 and is the transpose of . The task then becomes finding the proper with knowledge of . In the SVI-GP, the variational distribution ( )~( , ) is introduced to enable the usage of SVI which allows for online training [36].
There is difficulty in determining the optimised inducing inputs ( ) as it will introduce a dependence on the training data if the optimisation procedure is done w.r.t it. To solve this problem, a fine predifined inducing points can be used [33,37].
Parameters in ( ), including , , and hyperparameters of the kernel function, are optimised using the lower boundary of ( | ).The objection function can be written as with being the i th row of and the ̃ is the diagonal element of ̃. Once ( ) is determined, the posterior GP mean and covariance can be approximated as described in [38].

Standard Bayesian calibration (STD-BC) revisited
Consider observations = { ∈ ℝ} =1 1 at = { , ∈ ℝ } =1 1 , and simulation outputs The number of elements in is 1 = 1 1 . It is well known that numerical outputs are also subject to numerical error ( ). Therefore, we assume ~( , ), in which = + −1 and is a 1 × 1 matrix where the element is evaluated using In Eq. (7), is is a 1 × 1 cross-covariance matrix and is designed to describe the relationship between numerical outputs and observations with .
It should be noted that we eliminated = [ ] in Eq. (7). It is noted in [5] that marginalised the effect of and would be too complex. Therefore, the optimal hyperparameters * and * are suggested to be used. It should be noted that and are folded in * and * here.
The flow chart for STD-BC is illustrated in the Figure 1. STD-BC contains three major steps: 1. Determine a set of optimal hyperparameters ( * ) for w.r.t. the marginal likelihood of ~( , ).
2. Determine a set of optimal hyperparameters ( ) for with knowledge of * .

Solve
, and In Eq. (12), = and = . and = + ′ are the covariance matrix for ̃ and ̃, respectively. ′ is the covariance matrix with element calculated using kernel function for the model inadequacy related to ̃. is the cross-covariance matrix for ̃ and ̃ with element evaluated using and = .
We further write the variational distribution of ̃ and ̃ to be Expectation of ( | ,̃ * ,̃ * ) given ( ) is By eliminating the hyperparameters in the expression, we can now comfortably write 〈 ( | ,̃ * ,̃ * )〉 ( ) as ( | ) and its logarithm form is   [39]. The validity of the GP/SVI-GP models for ℎ( , ) and measurements ( ) are checked using the validation dataset. For the sake of simplicity, the GP model is expressed as and the SVI-GP model is expressed as . Several attempts to determine the optimal hyperparameters have been done.
In each attempt, the initial searching point of the adaptive Nelder-Mead simplex algorithm is carefully selected. Only the best-performed model is considered in the following discussion.
The quantity-to-quantity (Q-Q) plot and the mean square error (MSE) are illustrated in Figure   4. Outputs of ( , ) are very similar to those of ℎ( , ) except for Case 1, seeing Figure   4      Major efforts have been spent in the development of turbulence models which describe the effect of turbulent eddies sufficiently well. In this section, the FIPO-BC algorithm is used to calibrate a damping factor ( ) in a scale-resolving CFD turbulence model (SAS-SST). The SAS-SST model [40,41] aims to generate turbulent eddies using information obtained from the mean flow. Unrealistic turbulent eddies in the high-wave number region may be introduced into the simulated turbulent flow system. As a result, a specific term is designed to filter those unphysical turbulent eddies. The damping factor ( ) controls the size of the filter and is the key parameter in this term. It changes the balance of the modelled and resolved turbulence effect, and hence affecting the model accuracy.
The FIPO-BC algorithm is applied to calibrate the with regard to the streamwise velocity measurement of incompressible flow passing a prism bluff body [42,43]. The SAS-SST model is implemented within the commercial CFD solver, STAR-CCM+ 12.04, using user- It suggests the SAS-SST using between 0.9 and 1.0 would produce almost the same / at the location.
The FIPO-BC with several different sets of inducing points are performed of which details are included in Table 2. In the first three cases, the number of the inducing points for gradually increases from 6 to 11, while it is kept the same (36) on the spatial coordinate ( ).
Spatial inducing points increase to 41 and 46 in the last two cases, respectively. Again, the adaptive NM simplex algorithm is chosen to obtain the hyperparameters. But in this part of the study, the initial searching points are randomised. Dozens of attempts have been made and the hyperparameters resulting in the best results are chosen to show here.
Suitability checking of the surrogate models for CFD models and measurements are illustrated in Figure 9. As in the toy case, the surrogate models in the FIPO-BCs are checked by using the inducing data points to predicted the training dataset, whilst the LOO-CV test is adopted to assess the validity of surrogate models in the STD-BC. As shown in Figure 9  the SVI grid converged SVI-GP models, as its MSE is ~2.0E-02. This may be due to the sparsely distributed noisy data points and the LOO-CV test. Figure 10 illustrates the histogram of the posterior PDF of (obtained using the MH algorithm), as well as profiles of log marginal likelihood and maximums in FIPO-BC cases and STD-BC are plotted. As the inducing points refined, the shape of ( ( | )) becomes more like that of ( ( | )). Since the effect of hyperparameters is not marginalised, the range of the marginal likelihood greatly varies and causes differences in the posterior PDF. For instance, the range of ( ( | )) in FIPO-BC-2 and 4 are much larger than that in the other cases. Subsequently, the posterior distribution of in these three cases is more concentrated than others. To demonstrate the online FIPO-BC, the original data is split into four different data batches and fed into the FIPO-BC solver sequentially, whilst the inducing is the same as that used in FIPO-BC-5. The results of this online calibration are shown in Figure 12. Each pane of Figure 12 includes the used to generate the numerical outputs, response surface of absolute error (|numerical output -SVGP prediction|), as well as the optimal . Clearly, with more data fed into the FIPO-BC solver, the absolute error is reduced across the parameter space. Meanwhile, the optimal continuously improves. The optimal obtained using 'Batch-1, 2, 3 & 4' are 0.5, 1.0, 1.0, and 0.94, respectively. In particular, the absolute errors produced at the stage of 'Batch-3' and 'Batch-4' are very similar. We can confidently claim that the calibration can be achieved only using Batch-1, 2 & 3 when the FIPO-BC is used.      It has been demonstrated in both test cases that FIPO-BC algorithm provides the same calibrated parameter as the STD-BC algorithm, once the FIPO-BC algorithm is converged on the number of inducing points. Due to the online learning capability of the FIPO-BC algorithm, the calibration can be achieved using less data. Like the STD-BC algorithm, the posterior distribution of the parameter will be affected by the optimal hyperparameters. Both the STD-BC and FIPO-BC algorithms require marginalisation to eliminate this effect. In order to enable the algorithms to utilise marginalisation requires additional computational resources. However, the inclusion of marginalisation within the FIPO-BC algorithm would be more computationally efficient than the STD-BC algorithm with marginalisation. The FIPO-BC algorithm with marginalisation will be investigated in future work.