Rate-Informed Discovery via Bayesian Adaptive Multifidelity Sampling

Authors: aman  sinha
payam  nikdel
supratik  paul
shimon whiteson waymo
Venue: robot-learning_org_CoRL_2024_Conference
Published: 2024-06-07 00:33:42

Abstract: Ensuring the safety of autonomous vehicles (AVs) requires both accurate estimation of their performance and efficient discovery of potential failure cases. This paper introduces Bayesian adaptive multifidelity sampling (BAMS), which leverages the power of adaptive Bayesian sampling to achieve efficient discovery while simultaneously estimating the rate of adverse events. BAMS prioritizes exploration of regions with potentially low performance, leading to the identification of novel and critical scenarios that traditional methods might miss. Using real-world AV data we demonstrate that BAMS discovers 10 times as many issues as Monte Carlo (MC) and importance sampling (IS) baselines, while at the same time generating rate estimates with variances 15 and 6 times narrower than MC and IS baselines respectively.

Main paper: 
 introduction: 
 Commercial autonomous vehicle (AV) development usually follows an iterative process wherein issues with the current planner are identified, root causes are determined, and the issues are addressed for the next release. While these software updates are meant to improve performance, they can also lead to regressions. Evaluating the long tail of safety issues is particularly difficult: as the AV's performance improves, safety-related issues become rarer and more difficult to discover. In turn, this makes improving the planner increasingly challenging, as identifying failure cases is the first step towards addressing them. Together with discovering particular failure modes, it is also important to efficiently estimate the overall safety or failure rate of any potential release before deployment.Improvements to an AV planner can be tested on the road with human drivers in the vehicle ready to take over in the event of imminent danger. However, this approach is expensive and does not ensure sufficient coverage of driving scenarios since we cannot control the environment [1] . Furthermore, since this approach essentially produces a Monte Carlo estimate of the AV's performance, the scale required for evaluation near human-level performance is too large [2] to repeat at regular intervals during development. As such, simulation plays a key role in industrial AV software evaluation (see, e.g., [3, 4] ), as it enables testing software in diverse settings, at a fraction of the cost and time, and with no risk to the public. However, a naive simulation-based approach to evaluation also suffers from intractable scaling and can be prohibitively expensive.In this paper, we address the twin problems of (i) estimating the safety of a given AV planner and (ii) discovering failure cases that can then be used to drive further improvements. Concretely, given a distribution P of simulation parameters that describe the AV's operational design domain, the governing problem for evaluating safety is to estimate the rate of adverse events:EQUATIONwhere X ∼ P , the function f : X → R is a performance metric scoring the realization x ∈ X of the AV and its environment, and γ is a user-defined threshold below which behavior is undesirable. Estimating p γ is a rare-event simulation problem, so naive Monte Carlo methods are expensive and we must resort to more tractable IS techniques; these techniques place high probability on regions of X where f (x) is small. However, simply estimating p γ (the objective of importance sampling (IS)) is insufficient, as simulation testing must provide actionable feedback on how to improve AV design (Figure 1 ). This feedback comes in the form of concrete examples x i where f (x i ) ≤ γ. In particular, during the continuous Figure 1 : Rate-informed discovery loop -A: A diverse set of simulations are run with the latest version of the AV planner. B: Some simulations lead to adverse events which are attributable to the planner. C: The results from all simulations are used to compute an unbiased estimate of the rate of adverse events. D: The adverse events found in B are used to understand and address the failure cases. E: An improved planner is deployed, and the cycle repeats. feedback loop of simulation testing followed by design improvement, it is most useful for simulation to discover failure cases x i with a) relatively high likelihood under P and b) novelty with respect to previously discovered failure examples. From a practical perspective, prioritization by likelihood is obvious when trying to allocate resources to fixing problems, and novelty is important because previously-discovered failure modes typically have existing work underway to improve behavior.We define rate-informed discovery as the problem of returning novel, likely examples of undesirable AV behavior in service of confidently estimating the rate p γ . Our goal is to perform rate-informed discovery as efficiently as possible. Just as with purely estimating the rate p γ , naive Monte Carlo is intractably expensive for rate-informed discovery. Furthermore, IS techniques that rely heavily on previously determined failure modes are poor at discovering novel issues. In this paper, we develop Bayesian adaptive multifidelity sampling (BAMS), which adaptively samples from X in a targeted fashion and builds an importance sampler for estimating p γ . As efficiency is paramount, we exploit the availability of a simulation hierarchy: low-fidelity yet cheap simulators as well as high-fidelity but expensive simulators that model real-world environments as accurately as possible. As we show below, this multifidelity simulation setup enables BAMS to perform better rate-informed discovery for any given simulation budget than using the high-fidelity simulator alone. 
