{'fraction_non_alphanumeric': 0.04589914505144182, 'fraction_numerical': 0.012389508766845385, 'mean_word_length': 4.505584363781412, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-fidelity Bayesian experimental design to quantify extreme-event statistics

In this work, we develop a multi-fidelity Bayesian experimental design framework to efficiently quantify the extreme-event statistics of an input-to-response (ItR) system with given input probability and expensive function evaluations. The key idea here is to leverage low-fidelity samples whose responses can be computed with a cost of a certain fraction of that for high-fidelity samples, in an optimized configuration to reduce the total computational cost. To accomplish this goal, we employ a multi-fidelity Gaussian process as the surrogate model of the ItR function, and develop a new acquisition based on which the optimized next sample can be selected in terms of its location in the sample space and the fidelity level. In addition, we develop an inexpensive analytical evaluation of the acquisition and its derivative, avoiding numerical integrations that are prohibitive for high-dimensional problems. The new method is tested in a bi-fidelity context for a series of synthetic problems with varying dimensions, low-fidelity model accuracy and computational costs. Comparing with the single-fidelity method and the bi-fidelity method with a pre-defined fidelity hierarchy, our method consistently shows the best (or among the best) performance for all the test cases. Finally, we demonstrate the superiority of our method in solving an engineering problem of estimating the extreme ship motion statistics in irregular waves, using computational fluid dynamics (CFD) with two different grid resolutions as the high and low fidelity models.

network [21,22] as surrogate models, multi-fidelity sampling algorithms have been developed for the purpose of global optimization [23,24], function learning [22], and contour detection [25]. In terms of the extreme-event statistics, the only related work within the multi-fidelity framework, to our knowledge, is [26] which estimates the exceeding probability for reliability analysis. However, [26] employs the sub-optimal acquisition function developed in AK-MCS which has been shown to be not only (much) less efficient than many improved algorithms later, but also not applicable to our purpose of obtaining the overall extreme-value portion of the response PDF.
One of the key issues involved in the multi-fidelity sampling method is the determination of the fidelity level of each sample. Two types of methods have been considered regarding this issue, one to select the fidelity of next sample adaptively based on existing samples in order to reduce the overall computational cost (though in a heuristic manner) [25,17], and the other to follow a pre-defined fidelity hierarchy (say pre-defined ratio and sequence of high and low fidelity samples in a bi-fidelity context) [24,18]. While the former type of method is developed in the hope of outperforming the latter one, there are not sufficient evidences to support the idea as systematic comparisons between the two types of methods are not available. It is one purpose of the current study, in the context of capturing the extreme-value response PDF, to systematically compare these two types of methods in determining the fidelity level of samples, along with the identification of their improvements relative to the single-fidelity algorithm [27].
Specifically, in the present work, we develop a multi-fidelity sequential BED framework for the quantification of the response PDF of an ItR system, with emphasis on the extreme-value portion. In particular, we use a multi-fidelity Gaussian process as the surrogate model, and develop an acquisition function (as a substantial extension to the single-fidelity function [27]) which allows adaptive choice of both the location in the parameter space and fidelity level of the next sample. We also construct an analytical computation of the acquisition which avoids expensive numerical integration and enables high-dimensional implementation of the algorithm through gradient-based optimization. Our new method is systematically tested in a bi-fidelity context for a series of synthetic problems with varying dimensions, low-fidelity model accuracy and computational costs. We show that our bi-fidelity method outperforms the single-fidelity method in all test cases, and that our method for adaptive choice of fidelity level consistently performs among the best in all bi-fidelity runs with pre-defined fidelity hierarchy varying in a broad range. Finally, we demonstrate the coupling of our method with CFD to compute the PDF of extreme ship roll motion in irregular ocean waves. By using CFD simulations with two different grid resolutions as high and low fidelity models, we show that our bi-fidelity method achieves much faster convergence of the result (i.e., extreme response PDF) than the previous single-fidelity method.
The python code for the algorithm, named MFGPextreme, is available on Github 1 .

Problem setup
We consider a black-box ItR function f (x) : R d → R, with x ∈ R d representing the input parameters with known probability distribution p x (x). We assume that we have a hierarchy of models f 1∼s (x) (from low to high fidelity) to compute f (x), with f s (x) = f (x) and f i (x) having increasing deviations from f (x) for i = s − 1, s − 2, ..., 1. In addition, the models f i (x) are associated with fixed computational costs c i which increases for i = 1, 2, ..., s. The evaluation of these models are generally corrupted by noise, with observation y i defined as where i is the Gaussian noise with variance γ i (representing observation error).
Our quantity of interest is the PDF of the response p f (f ), focusing on the tail part. Specifically, we aim to obtain an estimation p f,est (f ) with minimized (see [12]) with p f,true (f ) the true PDF of the response. We note that the log function in (2) acts on the ratio p f,est (f )/p f,true (f ), which is amplified when p f,true (f ) is small, i.e., (2) emphasizes on the error in the small-probability portion (in many cases extreme-value portion) of the PDF.
To compute p f,est (f ), we can use a sequence of samples y i (x) with i and x varying for each sample. Our objective is to find an optimized sequence, in terms of both i and x, such that e is minimized under a given total computational cost c (i.e., summation of computational cost c i over all members in the sequence). In general, there is no solutions to this type of problem that can be guaranteed to be global optimal, and the method we propose in this paper should be considered as a greedy algorithm that looks one step ahead of the existing samples. In particular, 1 https://github.com/umbrellagong/MFGPextreme our method is based on a multi-fidelity sequential BED, which involves two basic components: (1) an inexpensive surrogate model based on the multi-fidelity Gaussian process which infuses the information of multi-fidelity samples; (2) a new acquisitive function measuring the benefit (i.e., reduction in e) per computational cost, through the optimization of which the next-best sample can be selected in terms of both i and x. The two components are next described in detail in §2.2 and §2.3. In addition, in §2.4, we develop an analytical formula to compute the acquisition function and its derivative with respect to x, enabling the gradient-based optimization that is suitable for high-dimensional problems.

Surrogate model
In this section, we briefly outline the multi-fidelity Gaussian process developed in [20] as our surrogate model. Assume we have a dataset D = {X , Y} consisting of s levels of model outputs sorted by increasing fidelity. The purpose of the multi-fidelity Gaussian process is to learn the underlying relation f i (x) from D. This can be achieved through an auto-regressive scheme, which models f i (x) as with f 1 (x) ∼ GP(0, k 1 (x, x )) and {d i (x) ∼ GP(0, k i (x, x ))} i=s i=2 pairwise independent Gaussian processes, ρ i−1 a scaling factor to quantity the correlation between f i and f i−1 . The kernels k i (x, x ) are defined as radial-basis functions with τ i and the diagonal matrix Λ i respectively representing the characteristic amplitude and length scales.
are hyperparameters in the model and can be determined by maximizing the likelihood function p(Y = Y), where Y is a random vector of the (noise corrupted) surrogate with input at X , satisfying a Gaussian distribution N (0, cov(Y)).
Here we apply the shorthand notation cov(Y) = cov(Y, Y) to represent a covariance matrix for each pairwise random variables in Y, which will be used throughout this paper.
The posterior prediction f i (x) given the dataset D can then be derived as a Gaussian process with analytically tractable mean and covariance (also defined across different fidelity levels) In (6) and (7) (as well as the likelihood function), the covariances are computed as (or can be derived from): where We finally summarize the bi-fidelity counterpart of (5) in Appendix A, which we will use in §3 for computation.

Acquisition for sample selection
Given the Gaussian process surrogate f s (x)|D as in (5) of the ItR, we can estimate the response PDF p f |D (f ). Our purpose is to select the next sample in terms of the fidelity level i and the locationx to significantly reduce the uncertainty in the extreme-value part of the response PDF (which is expected to lead to significantly smaller e in (2)). In particular, this uncertainty can be estimated by (see previous work in the single-fidelity context [12]) where, y i = E(f i (x)|D) is the mean response computed by the surrogate f i (x)|D from a hypothetical locationx and fidelity i, p f ± |D,y i (x) (f ) are PDF bounds generated by upper and lower bounds (say two standard deviations away from the mean) of f|D, y i (x).
Using U (D, i,x) directly as the acquisition, however, involves significant computational cost (e.g., building a new Gaussian process f|D, y i (x) for each hypothetical sample) even for singlefidelity problems. To address this issue, we extend the methodology developed for single-fidelity problems in [14,27] to the multi-fidelity context. The first step is to introduce an upper bound as a proxy to U , defined as (proven in [14] for single-fidelity applications) with where f (x) = E[f (x)|D] represents the mean prediction (6) with i = s. Q measures the model uncertainty with emphasis on positions of large w, i.e. rare (usually large) response regions with significant input probability. We are interested in the reduction in Q (i.e., the benefit) after adding the sample atx and i, formulated as The derivation of the result in (13) is summarized in Appendix B, which makes use of the recursive update of the Gaussian process that is simpler than the derivation in [27] for the single-fidelity method. We note that the expensive computations of the new posterior in (10) and (11) are not involved in (13). Further reduction of computational cost by avoiding the numerical integration in (13) will be discussed shortly in §2. 4.
In general, one may expect that adding a high-fidelity sample is more beneficial than adding a low-fidelity sample at the same x. While this is indeed generally true, we note that there exist some certain special situations in which adding a low-fidelity sample becomes more beneficial according to (13). Such situations can occur when the function d i (x) in (3) becomes uncorrelated, with a rigorous justification provided in Appendix C. To select the next best sample in terms of both location and fidelity level, we need an acquisition function taking into consideration both the benefit (13) and cost of the sample c i . In particular, we solve an optimization problem We remark that (14) provides the optimal next sample in terms of the uncertainty reduction (13) per computational cost. Nevertheless, there is no guarantee that successively applying (14) provides a global optimal solution although this type of fidelity-choice algorithm has been also applied for other purposes [25,17]. Therefore, the ultimate validity of (14) needs to be tested in a sufficiently wide range of examples, especially against algorithms with a fixed fidelity hierarchy, which is one of our purposes in §3.
In solving (14) as a combined discrete and continuous optimization problem, we first find the optimal location x for each fidelity i, i.e., . . , s, then we compare the benefit per cost for each fidelity level and find the optimal fidelity level, i.e., While the solution procedure outlined above seems straightforward, there still exists a difficulty for applying the method to high-dimensional problems. The reason lies in that the numerical integration in (13) can become prohibitively expensive for high-dimensional x. Furthermore, the high-dimensional optimization (14) needs to rely on gradient-based algorithm where the derivative of (13) is also expensive to compute. To address these issues, analytical formulae for (13) and its derivative are much preferable, which will be discussed in the next section.

Analytical formulae for (13) and its derivative
To develop an analytical formula for (13), we first substitute the expression of covariance function (7) into (13) and obtain with We see that every term in (15) is analytically tractable except the K function in (16) where the numerical integration is carried on. One idea to obtain an analytical form of K, which has been suggested in the single-fidelity cases [27,28], is to approximate the w(x) with a Gaussian mixture model [29] with n GM M Gaussian functions: This allows us to reformulate (16) as with The problem now boils down to developing an analytical formula for (19), which involves in the integrand the multiplication of two different multi-fidelity covariance functions and a Gaussian distribution function. This situation here is more complicated than that in the single-fidelity case [27] where the problem is simplified by only involving two same single-fidelity covariance functions (e.g., analytical result from the latter case is already available [30]). We summarize the detailed derivation of the analytical form of (19), as well as the derivative ∂B(i,x)/∂x that can be derived in a similar manner, in Appendix D.
With analytical computation of B(i,x) and ∂B(i,x)/∂x available, we can solve the optimization (14) using gradient-based algorithm (which is more suitable for high-dimensional problems). In our current work, a gradient-based quasi-Newton method [31] with multiple starting points is used to solve (14), which completes the algorithm of the multi-fidelity BED method for extreme-event statistics. We finally summarize the full algorithm in Algorithm 1 and note that the algorithm reduces to a single-fidelity BED method for s = 1.

Algorithm 1 Multi-fidelity BED for extreme-event statistics
Require: (12) and approximate it with GMM model (17) 3. Solve the optimization (14) to find the next-best sample {i * , x * } 4. Evaluate the i * −fidelity function to get Output: Compute the response PDF based on the surrogate model (6) In this section, we test our developed method in the context of bi-fidelity problems, i.e.,  (2) is obtained from a computation using a sufficiently large number of high-fidelity samples (since we assume f (x) = f h (x) discussed in §2.1). In all cases, we compare the results from our method (bi-fidelity optimal sampling in both location and fidelity level as in Algorithm 1, hereafter termed "BF-O") to those from single-fidelity optimal sampling [27] (Algorithm 1 with s = 1, hereafter termed "SF"). In addition, for the synthetic problems with d = 2 and 8, we further include the results from the bi-fidelity model with fixed ratio n (as well as sequence) of low and high fidelity samples with locations optimized (Algorithm 1 but with solved for fixed i on each sample, hereafter termed "BF-Fn"). In these cases, we also vary the number n, the ratio c h /c l and the low-fidelity accuracy level in a broad range, to assess the performance of our method BF-O in various situations for cases of intermediate to relatively high complexity.
For all the methods in comparison, an initial data set is required to start the sequential BED procedure. To create a fair situation for comparison, we keep the cost to generate the initial dataset the same for all methods. In particular, the initial set is generated by a space-filling Latin Hypercube sampling method [32], with 4d high-fidelity samples for SF method, and 2d high-fidelity and 2d c h /c l low-fidelity samples for BF-Fn and BF-O methods. Since the initial sample locations are not fixed, we will present the results from SF, BF-O and BF-Fn methods in terms of the average over 100 implementations with different initial datasets for the synthetic cases in §3.1, §3.2, and §3.3.

One-dimensional Forrester function
We start the method validation from a one-dimensional (1D) Forrester function f (x) that has been previously used to demonstrate the multi-fidelity global optimization [33]. The high- where the input x is assumed to follow a Gaussian distribution with p sample as the first sequential sample. This is consistent with an intuitive understanding that the algorithm tends to use three high-fidelity samples (2 in the initial dataset and 1 as selected) to learn the linear difference term 10x in (21) using a Gaussian process, which results in a significant reduction in e at c ≈ 5 in figure 1(b).
While this simple example demonstrates the advantage of using BF-O method, it is only for a single case with fixed accuracy level of the low-fidelity model (in terms of (21)) and the cost ratio c h /c l . We will next use a two-dimensional (2D) case to test a much broader range of situations in the next section.

Two-dimensional stochastic oscillator
We consider a 2D function constructed from the solution of a stochastic oscillator equation, which has been previously used for testing the single-fidelity BED method (i.e., SF method) in [12,27]. In particular, the oscillator equation is formulated as where u(t) is the state variable, F is a nonlinear restoring force defined by: The stochastic process ξ(t), with a correlation function σ 2 ξ e −τ 2 /(2l 2 ξ ) , is approximated by a twoterm Karhunen-Loeve expansion with λ i and φ(t) respectively the eigenvalue and eigenfunction of the correlation function, x ≡ The response of the system is considered as the mean value of u(t; x) in the interval [0, 25], which serves as our high-fidelity model: For our low-fidelity model, we construct a function f l (x) (to be varied later in this section) with a difference d(x) from (24):  To further test the performance of BF-O (as well as BF-Fn) in more diversified situations, in the following, we construct additionally 9 cases with different computational costs c l and lowfidelity functions f l (x) as summarized in Table 1. In addition to the linear difference function used in the previous case (now case 1 in Table 1), we also consider a nonlinear difference function d(x) = δ sin(x 1 + x 2 ) which is known to be a difficult situation to approximate by a Gaussian process [6]. i.e., cases 2, 3, 5, 8-10 in Table 1 with varying ρ and δ, with the range of δ corresponding to difference terms up to 20% of the maximum response f h (x) for |x| < 4. Together with figure   2, we see that the BF-O scheme consistently performs among the best of all tested methods.
More specifically, with the increase of complexity in the difference function (i.e., increasing δ but not much for increasing ρ), the performance of BF-O can deteriorate for smaller number of sequential samples (e.g., figure 3(f)), but still behaves close to the optimal at larger number of samples with cost c ≈ 80. These cases correspond to the situation where the difference terms are difficult to learn, which takes more cost to make low-fidelity samples useful for the final results.  Table 1. Similar to the results above, we see that the BF-O method consistently performs among the best. This indicates that the acquisition function (13) employing the ratio between benefit B(i, x) (uncertainty reduction) and cost c i effectively captures the optimal that balances the two factors.     we include in the table the average ratio of low and high fidelity sample numbers n l /n h in the BF-O method, as well as the value of n corresponding to the BF-Fn method with the best performance at c = 80. While the optimal n in the BF-Fn method does not necessarily correspond to the n l /n h in BF-O method, we find that the two numbers are close in most cases, with the latter automatically captured by the algorithm.

High-dimensional borehole hydrological model
We next consider an eight-dimensional borehole hydrological model, which has been used as an example of high-dimensional problems to quantify the extreme response statistics using single-fidelity methods [34,15,35]. The model physically computes the flow rate through a borehole, formulated as with an eight-dimensional input x = {r w , r, T u , H u , T l , H l , L, K w }, including their distributions, detailed in Table 2. We further construct a low-fidelity model as where we use a nonlinear (sinusoidal) form of the difference function and put more weights on the parameter x 1 = r w , as it is the most influential factor to the response f h (x) [34]. In particular, the coefficient 7.5 is chosen such that the term of sin(x 1 ) corresponds to 5% of the maximum response. We keep the computational cost ratio as c h /c l = 5 for this case.
The results from the SF, BF-O and BF-Fn methods are shown in figure 5. Comparing to the SF result, it is clear that the benefit of using bi-fidelity models is more evident than the low-dimensional cases, even for the initial data set without sequential samples. The reason may lie in the fact that the "value" of a high-fidelity sample becomes compromised with the increase of dimensions. In addition, the BF-O method again performs the best among all BF-Fn methods with varying n.
We finally remark that the computation in this eight-dimensional case is only enabled because of (i) the development of (13) (13) with numerical integration ( ) and analytical formula ( ) as objective functions for different sizes of the existing dataset. In the former two computations, we assume that 10 times (a conservative number) of acquisition evaluations are needed compared to the gradient-based optimization in the third case, as a common practice found in [36]. more times in gradient-free optimization than that in the gradient-based method, resulting in prohibitive computational costs for large number of samples. In contrast, using the analytical formula combined with the GMM model, the computation takes only O(100)s even for a dataset of 1000 sample points, which is supposed to be negligible compared to the evaluation of the output of the high-fidelity physical model.

Coupling to CFD to compute extreme ship motion statistics in waves
We further consider an application of our method to evaluate the PDF of extreme ship roll motion in irregular waves. More specifically, we study the motion of a two-dimensional, squareshaped hull geometry with 40m × 40m cross section and density ρ h = 0.5ρ w (with ρ w the water density) subject to beam waves. The input to this problem is considered as x = {A, L}, with A and L the wave group amplitude and length, as a reduced-order description of an uni-directional irregular wave field [37,13]. Figure 7 shows an example to evaluate parameters A and L from a given wave field as well as the resulted probability distribution p x (A, L). The irregular wave field is described by a Gaussian spectrum in the form with the significant wave height H s = 12m, peak (carrier) wavenumber k 0 = 0.018m −1 (corresponding to peak period T p = 15s), and K = 0.05k 0 . The response in this case is considered as  are shown in figure 9 for two examples with different input wave parameters, showing that the difference of the results from the two models (in general) increases with the group amplitude A.
We use as the initial dataset 4 high-fidelity and 16 low-fidelity samples in the BF-O method, and 9 high-fidelity samples in the SF method (resulting in almost the same cost). Since the exact PDF in this case is unknown, we directly use the uncertainty level of the extreme-value PDF, i.e., U as in (10), as an evaluation of the quality of results, which is plotted in figure 10(a) as a function of the computational cost measured by computation time in hours (only CFD times).
It can be seen that the BF-O method results in a faster convergence of the extreme-value PDF compared to the SF method. This point is further illustrated in figure 10(b), which plots the PDFs from SF and BF-O method for the same computation time of 12 hours together with the upper and lower bounds (in terms of two standard deviations two-sided from the mean). It is clear that the PDF from the BF-O method is associated with appreciably lower uncertainty.

Conclusion
In this paper, we develop a multi-fidelity sequential Bayesian experimental design framework for efficient evaluation of the extreme response PDF. Our method leverages the multi-fidelity Gaussian process as a surrogate model, and a new acquisition function which allows the selection of the next-best sample in terms of both the location and fidelity level. We also construct an analytical formula for the acquisition function, which enables implementation of the method (e.g., gradient-based optimization) for high-dimensional problems. Our new method is first tested in a bi-fidelity context for a series of synthetic problems. With a broad range of low-fidelity accuracy level and computational cost, we show that the bi-fidelity method always outperforms the single-fidelity method, and that the BF-O method consistently shows advantage over the BF-Fn method, i.e., the bi-fidelity method with pre-defined fidelity hierarchy. We finally demonstrate the effectiveness of our BF-O method (relative to the single-fidelity method) in an engineering problem to evaluate the extreme ship roll PDF in irregular waves, where CFD of two resolutions serve as the high and low fidelity models. Further improvements on this framework may be achieved by considering reinforcement learning, which will be a topic of future study.

A Bi-fidelity Gaussian process
The bi-fidelity Gaussian process conditioned on D = {{X 1 , X 2 }, {Y 1 , Y 2 }} can be considered as the counterpart of (5) for s = 2, with its mean and covariance formulated as cov where cov( cov( cov B Derivation of (13) using recursive update For the derivation of (13), we consider the following Bayes' theorem: where f (x)|D, y i (x) can be seen as the posterior of f (x)|D (as a prior) after adding one sample y i (x). One can then get the mean and variance of f (x)|D, y i (x) using the standard conditional Gaussian formula: The formula in (13) is a direct result of (36).

C Example of a more informative low-fidelity sample
We show a special case in Theorem 1 where a high-fidelity sample is less informative than a low-fidelity sample in a bi-fidelity model (s = 2).
With the analytical formula for G t (f i (x 1 ), f j (x 2 )) available in (42), the benefit B(i,x) in (15) can be derived accordingly. To further obtain the analytical derivative of B(i,x), we rewrite B(i,x) = T (i,x)/var(y i (x)|D) (comparing to (13)), and derive the derivative of T as The analytical computation of (49) requires the formula for the derivative of the covariance function and K. For the former, we have: with ∂k(x 1 , x 2 ; Λ) ∂x 1 = k(x 1 , x 2 ; Λ)Λ −1 (x 2 − x 1 ).