background and related work: 
 Gaussian processes A Gaussian process (GP) [5] imposes a prior belief on an unknown random function θ : X ⊆ R d → R through its mean function µ(x) (for convenience usually set to zero) and covariance function k(x, x ), such that any finite collection of function values {θ(x i )} has a joint Gaussian distribution. A popular choice of the covariance function is the Matérn kernel:k(x, x ) ∝ t ν K ν ( √ 2νt), where ν is a lengthscale parameter (typically ν = 2.5), K ν is the modified Bessel function of the second kind, and t := (x -x ) T Σ -1 (x -x ) is a distance scaled by a diagonal matrix Σ of lengthscales for each dimension. Letting X n ∈ R n×d be the matrix of observed inputs and θ(X n ) be the vector of associated function values, the posterior prediction θ(x * ) for an input x * is a Gaussian distribution N (µ n (x * ), σ 2 n (x * )):µ n (x * ) = k(x * , X n )k(X n , X n ) -1 θ(X n ), σ 2 n (x * ) = k(x * , x * ) -k(x * , X n )k(X n , X n ) -1 k(X n , x * ). The posterior covariance between any two inputs (x, x ) is cov n (x, x ) = k(x, x ) - k(x, X n )k(X n , X n ) -1 k(X n , x ).There is a vast literature on speeding up GP inference, namely assuaging the O(n 3 ) time complexity of inverting k(X n , X n ). In particular, inducing point approximations [6, 7, 8] impose structure on the GP prior, which in turn induces (block) sparsity onto k(X n , X n ). Whereas inducing point methods make posterior mean inference more tractable, other approaches make predictive variances cheaper by exploiting sparsity in the underlying data and employing iterative techniques for approximations [9, 10] . All of these approaches are motivated by the regime wherein the training set size n is typically much larger than the test set. Our setting is exactly the opposite; the discovery process occurs frequently in AV development cycles, so it is highly budget-constrained and n is small. In our experiments, we simulate and use only tens of samples for the GP posterior, whereas our predictive covariance is computed on a much larger set of available points. Thus, while our approach is complementary to the existing literature and can be used in conjunction with these methods, we do not benefit substantially from them.Bayesian optimization and quadrature Bayesian optimization (BO) [11, 12] and quadrature (BQ) [13, 14, 15] use GPs for targeted function evaluation to perform minimization for expensive functions or to compute expensive integrals over X . BO and BQ methods utilize an acquisition function to select points to evaluate. Acquisition functions typically balance exploration of uncertain regions in X with exploitation of previous evaluations f (x) (see, e.g., Wilson et al. [16] for an overview).There is a rich literature on using multiple fidelities for evaluation, also known as surrogate modeling (see, e.g., [17, 18, 19, 20, 21] ). Specializing to GPs, Bonilla et al. [22] introduce the concept of multi-task GP prediction, which has subsequently been employed in BO use cases [23, 24, 25] . The multi-task covariance kernel is highly structured as a tensor product of a base kernel and a covariance between discrete tasks [22] . By contrast, Marco et al. [26] and Poloczek et al. [27] consider modeling the dependence between fidelities via additive noise GPs. We employ this modeling technique in BAMS due to its generality; it requires little prior knowledge about the relationship between multiple simulators and allows us to easily incorporate costs of each fidelity into our acquisition function, especially when multiple fidelities are available for the same data point x.AV safety-measurement techniques AV evaluation literature generally considers approaches for formal verification [28, 29, 30, 31] , (probabilistic) falsification [32, 33, 34] , and rate estimation [35] . Formal verification aims to define and then prove correctness of behavior; this is both NP hard [36] and rife with subjectivity on notions of correctness. Falsification aims to find any problems with the AV system under test, and probabilistic variants aim to find high-likelihood failures [33] . Finally, rate estimation aims to measure the probability of adverse events under a (possibly empirically defined) distribution of driving environments [35] . Previous papers note a dichotomy between falsification and rate estimation [37, 38] , which our paper aims to bridge; our rate-informed discovery loop first evaluates targeted samples solving an optimization problem (similar to falsification), and we follow suit with updating an importance sampler to measure rates. Unlike many probabilistic falsification approaches, which require whitebox access to the system [34] , our formalism performs optimization with entirely blackbox access. See Appendix F for a detailed discussion on how our Bayesian approach connects to existing rare-event simulation techniques, addresses key challenges within this domain, and contributes to the broader field of AV safety evaluation.We train an encoder-decoder model to reduce the embedding to a 12-dimensional representation, which is more computationally tractable for GPs.Since the 32 AV planner rollouts in WOSAC are generated based on one second of history, our embedding model utilizes the embedding at the start of the run segment.Performance metric We utilize the 32 rollouts to calculate a single time-to-collision (TTC) value for each run segment using the formulaf (x) = Q 0.25 (mT T C [1] (x) . . . mT T C [32] (x)), where mT T C [i] (x) is the minimum TTC across the i th rollout for run segment with embedding x and Q 0.25 (. . .) is the 25 th percentile of the arguments. We use the 25 th percentile to be robust to outliers in the rollouts. Finally, we set the rate p γ at 1% by using the threshold γ = 0.43 for this metric.Experimental procedure We utilize BAMS with S = 6 clusters, based on the parameter-tuning study presented in Appendix C.1. We compare BAMS, which leverages both fidelity levels, with its ablated counterpart Bayesian adaptive sampling (BAS), which only employs the high-fidelity simulator. Additionally, we measure performance of BAMS and BAS against the following baselines:• MC: Standard Monte Carlo sampling.• MC-GP: Our method with random sampling as the acquisition function.• MCM-GP: Multifidelity version of MC-GP.• DS: IS with difficulty-model scores.• DS-GP: Our method with IS based on difficulty-model scores as the acquisition function.• CE: Cross-entropy method [47] . See Appendix E for implementation details.DS and DS-GP leverage IS based on a difficulty model. Inspired by Bronstein et al. [39] , this model predicts the intrinsic difficulty of a run segment based on previous simulations conducted under diverse planner versions. Our experimental procedure for all models is as follows:A. Initialization: All methods begin with a first batch of random samples with a budget of m 1 = 20 evaluations. This ensures unbiased exploration of the search space. Subsequently, the GP-based methods leverage this initial data to train their hyperparameters by maximizing the marginal loglikelihood (MLL) using the Adam optimizer [48] .B. Two batches with budget m b : Each method samples two additional batches, each with a budget of m b = 15 evaluations. We retrain the GP hyperparameters after each batch.Multifidelity setup For BAMS and MCM-GP, we utilize two fidelity levels:• Level 0 (high fidelity): This level uses the 25 th percentile TTC calculated from all 32 rollouts Q 0.25 (mT T C [1] . . . mT T C [32] ) with cost c(y 0 ) = 1.0. • Level 1 (low fidelity): This level uses the 25 th percentile TTC calculated from only 5 randomly selected rollouts, denoted Q 0.25 (mT T C [1] . . . mT T C [5] ), whereby c(y 1 ) = 5/32. This parameter was chosen based on a parameter-tuning study (see Appendix C.3) conducted overQ 0.25 (mT T C [1] . . . mT T C [i]) with corresponding cost i/32, evaluating the model's recall at a retention budget of 5p γ N . 
bayesian adaptive multifidelity sampling: 
 In typical AV development environments, the distribution P over which we compute the rate of adverse events (1) is an empirical distribution over run segments-specified time intervals (e.g., 30 seconds) of logged real-world driving over which simulations are performed [39, 40] . We consider the empirical distribution over N run segments, and each run segment is encoded by an embedding x ∈ R d . There are a variety of ways to create such an embedding (see, e.g., [41, 42, 43, 39, 44, 45] ); BAMS does not have any restrictions for this embedding, although smaller dimensionality d improves computational efficiency. We provide details for our embedding model in Section 4.Our formalism begins by first considering a Bayesian estimator for the rate of adverse events (1) . Specifically, we define a random function θ : X → R and impose a GP prior over it such that the collection {θ(x i )} for any subset of run segments {(x i )} has a joint Gaussian distribution. Considering a set of n evaluated run segments and their associated performance metrics X n := {(x i , f (x i ))} n i=1 , we define the shorthand θ n and E θn [•] as the posterior GP and the expectation with respect to this posterior induced by conditioning on the σ-algebra generated by X n . We can then define a function which models the presence of adverse events based on our GP posterior:EQUATIONwith marginals pγ (θ n ) := E X [g(X, θ n )] and p n (x) := E θn [g(x, θ n )]. The marginal pγ (θ n ) is an estimator of the rate of adverse events (1) .Using θ to model the performance metric f (x) of a run segment rather than the binary indicator 1{f (X) ≤ γ} has two benefits: (i) modeling the real-valued function provides better signal for the GP to interpolate, and (ii) in typical AV development paradigms, we want to understand the safety performance with respect to multiple values of the threshold γ, so we can reuse a single trained GP model. Next we describe various aspects of BAMS, motivated by computational efficiency.Acquisition function As noted in Section 2, the acquisition function defines a minimization objective that selects the next point(s) to sample the performance metric f (x) and compute the new GP posterior. Our goal of rate-informed discovery implies that choosing evaluation points to minimize the variance of the estimator, Var(p γ (θ)), is beneficial. In particular, minimizing Var(p γ (θ)) provides a confident estimate of the probability (1) . Furthermore, the calculation of Var(p γ (θ)) is dominated by regions in X that have both high likelihood under P and high uncertainty under θ n . Thus, minimizing Var(p γ (θ)) prioritizes sampling points x that are clustered near other points (high likelihood) and are far away from previously sampled points (high uncertainty)-precisely our goal for rate-informed discovery. However, this variance requires an expectation over the joint distribution X, X i.i.d ∼ P (see Appendix B.1); because P is the empirical distribution over run segments, this O(N 2 ) computation is too expensive. Instead, we consider a cheaper upper bound.First we define the point-variance function h n (x) := p n (x)(1 -p n (x)). The average of h n (x) over P , which requires only O(N ) computation, is greater than the variance of the GP estimator: Proposition 3.1. The variance of the estimator pγ (θ n ) is upper-bounded by the average point variance:Var(p γ (θ n )) ≤ E X [h n (x)].See Appendix B.2 for the proof. Since we evaluate the next m points in a batch setting, the order of these points does not matter, and our acquisition function is a function of the set {x j } n+m n+1 . In particular, consider the expectation of the point-variance function conditioning on the set of locations of these m future evaluations:EQUATIONThis forward-looking point variance has a nested expectationE θn E θn+m • |{X j = x j } n+m n+1, the average over all future posteriors conditioned on the locations {x j } n+m n+1 ; this is an analytic computation in GP models. Then, our acquisition function to choose the next m evaluation locations is the average over P of this forward-looking point variance:EQUATIONAn immediate corollary of Proposition 3.1 is that this acquisition function is an upper bound on the expected variance conditioned on the same points (see Appendix B.3 for the proof): Corollary 3.2. The expected variance of pγ (θ n ) conditioned on locations {x j } n+m n+1 of the next evaluation points is upper-bounded by the acquisition function:E θn Var(p γ (θ n+m ))|{X j = x j } n+m n+1 ≤ J {x j } n+m n+1 . Sequential selection For m, nN , the acquisition function ( 4) can be efficiently calculated in O(m 2 N ) time [46] (see Appendix A.2 for details). Because our distribution P is the empirical distribution over N samples, jointly choosing the locations for m > 1 points is a combinatorial problem over N m combinations. Since N m = O(N m /m m ), this combinatorial problem has an overall time complexity of O(m 2-m N m+1 ), which is too expensive. As a consequence, rather than choose all m locations jointly, we instead choose each one sequentially. Naively, this sequential selection takes O(m 3 N 2 ) time, but exploiting the recursive structure of the problem allows us to reduce this cost to O(mN 2 ) (see Appendix A.3 for details). In contrast, sequential selection using expected variance as an acquisition function requires O(mN 3 ) time.In the sequential setup, our minimization objective for each iteration n + 1 ≤ i ≤ n + m is: minimizex J {x j } i-1 n+1 ∪ {x} -J {x j } i-1 n+1 , where {x j } n n+1 := ∅ and J(∅) := E X [h n (x)].The subtraction of the objective from the previous iteration is not necessary for this minimization problem but it is useful to define our problem this way to motivate the corresponding multifidelity problem below.Multifidelity sampling Thus far, we have considered a model wherein there is only one way to perform evaluation on a selected point x via a single simulator. However, in our AV use case, we often have multiple simulators whose computational cost correlates with fidelity. For example, instead of simulating the entire AV planner and other traffic participants, we can simulate purely with non-interactive traffic participants or a distilled version of the AV planner. Intuitively, a noisier but cheaper simulator that correlates well with the high-fidelity simulator offers an information advantage under budget constraints: we can learn more about the ground truth with the same simulation cost. Of course, we also need to learn the correlation between simulators without incurring extra overhead. To do so, we augment our Bayesian generative model to include multiple fidelities and modify the acquisition function (4) to account for the trade-off between fidelity and computational cost.For l ∈ {0, 1, . . . , L}, let the fidelities of simulation be denoted by f l , with the original high-fidelity (l = 0) simulator denoted f 0 . Consider the augmented input y l := (x, l), which then allows us to define the augmented simulation function f (y l ) := f l (x). Next, let the model for each augmented input θ : X ×N → R be defined as θ(y l ) := θ 0 (x)+α l (x), where θ 0 is the GP model over the function f 0 , α 0 := 0, and for l ≥ 1, α l is an independent zero-mean GP with covariance function k l . Then, the augmented function θ is a GP with mean θ(y l ) = µ(x); the covariance between points y a = (x, a) andy b = (x , b) is: k(y a , y b ) = Cov(θ 0 (x) + α a (x), θ 0 (x ) + α b (x )) = k(x, x ) + δ ab k a (x, x ),where δ(a, b) is the Kronecker delta function.To update the acquisition function (4), we first define c(y l ) as the cost associated with evaluating each simulator f l . In typical AV settings, this cost is expressed in units of money or flops of computation. We define c(y 0 ) := 1 and by construction c(y l ) < 1 for l ≥ 1. Whereas in the single-fidelity setting we select a batch of m points, we now consider selecting r ≥ m points whose total cost is m. Then, the new multifidelity acquisition problem is to minimize the cost-normalized expected forward-looking point-variance function. For clarity, we write y in its expanded form (x, l):EQUATIONwhere n + 1 ≤ i ≤ n + r. Unlike in the single-fidelity case, the subtraction of the previous iteration's objective is now required-we normalize the decrease of the forward-looking point-variance with respect to the price of simulation.The theoretical results above extend directly to the mutifidelity setting where we consider P as the empirical distribution over all inputs and all fidelities. In general, we may not have the ability to evaluate every input x i at every level l, in which case we simply define our distribution P over the empirical distribution of locations and corresponding available simulation fidelities. In Section 4, we assume that all simulators are available for all N points.Clustering As noted in Section 2, our problem setting is different from most work on GPs and BO/BQ that are concerned with computational efficiency. Specifically, the number of training data points is extremely small (typically tens in our experiments) compared to N . Instead, our bottleneck is the sequential selection of points to evaluate via minimization of the acquisition function, an O(mN 2 ) computation. To decrease this cost, we employ clustering. With S clusters, the complexity of our sequential selection strategy is O(mN 2 /S 2 ). We can trivially take advantage of parallel computation in this approach as well. See Appendix A.4 for details.The complete BAMS algorithm with clustering is detailed in Algorithm 1 (Appendix A.1). Figure 2 illustrates our approach on a synthetic two-dimensional dataset. We show a schematic of Algorithm 1 in Figure 2a . Figure 2b visualizes the input distribution P and Figure 2c shows samples selected by our algorithm. For a detailed analysis of the synthetic experiment, see Appendix D. 
experiments: 
 To showcase the effectiveness of BAMS, we apply it to a realistic AV context using the Waymo Open Sim Agent Challenge (WOSAC) dataset [40] . The WOSAC dataset contains N = 44, 911 logged run segments and 32 rollouts per run segment generated via AV planner simulation. The AV planner resimulates each segment by predicting the future states of all agents, conditioned on the initial states of the logged run segment. Our goal is to estimate the rate of near-accident events. We describe all aspects of our setup, including the embedding model to generate x for each run segment, the performance metric f (x), experimental procedure, multifidelity setup, and evaluation criteria.Embedding model The input embedding is generated in two stages:1. 512-dimensional embedding: Using a dataset of run segments, we extract features such as road layouts, other road users' positions, and AV kinematics. These features are used to generate top-down images, encoded into 512-dimensional vectors using a convolutional neural network.Training uses contrastive learning on image pairs labeled as being from the same run segment or not. See Bronstein et al. [39] for details.We evaluate all approaches using two criteria: retention-recall curves and rate estimation/discovery via IS, detailed below.Retention-recall curves At the end of each batch, we use retention-recall curves to compare how effectively different methods identify novel failure modes. These curves offer insights into how efficiently each method recalls failure modes while minimizing the number of retained samples. For GP methods, the curve is generated by sorting the high-fidelity data points (l = 0) based on p n (x i ) estimates. For MC, samples are sorted randomly. Recall is then calculated at different retention levels, representing the fraction of failures found within the top-ranked samples.To facilitate interpretability, we scale the x-axis (retention) of retention-recall curves by p γ N (number of failures). Thus, a value of 1 on the x-axis corresponds to retaining the p γ N highest-scoring samples. This scaling also allows us to infer precision. For example, if recall is 100% at a retention of 2p γ N , all failure modes are captured within the top 2p γ N sample scores p n (x i ), or precision is 0.5.Retention-recall curves illustrate the trade-off between recall and the number of retained samples, which reveals a method's ability efficiently discovering failure modes. However, generation of this curve requires evaluating all x i in the high-fidelity simulator, so it is not a realistic way of measuring rate-informed discovery in AV development settings. We consider tractable statistics next. Rate estimation and discovery via IS After completing the third batch of simulation, we conduct a combined evaluation of both rate estimation and discovery. We do 200 repeated trials using model scores to perform IS with a budget of K = 5p γ N high-fidelity samples (l = 0). We calculate the resulting recall and relative variance (RV) of the rate estimate, defined as the variance divided by p 2 γ (e.g., the square of the coefficient of variation). For GP methods, we define our importance sampler as p n (x) α , with α = 2.5 chosen to minimize relative variance (see Appendix C.2). For MC, sampling is based on random scores.This IS approach guarantees that our estimate of the failure rate p γ is unbiased by construction, even if our underlying GP model θ is biased. Together, these two statistics-relative variance and recall-test each method's ability to perform rate-informed discovery: a good sampler has high recall and low variance for the rate estimate. Most importantly, generation of these statistics is a natural and tractable component of iterative AV development loops.Figure 3 presents retention-recall plots comparing BAMS, BAS, and baseline methods after each batch. BAMS achieves the highest retention at every retention level. Analyzing the results by batch reveals interesting trends. Initially, in Batch 1 methods within the same fidelity level exhibit similar performance: BAS, MC-GP, and DS-GP perform similarly, as do BAMS and MCM-GP. This is due to the use of the same random samples for initialization. However, the impact of the acquisition function becomes evident in the second batch, where BAMS significantly outperforms MCM-GP, and BAS outperforms other single-fidelity methods. This highlights the importance of the acquisition function in guiding the search towards potential failure modes.By the third batch, BAMS achieves 80% recall at a retention of 5p γ N , exceeding other methods by over 10% and highlighting its superior ability to discover failure modes. BAS achieves the best recall among single-fidelity methods, demonstrating the effectiveness of our acquisition function. Additionally, MCM-GP performs comparably to BAS, showcasing the benefit of the multifidelity approach. These results also reinforce the advantages of leveraging multiple fidelity levels for efficient discovery, as both MCM-GP and BAMS outperform their single-fidelity counterparts. Notably, CE displays worse recall than all of the GP and multifidelity methods, but, as CE is still an adaptive method, it outperforms the non-adpative MC and DS baselines.To illustrate the diversity of samples selected by each method, we use a determinantal point process (DPP) [49] to sample five diverse run segments from the collection of selected points at the end of the third batch. Figure 4 compares these run segments for BAMS and MCM-GP. The BAMS run segments include features such as busy intersections, pedestrians, and cyclists, whereas the MCM-GP run segments just show diversity in the road type (single lane vs. multi-lane). Table 1 demonstrates the evolution of focus for each method across batches, similar to the synthetic example (Figure 2c ). Notably, BAMS concentrates on low TTC areas, facilitating the discovery of diverse failure modes. While CE shows low average f (x) values in the third batch (cf. Table 1 ), the low recall performance in Figure 3 illustrates a common pitfall of CE: it collapses to sampling a subset of failure regions. Now we turn to IS statistics: Figure 5 shows rate estimates and Table 2 shows recall and relative variance (RV). BAMS achieves the best recall while exhibiting a similar rate precision to BAS. These results reinforce the superior performance of multifidelity versions of each model. Additionally, GP-based methods demonstrate better recall and rate estimation compared to DS, CE, and MC. 2c , the transition from exploration to targeted selection is most evident for BAMS.  
limitations: 
 While BAMS demonstrates substantial improvements for efficiently discovering novel, highlikelihood scenarios and estimating the rate of adverse events, our approach has limitations. Specifically, to ensure tractable GP inference, we must perform dimensionality reduction on our input embeddings; this pre-training step adds computational overhead. Furthermore, the effectiveness of BAMS hinges on the availability of a continuous (or at least real-valued) metric function f (x); this limits applicability in scenarios with purely binary labels. Finally, despite using an adaptive approach for selecting evaluation points across batches, our final IS step for rate estimation is not adaptive. Incorporating adaptive mechanisms within the IS stage, potentially leveraging techniques from sequential Monte Carlo (SMC) estimation, promises further enhancements to both the efficiency and accuracy of rate estimation. 
conclusion: 
 This paper introduced rate-informed discovery to address the critical need for efficient estimation of AV safety and discovery of high-impact failure cases. We proposed BAMS, a method that utilizes Bayesian adaptive sampling and multifidelity simulation. Our results demonstrated the effectiveness of BAMS in discovering novel, high-likelihood scenarios of undesirable AV behavior while simultaneously estimating the rate of adverse events. In future work, we plan to "close the loop" by embedding rate-informed discovery via BAMS in an iterative AV development cycle and measure its effectiveness at reducing the end-to-end cost and latency of planner improvement. 
 
 Appendix: For notational brevity in Algorithm 1, we use ∆J j (y l ) := J {(x, l)j } i-1 n+1 ∪ {(x, l)} - J {(x, l) j } i-1 n+1 .Algorithm 1 BAMS input: data points {(y l )i} N i=1 for all l ∈ {0, 1, . . . L}, posterior θn, number of clusters S, budget m, overbudget parameter η Create S clusters of points (Algorithm 2)for s = 1 to S costs ← 0 Queue s ← [∅] while costs < ηmNs/N y l ← arg min of optimization problem (5) in cluster s costs ← costs + c(y l ) Store (y l , ∆Jj(y l ), c(y l )) in Queue s costg ← 0 Setg ← {∅} while costg < m Candidates ← [ ∅ ] for s = 1 to S (y l , ∆Jj(y l ), c(y l )) ← next item in Queue s if costg + c(y l ) < m Candidates ← Candidates + (y l , ∆Jj(y l ), c(y l )) if Candidates = [ ∅ ] break (y l , ∆Jj(y l ), c(y l ))next ← arg min Candidates ∆Jj(y l ) costg ← costg + (c(y l ))next Setg ← Setg + {(y l )next}return the set of next points to evaluate SetgAs noted in Proposition 1 of Chevalier et al. [46] , we haveβ n x; {x j } n+m n+1 := E θn h n+m (x)|{X j = x j } n+m n+1 = Φ 2 s(x) -s(x) , t(x) 1 -t(x) 1 -t(x) t(x) ,wheres(x) := γ -µ n (x) σ n+m (x) , t(x) := σ 2 n (x) σ 2 n+m (x) = 1 + cov n (x, X m )cov n (X m , X m ) -1 cov n (X m , x) σ 2 n+m (x), and X m ∈ R m×d is the matrix of m new evaluation points. Departing from the form in Chevalier et al. [46] , we further simplify this result to a form that eliminates σ n+m . In particular, note that we can rewritet(x) := 1 t(x) = 1 - cov n (x, X m )cov n (X m , X m ) -1 cov n (X m , x) σ 2 n (x).We can then rewriteΦ 2 s(x) -s(x) , t(x) 1 -t(x) 1 -t(x) t(x) = Φ 2 ŝ(x) -ŝ(x) , 1 t(x) -1 t(x) -1 1 ,where ŝ(x) := γ-µn(x) σn(x) . By computing the Cholesky decomposition of cov n (X m , X m ) once (an O(m 3 ) operation), we can then compute t(x i ) for all i in O(m 2 N ) time by solving the triangular linear system. Then computing the expectation E X β n X; {x j } n+m n+1 is O(N ), resulting in an overall complexity of O(m 2 N ).We first outline the procedure for naive sequential selection followed by how we improve upon it.Naive sequential selection goes as follows: given q < m previously selected points, we select consider all N -n -q = O(N ) candidates for the next point. We compute the acquisition function (4) for all of them and choose the argument of the minimum, resulting in O(N ) iterations each of an O(q 2 N ) computation for a total of O(q 2 N 2 ) cost. Since we must do this for 1 ≤ q ≤ m -1, the resulting naive sequential selection strategy is O(m 3 N 2 ).We increase efficiency over the naive strategy by exploiting the fact that we are simply adding a single point to the selected candidates, so we do not need to blindly perform matrix inversion on the entire cov n (X q+1 , X q+1 ) -1 k(X q+1 , x). Instead, we employ identities for block matrix inversion.We introduce notation in this section purely for brevity; we do not use this notation anywhere else in the paper. Define the shorthand Σ n,q+1 := cov n (X q , X q ). Then we can write (with the right hand side as shorthand),Σ n,q+1 = Σ n,q cov n (X q , x q+1 ) cov n (x q+1 , X q ) cov n (x q+1 , x q+1 ) =: A f f gDefining shorthand cov n (X q+1 , x) := (k T q |k) T , we have via the block matrix inversion identitycov n (x, X q+1 )Σ -1 n,q+1 cov n (X q+1 , x) = k T q A -1 k q + 1 h (k T q A -1 f ) 2 - 2k h (k T q A -1 f ) + k 2 h where h := g -f T A -1 f. Note that the first term on the right hand side is precisely cov n (x, X q )Σ -1 n,q cov n (X q , x), so we can define a recursive relationship as long as we can also write a recursive relationship for k T q A -1 f . For an arbitrary vector (e T q |e) T we can utilize the block matrix inversion identity again to show that(e T q |e)Σ -1 n,q+1 cov n (X q+1 , x) = e T q A -1 k q + e T q A -1 f h - e h k T q A -1 f -k ,where the first term on the right-hand-side is precisely e T q Σ -1 n,q cov n (X q , x). Then we have a recursive formula for cov n (x, X q+1 )Σ -1 n,q+1 cov n (X q+1 , x) which allows us to efficiently compute t(x) iteratively. In particular, we can compute t(x i ) for all i in O(N 2 ) time, after which we must also compute the acquisition function ( 4) for all N -n -q candidates and choose the argument of the minimum, an O(N 2 ) computation. Thus, the overall time complexity is O(mN 2 ).First we divide P into S independent clusters. A natural choice for this task is spectral clustering [50, 51] , but this requires both computing and storing the entire posterior covariance matrix, which is expensive. Instead, we use the trained Matérn lengthscales at l = 0 from the GP model to rescale the input points x i and run K-means clustering on the N points; Euclidean distance with the rescaled dimensions is a cheap approximation for the Matérn kernel. Since each point x corresponds to L points y l , each resultant cluster has N s points at each of the L levels, with S s=1 N s = N (for full details, see Algorithm 2).Then, we solve the minimization problem (5) over the independent clusters, with each cluster given a budget proportional to its size ηmN s /N . Here, η ≥ 1 is an "overbudget" parameter such that we draw more points than necessary after pooling together outputs from all clusters (cf. Algorithm 1).The average complexity of K-means via Lloyd's algorithm is O(SN ) [52] ; it is significantly faster than computing the acquisition function, yielding an overall speedup as shown in Section A.4. When there are S equal-size clusters, the complexity of our sequential selection strategy is O( ηm/S N 2 /S). We can trivially take advantage of parallel computation in this approach as well.As a heuristic to prevent extremely small cluster sizes, we perform a slightly modified version of K-means clustering. In particular, we first perform K-Means with Ŝ > S clusters. followed by iterative removal of the smallest clusters by merging them with another cluster with the smallest Haussdorff distance. For shorthand, let the point z be the original data point x scaled by the Matérn lengthscales at l = 0 from the GP model. Let C i for the set of points in cluster with index i. Then, for points z a in cluster i and z b in cluster j, let the Euclidean distance be denoted as d(a, b). The Proof. Let the probability distribution function be ρ(x) such that X ρ(x)dx = 1. Then the variance satisfies:Var(p γ (θ n )) := E θn (p γ (θ n ) -E θn [p γ (θ n )]) 2 (a) = E θn (p γ (θ n ) -E X [p n (X)]) 2 = E θn (E X [g(X, θ n ) -p n (X)]) 2 = E θn X (g(x, θ n ) -p n (x))ρ(x)dx 2 ≤ E θn X (g(x, θ n ) -p n (x)) 2 ρ(x)dx X ρ(y)dy = E θn X (g(x, θ n ) -p n (x)) 2 ρ(x)dx = E θn E X (g(X, θ n ) -p n (X)) 2 (b) = E X E θn (g(X, θ n ) -p n (X)) 2 = E X [p n (X)(1 -p n (X))] = E X [h n (X)] ,where the inequality is the Cauchy-Schwarz inequality, and we assume Fubini's Theorem applies so that we may interchange expectations in equalities (a) and (b).Proof. The proof is immediate following Proposition 3.1. Namely,E θn Var(p γ (θ n+m ))|{X j = x j } n+m n+1 (a) ≤ E θn E X [h n+m (X)] |{X j = x j } n+m n+1 (b) = E X E θn h n+m (X)|{X j = x j } n+m n+1 = J {x j } n+m n+1, where inequality (a) follows from Proposition 3.1, and again we assume Fubini's Theorem applies so we may interchange expectations in equality (b).Here we provide experiments which explore how we chose (i) the number of clusters, (ii) IS parameter α and (iii) cost for the cheap fidelity for simulation in our main experiments.We explore the impact of varying the number of clusters on both computational time and recall achieved at a retention budget of 5p γ N after Batch 3. Figure 6 illustrates the relationship between the number of clusters and the recall achieved at this retention level, averaged over 6 random seeds.Based on these findings, our clustering approach can significantly speed up the process while maintaining an acceptable recall. We select S = 6 clusters for BAMS, achieving a substantial speedup of 17× compared to the non-clustering approach while maintaining competitive recall performance.In BAMS, the importance sampler is defined as p n (x) α . The parameter α plays a crucial role in aligning the IS with the objective of rate-informed discovery. We select α to minimize the relative variance of the rate estimate while maintaining acceptable recall for importance-sampled items. Figure 7 showcases the impact of changing α on both the relative variance of the rate estimate and the recall achieved by IS with budget of 5p γ N samples. Based on this analysis, we choose α = 2.5 as it minimizes the relative variance and achieves a satisfactory recall. To determine the optimal level l = 1 cost in the AV setting, we conducted a parameter-tuning study. We varied the number of AV planner rollouts used to calculate the Q 0.25 (mT T C [1] . . . mT T C [i] ) and evaluated BAMS recall at a retention budget of 5p γ N after Batch 3. The cost factor for level 1 was set to i/32, where i represents the number of rollouts used. We employed the same parameters as our main experiment, with m b = 15 and m 1 = 20, and averaged the results across 6 random seeds.Figure 8 depicts the recall achieved at retention level of 5p γ N for varying level 1 costs. The results demonstrate that BAMS benefits from utilizing more samples, even with the added noise introduced by using fewer rollouts at level 1. The model can effectively learn the noise characteristics and achieve improved performance. Based on these findings, we selected i = 5 (and c(y 1 ) = 5/32).Objective We consider a two-dimensional input X = R 2 and P is the empirical distribution of N = 20000 samples from a Gaussian distribution N (0, I). Denoting x [i] as the i th dimension of x ∈ X , our objective function is f(x) = x [0] -c, x [1]-c 1 , where c = 1.95 determines the centers (±c, c) of two diamonds for the failure modes. These modes are areas where f (x) ≤ γ with γ = 0.56, resulting in a rate p γ = 0.005. Figure 2b visualizes P and highlights failure modes.Multifidelity setup Level 0 represents f (x) without noise. Level 1 has a cost c(y 1 ) = 0.10 and adds Gaussian noise N (0, 0.1 2 ) to f (x). For the initial batch (Batch 1), we employ a budget of m 1 = 10 evaluations. Subsequently, for batches 2 and 3, we allocate a budget of m b = 5 evaluations across both fidelity levels. Figure 9 showcases BAMS's superior performance in discovering failure modes. It consistently achieves the highest recall across all batches, surpassing other methods. Multifidelity approaches (BAMS and MCM-GP) outperform their single-fidelity counterparts (BAS and MC-GP), highlighting the advantage of leveraging different fidelity levels for efficient discovery. While BAMS and BAS achieve comparable performance in rate estimation (Figure 10 ), BAMS stands out with the highest recall (Table 3 ), signifying its exceptional ability to discover failure modes while accurately estimating the rate. This demonstrates BAMS's superior performance in rate-informed discovery. In Figure 10 we utilize K = 2p γ N samples for IS rather than 5p γ N as in the AV experiments.Our setup is as follows: we implemented the Cross-Entropy (CE) method using a multivariate Gaussian sampling distribution with a diagonal covariance matrix. At every iteration, we drew samples from the Gaussian and found the nearest neighboring run segments according to the Euclidean embedding distance. From the samples in this batch, we used the samples with the 5 most dangerous TTC values to update the mean and diagonal covariance of the sampling distribution. After 3 batches (like all of our other methods), we used the final multivariate distribution function in the third batch to generate importance scores for all run segments and performed rate estimation with a sample size of 5p γ N like the other methods. Recall ±SE 100(RV ± SE) BAMS 1.00 ± 0.000 2.00 ± 0.003 BAS 0.917 ± 0.003 3.85 ± 0.366 MC 0.010 ± 0.001 710. ± 214. MC-GP 0.394 ± 0.004 39.4 ± 3.88 MCM-GP 0.537 ± 0.005 29.76 ± 5.08 Our work on rare event discovery and rate estimation complements several areas within the field of AV safety evaluation. For instance, while we focus on selecting informative run segments from a given set, other works like Feng et al. [55] address the orthogonal problem of efficient IS within a single run segment, often employing techniques like adversarial agents. Similarly, while methods like Arief et al. [56] rely on specific assumptions about the failure region's topology (e.g., orthogonal monotonicity), our approach remains agnostic to such constraints, enhancing its applicability to a wider range of scenarios.Our work aligns with the broader themes of scenario generation and safety evaluation. Existing surveys offer comprehensive overviews of scenario generation methods [57, 58, 59] , categorizing our approach under data-driven or learning-based techniques. Importantly, our method addresses the crucial challenge of transferability, highlighted in Ding et al. [57] , by focusing on efficiently identifying critical scenarios within a given logged dataset rather than relying on potentially brittle synthetic scenario generation.This work was motivated by the desire to develop a Bayesian formulation for Sequential Monte Carlo (SMC). While non-Bayesian SMC approaches exist [35, 37] , a Bayesian perspective offers several benefits. Our primary contribution, BAMS, embodies the novel Bayesian component of this broader research direction and is designed for seamless integration into SMC frameworks.A key advantage of a Bayesian formulation lies in reducing the sample complexity of Markov Chain Monte Carlo (MCMC) steps within the sequential levels of SMC. Specifically, by leveraging a GP, BAMS performs these MCMC steps directly within the embedding space, completely bypassing the need for simulations during this phase. This approach stands in contrast to previous suggestions in [35] like using surrogate simulations for a portion of the MCMC steps followed by a final simulationbased Metropolis-Hastings (MH) acceptance step. In particular, the Bayesian formulation requires no simulation for any of the MH steps, contributing significantly to its efficiency.Theoretical bounds on multilevel splitting performance We can quantify the intuition of MCMC being expensive by comparing the theoretical lower bound of simulations needed by (single fidelity) multilevel splitting to achieve the same relative variance as BAMS and BAMS. Using the notation in Norden et al. [37] , the number of iterations is K = log(p γ )/ log(1 -δ), the total simulation cost is N + T δN K, and the resulting relative variance is K δ N (1-δ) . Using the standard value of δ = 0.1 for splitting and the desired relative variance of 0.00851 (the value for BAMS after Batch 3), we can calculate K = 43, N = 561, and the total number of simulations is 561 + 2412T ≥ 2973, where the lower bound (T = 1) is achieved only if we use surrogate simulation to eliminate all but one of the iterations of simulation during MCMC. BAMS only has 2296 simulations, so vanilla multilevel splitting requires at least 30% more simulation cost for the same precision. Carrying out the same math for BAS (with relative variance 0.00969), we get N = 493 and multilevel splitting has to have at least 2613 simulations for the same precision. Fixing the budget for multilevel splitting to 2296 simulations, multilevel splitting has a lower bound of 100RV ≥ 1.10, which is worse than BAMS and BAS.How loose are these lower bounds for simulation cost given a fixed precision, and a fixed simulation cost? The discrete structure of our problem space (since we consider the base distribution P as a discrete distribution over N run segments) creates difficulty for implementing vanilla SMC baseline approaches. In particular, efficient MCMC techniques such as Hamiltonian Monte Carlo cannot be easily utilized. In this sense, the theoretical lower bounds are quite loose because it is likely that we need to perform more than one MH step using the real simulator. Furthermore, creating competitive SMC baseline approaches would require developing a novel graph-based MCMC kernel (to efficiently deal with the discrete structure); this is out of scope for this paper and may be a novel paper itself. Part of the motivation for this paper is to develop a Bayesian approach that sidesteps the necessity of performing difficult MCMC on a discrete space.Comparison with the cross-entropy method In addition to the advantages of BAMS over SMC methods, it also demonstrates superior performance compared to the cross-entropy method [47] , as detailed in Section 4.2. BAMS bears many similarities to the cross-entropy method. In particular, we have a sampling distribution (built upon the GP in our case) that updates as we iteratively sample. In BAMS, we utilize an optimization problem to perform sampling rather than simply drawing from a parametric sampler. 


Decision: The paper was initially evaluated with mixed reviews, but after the rebuttal, we the reviewers suggest accepting the paper. 

Here is a high level overview of the reviews. (before rebuttal)

### Strengths:
1. The multifidelity approach of BAMS is novel and provides a balance between computational cost and accuracy by combining different levels of simulation fidelity. (Reviewer EQTe)
2. The proposed acquisition function effectively guides the exploration and discovery of failure modes, improving recall and reducing variance in rate estimates. (Reviewer EQTe)
3. The experimental results show that BAMS outperforms existing methods in terms of recall and rate estimation accuracy in the multifidelity setting. (Reviewer EQTe)
4. The motivation for estimating AV performance under risky scenarios and discovering potential failure cases is promising and crucial for the autonomous driving field. (Reviewer GzcW)

### Weaknesses:
1. **Lack of Implementation Details**: Additional implementation details and open-source code availability would improve reproducibility. (Reviewer EQTe)
2. The paper does not provide a limitations section, which is essential for a comprehensive understanding of the work. (Reviewer EQTe)
3. **Insufficient Related Work Discussion**: The paper lacks a thorough discussion of related works in the literature, especially on rare-event simulation for autonomous driving and safety-critical driving scenario generation. (Reviewer GzcW)
4. There are many ambiguous terms and missing details in the paper, such as the definitions of certain concepts and the specifics of the multifidelity setting implementation. (Reviewer GzcW)

The reviewers find their concerns to be sufficiently addressed. The authors also improved organization an readability in the latest version of the paper. I suggest accepting the paper.

Reviews: review: title: title: Official Review of Submission379 by Reviewer EQTe
paper_summary: paper_summary: This paper introduces the Bayesian Adaptive Multifidelity Sampling (BAMS) algorithm for improving rate-informed discovery in the context of autonomous vehicle (AV) simulations. BAMS uses multifidelity simulations to efficiently identify diverse failure modes by combining high-fidelity and low-fidelity data to improve the discovery and estimation processes.
main_review: main_review: 1. Strengths/Weaknesses:

Strengths:

Innovative approach: The multifidelity approach of BAMS is novel and provides a balance between computational cost and accuracy by combining different levels of simulation fidelity.

Effective acquisition function: The proposed acquisition function effectively guides the exploration and discovery of failure modes, improving recall and reducing variance in rate estimates.

Strong empirical results: The experimental results show that BAMS outperforms existing methods in terms of recall and rate estimation accuracy in the multifidelity setting.

Weaknesses:

Additional implementation details and open-source code availability would improve reproducibility.

Limitations Section is not provided.

2. Clarity, Quality, Novelty And Reproducibility:

Clarity: The paper is well-written and clearly organized.

Quality: High-quality in terms of theoretical development and empirical validation.

Novelty: The multifidelity approach, the specific acquisition function, and the practical (clustering) relaxation used in BAMS are novel contributions to the field.

Reproducibility: While the theoretical aspects are well documented in the appendix, additional details on implementation and code availability would improve reproducibility.

3. Correctness:

The paper appears to be correct in its theoretical foundations and empirical validations.

The proofs and propositions presented support the claims made about the algorithm's performance.

4. Novelty And Significance:

Technical: BAMS represents an advancement in the field of autonomous vehicle simulation by introducing a method that efficiently combines multiple fidelities of data.

Empirical: The empirical results are compelling, showing that BAMS consistently outperforms existing methods in detecting diverse failure modes and providing accurate rate estimates.

This empirical strength underscores the practical utility and potential impact of the proposed method in AV development.
questions: questions: None
limitations: limitations: 1
review_summary: review_summary: Accept: The paper should be accepted for its novel approach, strong theoretical foundations, and impressive empirical results.
title: title: Review of Submission379
paper_summary: paper_summary: This paper aims to accurately estimate the performance of AV under risk scenarios and efﬁciently discover potential failure cases. To do so, the authors introduce Bayesian adaptive multiﬁdelity sampling (BAMS), which prioritizes exploration of regions with potentially low performance, leading to the identiﬁcation of novel and critical scenarios that traditional methods might miss. Previous work such as naive Monte Carlo (MC) is intractably expensive for rate-informed discovery. Furthermore, importance sampling (IS) that rely heavily on previously determined failure modes are poor at discovering novel issues. Compared with these two methods, BAMS discovers 10 times as many issues as MC and IS methods. In addition, BAMS also generates rate estimates with variances 15 and 6 times lower than MC and IS respectively.
main_review: main_review: **Strength**

1. The motivation for estimating AV performance under risky scenarios and discovering potential failure cases is promising. This is definitely an unsolved but crucial problem in the autonomous driving field.
2. This paper is mostly well-written and well-organized. It is easy to understand the core idea and main contribution.

**Weakness**

1. A lot of details are missing. Some examples can be found in my questions. Overall, I think some of the concepts mentioned in the paper may be familiar to audients with a statistical background rather than a robotics background. Some terms also have different meanings in these two domains, for example, embedding model and simulator.
2. Novelty is a little bit limited as estimating the rate of rare events in autonomous driving has been investigated in a lot of previous work. MC and IS are two well-known methods but I think there are some more recent works in this direction.
3. Based on the previous point, there is a lot of missing literature on rare-event simulation for autonomous driving, for example, [1-6]. I think the authors should thoroughly discuss the connection between the proposed method and these existing methods. Comparison in the experiment section is also necessary. 
4. In addition to the estimation of the rate of rare events, I think this paper also proposes to discover failure scenarios. Then, safety-critical driving scenario generation is also a large area that needs to be included in the related work section. For example, the authors can start with some survey papers [7-9].

---
[1] Feng, S., Sun, H., Yan, X., Zhu, H., Zou, Z., Shen, S., & Liu, H. X. (2023). Dense reinforcement learning for safety validation of autonomous vehicles. Nature, 615(7953), 620-627.

[2] Sinha, Aman, et al. "Neural bridge sampling for evaluating safety-critical autonomous systems." Advances in Neural Information Processing Systems 33 (2020): 6402-6416.

[3] Norden, Justin, Matthew O'Kelly, and Aman Sinha. "Efficient black-box assessment of autonomous vehicle safety." arXiv preprint arXiv:1912.03618 (2019).

[4] O'Kelly, M., Sinha, A., Namkoong, H., Tedrake, R., & Duchi, J. C. Scalable end-to-end autonomous vehicle testing via rare-event simulation. NeurIPS 2018.

[5] Arief, M., Huang, Z., Kumar, G. K. S., Bai, Y., He, S., Ding, W., ... & Zhao, D. (2021, March). Deep probabilistic accelerated evaluation: A robust certifiable rare-event simulation methodology for black-box safety-critical systems. In International Conference on Artificial Intelligence and Statistics (pp. 595-603). PMLR.

[6] Arief, M., Huang, Z., Kumar, G. K. S., Bai, Y., He, S., Ding, W., ... & Zhao, D. (2020). Deep probabilistic accelerated evaluation: A certifiable rare-event simulation methodology for black-box autonomy. arXiv preprint arXiv:2006.15722, 1(4).

[7] Ding, W., Xu, C., Arief, M., Lin, H., Li, B., & Zhao, D. (2023). A survey on safety-critical driving scenario generation—A methodological perspective. IEEE Transactions on Intelligent Transportation Systems, 24(7), 6971-6988.

[8] Zhong, Z., Tang, Y., Zhou, Y., Neves, V. D. O., Liu, Y., & Ray, B. (2021). A survey on scenario-based testing for automated driving systems in high-fidelity simulation. arXiv preprint arXiv:2112.00964.

[9] Schütt, B., Ransiek, J., Braun, T., & Sax, E. (2023, June). 1001 ways of scenario generation for testing of self-driving cars: A survey. In 2023 IEEE Intelligent Vehicles Symposium (IV) (pp. 1-8). IEEE.
questions: questions: 1.	Figure 1 shows the pipeline of a loop of finding adversarial scenarios and using it to improve the planner. However, I couldn’t find the stage of the planner improvement in this paper. 
2.	What does “reducing the end-to-end cost” mean in the last sentence of the conclusion section?
3.	In line 227, what is the planner used in the sentence “32 rollouts per run segment generated via AV planner simulation”? Why are these 32 rollouts different?
4.	In line 230, why use an embedding model to generate x? What is the specific definition of x rather than “the realization of the AV and its environment”?
5.	I am not sure about how the multi-fidelity setting is implemented during the experiment. Does it mean we use different planners to simulate the same initial positions, or do we use different simulators (e.g., Carla, highway-env, Waymax, …)? I think the terminology may lead to some ambiguity here. 
6.	What is the connection between “simulators” and “ﬁdelities”?
7.	In line 246, why use $\gamma=0.43$ as the threshold?
8.	In Figure 4, the authors show some examples of the discovered scenarios. The authors say that the scenarios discovered by BAMS showcase busy intersections, pedestrians, and cyclists. However, I think it is not enough to show that these scenarios are adversarial and interesting. It is also easy to use a rule-based filter to find those scenarios.
limitations: limitations: 2
review_summary: review_summary: My major concern is that the novelty of this paper is limited due to omit of a lot of related works. In addition, the paper needs to be improved to include more details about some design choices and terminologies.