Rate-Informed Discovery via Bayesian Adaptive Multifidelity Sampling
Authors: Aman Sinha
Payam Nikdel
Supratik Paul
Shimon Whiteson
Venue: CoRL 2024 (2024)
Keywords: Autonomous Driving
Rare-event Simulation
Adaptive Sampling
TL;DR: We present a method that efficiently estimates autonomous vehicle safety and performs targeted discovery of high-impact failure cases.

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.

Rate-Informed Discovery via Bayesian Adaptive Multifidelity Sampling
Aman Sinha', Payam Nikdel", Supratik Paul, and Shimon Whiteson Waymo, L.I.C
(thisisaman, payamn, supratikpaul, shimonu @waymo. com
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 simultancously estimating the rate of adverse events. BAMS prioritizes exploration of regions with potentially low performance, leading to the identifi- cation 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 Monte Carlo (MC) and importance sampling (IS) baselines, while the same time generating rate estimates with variances 15 and 6 times narrower than MC and IS baselines respectively.
Keywords: Autonomous Driving, Rare-event Simulation, Adaptive Sampling
1 Introduction
Commercial autonomous vehicle (AV) development usually follows an iterative process wherein issues with the current planner identified, root causes determined, and the issues 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: 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 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 expensive and does not ensure sufficient coverage of driving scenarios since we cannot control the environment [1]. Furthermore, since this upproach essentially produces 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 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, naive simulation-base 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 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:
Py := (X) S >) = Fx (X( S 7)]
where X ~ P, the function ʃ X + IR is performance metric scoring the realization x E X' of the AV and its environment, and 7 is 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
tractable IS techniques; these techniques place high probability on regions of X' where f(xx) ) is small.
However, simply estimating Phy (the objective of importance sampling (IS)) is insufficient, as simula- tion testing must provide actionable feedback on how improve AV design (Figure 1). This feedback comes in the form of concrete examples 24 where fxx)) < 7. In particular, during the continuous
Equal contribution
8th Conference on Robot Learning (CoRL. 2024), Munich, Germany.

---

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 aftributable 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 are used to understand and address the failure cases. E: An improved planner is deployed, and the cycle repeats.
B
feedback loop of simulation testing followed by design improvement, it is most useful for simulation to discover failure cases 2; with a) relatively high likelihood under P and b) novelty with respect to previously discovered failure examples. From a practical perspective, prioritization by likelihood obvious when trying to allocate resources fixing problems, and novelty is important because previously- failure modes typically have existing work underway to improve behavior.
We define aateiinormed 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 py, 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 targeted fashion and builds an importance sampler for estimating Py. 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.
2 Background and related work
Gaussian processes A Gaussian process (GP) [5] imposes prior belief on unknown random function 0 X C Rd + through its mean function p(x) (for convenience usually set to zero) and covariance function k(x, x²), such that any finite collection of function values {0(2,)} has a joint Gaussian distribution. A popular choice of the covariance function is the Matérn kernel: (xxxx') x') x t" (V2rt), where v lengthscale parameter (typically v 2.5), K2 the modified Bessel function of the second kind, and t {(x -')TE=' (2 x²)} is a distance sealed by diagonal matrix Σ of lengthscales for each dimension. Letting X,, E 1xxd be the matrix of observed inputs and 0 ((x) be the vector of associated function values, the posterior prediction ((x") for an input x² is Gaussian distribution N ((mn(" of"""""))))))))
(n(x") = k(x", Xx))(X,, (0)_))(Xn),
-z0())= = k(a",
k(0", (x), Xn)=x)(Xxx.2
The posterior covariance between any two inputs (x,) is cov, (2r, k(x,,²) k(x, X Xx,) - z'). There is a vast literature on speeding up GP inference, namely assuaging the O(n²) time complexity of inverting k(Xn, X... In particular, inducing point approx- imations [6, 7, 8] impose structure on the GP prior, which in turn induces (block) sparsity onto k(X,, X.). 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. 10J. 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 is highly budget-constrained and is small. In our experiments, 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.
simulate and use only tens of
Bayesian optimization and quadrature Bayesian optimization (BO) [11, 12] and quadrature (BQ) [13, 14, 15] use GPs for targeted function evaluation 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 fxz (see, e.g., Wilson al. [16] for an overview).
Multifidelity simulation There is a rich literature on using multiple fidelities for evaluation, also known as surrogate modeling (sce, 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 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 our acquisition function, especially when multiple fidelities are available for the same data point x.
easily incorporate costs of cach fidelity into
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 detailed discussion on how our Bayesian approach connects to existing rare-event simulation techniques, addresses key challenges within this domain, and contributes the broader field of AV safety evaluation.
3 Bayesian adaptive multifidelity sampling
In typical AV development environments, the distribution P over which compute the rate of adverse events (1) is an empirical distribution over run segments- segments-speeified time intervals (c.g.. 30 seconds) of logged real-world driving over which simulations are performed [39, 40J. We consider the empirical distribution over N run segments, and each run segment is encoded by an embedding
E Rdd There are 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 Bayesian estimator for the rate of adverse events (1). Specifically, we define random function 0 : X - I. and impose a GP prior over such that the collection {0(z:)} for any subset of run segments {xx)} has a joint Gaussian distribution. Considering set of evaluated run segments and their associated performance metries X, := {(x., f(x;)) { we define the shorthand On and Eg,, [-] as the posterior GP and the expectation with respect to this posterior induced by conditioning on the o-algebra generated by X,, We can then define a function which models the presence of adverse events based on our GP posterior:
g(X,0,, 0,) =- (X) < 7) {0, 1},
with marginals po(0)) := Ex [x(, 0,)] and (x) == Es, [9(2, 0,)]. The marginal (0,) is an estimator of the rate of adverse events (1).
Using 9 to model the performance metric f(x) of a run segment rather than the binary indicator 1[f(X) S > 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 7, so we can reuse 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(z) 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(ū (0)), is beneficial. In particular, minimizing (ar(p., (0) provides confident estimate of the probability (1). Furthermore, the calculation of Var(ý, (0) is dominated by regions in X that have both high likelihood under P and high uncertainty under 0n. Thus, minimizing Var(p, (0)) prioritizes sampling points x that are clustered near other points (high likelihood) and are far away from previously sampled points (high (nreetaintt)-rrciiisst our goal for rate-informed
discovery. However, this variance requires an expectation over the joint distribution X, x' Lid (see Appendix B));; because P is the empirical distribution over run segments, this O(N2) computation too expensive. Instead, we consider a cheaper upper bound.
3

---

First we define the point-vanne variance function (z) := (x(x) - (z)). The average of (x) over P. which requires only ((N) computation, is greater than the variance of the GP estimator:
Proposition 3.1. The variance of the estimator pyy (0,) ] is upper-bounded by the average point variance: Var(e, (0.)) I
See Appendix B.2 for the proof. Since we evaluate the next m points in batch setting, the order of these points does not matter, and our acquisition function is function of the set {2z]n4]* // In particular, consider the expectation of the point-variance function conditioning on the set of locations of these m future evaluations:
3r (z; (atr)) == Ee,, (x))(( xz]nti1)
This forward-looking point variance has a nested expectation Es, [Eon.m HX, I Jatr] ]]. the average over all future posteriors conditioned on the locations { kny]nti computation in GP models. Then, our acquisition function to choose the next m evaluation locations the average over P of this forward-looking point variance:
Jntm, this is an analytic
J (f2z) )= == Ex [&,, (X; (2y (et))t)))
(4)
An 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 pr, (0%) conditioned on locations {r))atT of the next evaluation points is upper-bounded by the acquisition function:
Eo., (Oxt)=))) = =z}nff] <J ((z))att)) 1=tm)
Sequential selection For m, & N. the acquisition function (4) can be efficiently calculated in O(m²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 > points is a combinatorial problem over (M) combinations. Since (w) = O(N""
this combinatorial problem has an overall time complexity of 0mm2~m Nm++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² N²2 time, but exploiting the recursive structure of the problem allows us to reduce this cost to O (m (see Appendix A.3 for details). In contrast, sequential selection using expected variance as an acquisition function requires Om..33) time.
In the sequential setup, our minimization objective for each iteration n +1+ +1 $ n m
minimize ((xx) U {x))= - J ((2))++1)
where {2, 3++1 := 0 and J(Ø) == Ex (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 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.
distilled version of the AV planner. Intuitively, a noisier but
For / 1.. L}. let the fidelities of simulation be denoted by fi, with the original high- siglliidliity (1=0) 0) simulator denoted fo- Consider the augmented input y1- := (xr I), which then allows to define the augmented simulation function f(v)) == fi(x). Next, let the model for each augmented input 0: X N + R be defined as (zx) := 6z(2)) ( (z) where 9, is the GP model over the function Jo. ogg : 0, and for 1 > 1, oz is an independent zero-mean GP with covariance function ky. Then, the augmented function 0 is GP with mean ((()) = u(z); ; the covariance between points ya = (I. a) and v/, (z', b) is: k(ya. 26) Cov(Ou(z) (z. (o((') + as(z') k(x, + (a,
where (a, 6) is the Kronecker delta function.
To update the acquisition function (4), we first define e(ur) as the cost associated with evaluating each simulator fi. In typical AV settings, this cost is expressed in units of money flops of computation.

---

Imputa
Clustering
Candidatee within clusters
Selected points
08
Update OP
Benca
postortor
(a) Schematic BAMS algorithm
(h) Synthetic data
(c) Selected points
Figure Illustration of BAMS. a) Iterative loop of our approach, defined in Algorithm (Appendix A. each iteration, we separate the inputs into clusters, solve problem (5) over each cluster. and then select the final points from the candidates in all the clusters. After performing simulations over these points, we update the GP posterior. For simplicity, we neglect illustrating multiple fidelities in this graphic. b) Distribution of samples P for simple synthetic setting, Samples are colored red if f(x) S Y and blue otherwise. c) Points selected over three iterations (batches) of our algorithm in the synthetic setting. Initial exploration over the state space in the first iteration transitions to targeted uncertainty reduction around the regions of interest in the third iteration.
We define c(yo) := 1 and by construction c(yi) < 1 for 2 1. Whereas in the single-fidelig setting we select a batch of m points, we now consider selecting 2 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 (z, l):
minimize (a,1)
J... (((=,)))77 O ((,)))) J (((z, 11)
(5)
where n n-+
Unlike in the single-fidelity case, the subtraction of the previous iteration's objective is now required-we we normalize the decrease of the forward-looking point-variance with respect the price of simulation.
selsiintee n n++1iss
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 2; at every level I, in which case simply define our distribution / over the empirical distribution of locations and corresponding available simulation fidelities. In Section 4, we assume that all simulators available for all N points.
Clustering As noted in Section 2, our problem setting 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 N. Instead, our bottleneck is the sequential selection of points to evaluate via minimization of the acquisition function, an OmmN²) computation. To decrease this cost, we employ clustering. With S clusters, the complexity of our sequential selection strategy is O(mN2 )52). 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 in Figure 2a. Figure 2b visualizes the input distribution /' 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 segments and 32 rollouts per 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 fir) experimental procedure, multifidelity setup, and evaluation criteria.
Embedding model The input embedding is generated in two stages:
1. 51 512-dimensional embedding: Using dataset of segments, we extract features such road layouts, other road users' positions, and AV kinematics. These features are used to generate top-down images, encoded into 312-dimensional vectors using 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.
5

---

2. Dimensionality reduction: We train an encoder-decoder model to reduce the embedding to 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 single time-to-collision (TTC) value for each run segment using the formula f(x) = (m.2s(mITT[n ... mw.C (z)), where m1TT (I) is the minimum TTC across the zth rollout for run segment with embedding x and (o.25 is the 25th percentile of the arguments. We use the 25th percentile to be robust to outliers the rollouts. Finally, we set the rate Pyy at 1 % by using the threshold 7 = 0.43 for this metric.
Experimental procedure We utilize BAMS with S = 6 clusters, based on the parameter-tunil study presented in Appendix C.I. We compare BAMS, which leverages both fidelity levels, with its ablated counterpart Bayesian adaptive sampling (BAS) which only employs the high-fidelit 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 difficuly-moodel scores.
DS-GP: Our method with IS based on dificuuly-model scores as the acquisition function.
CE: Cross-entropy method [47]. See Appendix 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² = 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 log- likelihood (MLL) using the Adam optimizer [48].
B. Two batches with budget my: Each method samples two additional batches, each with budget of m, = 15 evaluations. We retrain the GP byperparameters after each batch.
Multifidelity setup For BAMS and MCM-GP, we utilize two fidelity levels:
Level 0 (high fidelity): This level uses the 25th percentile TTC calculated from all 32 rollouts Q0.25 (m:TC .. mTTTC(32) with cost c(yo) = 1.0.
Level (low fidelity): This level uses the 25th percentile TTC calculated from only ran- domly selected rollouts, denoted Q0.25 (m]7CC mm.T'))) whereby c(y1) 5/32. This parameter was chosen based on parameter-tuning study (see Appendix C.3) conducted over 20.255mmTTC ... mTTC( ) with corresponding cost i/32, evaluating the model 's recall at retention budget of 5p
4.1 Evaluation criteria
We evaluate all approaches using two criteria: retention-recall curves and rate estimation/discovery via IS, detailed below.
Retention-recal 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 (I = 0) based on Pa (x;) estimates. For MC, samples are sorted randomly. Recall is then calculated different retention levels, representing the fraction of failures found within the top-ranked samples.
To facilitate interpretabilityy we scale the z-axis (retention) of retention- recall curves by N.N (number of failures). Thus, a value of on the z-axis corresponds 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 N,, all failure modes are captured within the top N sample scores Pn. (I;), precision is 0.5.
Retention-recaul 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; in the high-fidelity simulator, so it is not realistic way of measuring rate-informed discovery in AV development settings. We consider tractable statistics next.

---

BAMS
BAS
Recall
0.6
MC
MC-GP
Recall
0.6
Recall
0.6
MCM-GP
0.4
0.4
DS
DS-GP
CE
0-
0
10
10' Retention (P,N)
10"
10
10
16' Retention (p,N)
16f
10
10a
10'
10²
Retention
(p,N)
(a) Batch
(b) Batch 2
(c) Batch
Figure Retention-recall curves for AV experiments. The r-axis represents the retention budget scaled by Ny while the y-axis shows the recall at each corresponding retention budget.
Rate estimation and discovery via IS After completing the third batch of simulation, conduct combined evaluation of both rate estimation and discovery. We do 200 repeated trials using model scores to perform IS with budget of K = 5 high-fidelity samples (1 = 0). We calculate the resulting recall and relative variance (RV) of the rate estimate, defined as the variance divided by p² (e.g., the square of the coefficient of variation). For GP methods, define our importance sampler as Pn (x)", with a = 2.5 chosen to minimize relative variance (see Appendix C.2). For MC, sampling based on random scores.
This IS approach guarantees that our estimate of the failure rate py unbiased by construction. even if our underlying GP model 0 is biased. Together, these two statistics-relatiee variance and recall-test each method's ability to perform rate-informed discovery: 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.
4.2 Results
Figure 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 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 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, both MCM-GB GP and BAMS outperform their angle-fidelity counterparts. Notably, CE displays worse recall than all of the GP and multifidelity methods, but, CE 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 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 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 fz) values in the third batch (cf. Table 1), the low recall performance in Figure illustrates common pitfall of CE: 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 similar rate precision to BAS. These results reinforce the superior performance of multifidelity versions each model. Additionally, based methods demonstrate better recall and rate estimation compared to DS, CE, and MC.
7

---

Do00g 90300
BAMS
00
MCM-GP
00
0
60
We select five samples from Batch 3 data using DPP sampler maximize diversity. The vehicle is blue and the snapshot corresponds to the timestamp which the minimum TTC occurs; large transparent circles indicate the trajectory history, and small circles represent the future trajectory. The BAMS samples showcase busy intersections, pedestrians, and cyclists, whereas MCM-GP samples have less diversity.
Figure 4:
Methol
Batch
Batch
Batch
Table 1
The average f (x) values and standard errors (SE) for selected samples across three batches. Analogous to Figure 2c, the transition
BAMS
4.99 0.51
0.21
BAS
6.55 0.78
2.27 - 0.58
3.19 = 0.89
MC-GP
6.55
0.83
from exploration to targeted selection is most evi-
MCM-GP
1.99
6.09 = 0.51
5.86
+ 0.35
dent for BAMS.
DS-GP
6.55
7.37
CE
6.55 1: 0.78
2.91 1 0.26
2.20
0.12
Method
Recall LSE
10D(RV SE)
0.511 0.0013
evizaate
0.012
MC
BAS
0.388 + 0.0013
0.969 + 0.037
MC-GP
0.368 1. 0.0012
5.70 1 2.39
0.01
MCM-GP
0.0019
DS
0.051 1. 0.0007
5.40 0.096
0.008
CE
0.096 1: 0.0007
4.43 0.063
BAMS BAS
Mc
Mc-OP
MCM
D5-OP
Figure 5: AV setting 90% confidence intervals for IS rate estimate with K = 5p,, samples. The ground truth is py = 0.01. Smaller emror bars indicate a more precise estimator.
Table 2: evaluation
experiments. Recall relative variance of rate estimates are calculated after Batch 3 with K = 5p..N samples. Standard етror (SE) is over 10 seeds.
Limitations
While BAMS demonstrates substantial improvements for efficiently discovering novel, high- likelihood scenarios and estimating the rate of adverse events, our approach has limitations. Specif- ically, ensure tractable GP inference, we must perform dimensionality reduction our input embeddings; this pre-training step adds computational overhead. Furthermore, the effectiveness of BAMS hinges the availability of continuous (or at least real-valued) metric function f()); 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.
6 Conclusion
This paper introduced rate- vate-informed discovery to address the critical need for efficient estimation of AV safety and discovery of high-impact failure cases. We proposed BAMS, 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 simultancously estimating the rate of adverse events. In future work, we plan to "close the loop" by embedding inneoffeeeddd discovery via BAMS in iterative AV development cycle and measure its effectiveness at reducing the end-to-end cost and latency of planner improvement.
8

---

References
[1] Koopman and M. Wagner. Challenges in autonomous vehicle testing and validation, SAE International Journal of Transportation Safery, 4(1):15-24, 2016.
[2] N. Kalra and S. M. Paddock. Driving safety: How many miles of driving would it take to demonstrate autonomous vehicle reliability? Transportation Research Part A: Policy and Practice, 12211829199 2016.
[3] M. Scanlon, K. D. Kusano, Daniel, C. Alderson, A. Ogle, and T. Victor. Waymo simulated driving behavior in reconstructed fatal crashes within autonomous vehicle operating domain. Accident Analysis & Prevention, 163: 16:110445, 2021.
[4] Z. Yang. Y. Chen, J. Wang, S. Manivasagam, W.-C. Ma, J. Yang, and R. Urtasun. Unisim: A neural closed-loop sensor simulator. In Proceedings of the IEEE/CVF Conference on Computer Vision and Puttern Recognition, pages 1389-1399, 2023.
[5] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005.
[6] V. Tresp. A bayesian committee machine. Neural computation, 12(11)2271--244) 2000. 2000.
[7]
Ghinonero-Cander and E. Rasmussen. A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Research, 6:1999-1959, 2005.
[8] M. Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR, 2009.
[9] A. Wilson and H. Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In International conference on machine learning, pages 1775-1784. PMLR, 2015.
[10] G. Pleiss, J. Gardner, K. Weinberger, and A. G. Wilson. Constant-time predictive distributions for gaussian processes. In International Conference on Machine Learning, pages 4114-4123. PMLR, 2018.
[11] B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. De Freitas. Taking the human out of the loop: A review of bayesian optimization. Proceedings of the IEEE, 148, 2015.
[12] P. I. Frazier. A tutorial on bayesian optimization. arXiv preprint arXiv: 1807.0280 2018.
[13] A. O'Hagan. Bayes-hermite quadrature. Journal of Statistical Planning and Inference, 1991.
[14] M. Kennedy. Bayesian quadrature with non-normal approximating functions. Statistics and Computing, 8(4), 1998. URL https: http//ii..og///////// 1023/1:2332// 1011:::2200889220000
[15] C. E. Rasmussen and Z. Ghahramani. Bayesian monte carlo. Neural Information Processing Systems, 2003.
[16] Wilson, Hutter, and M. Deisenroth. Maximizing acquisition functions for bayesian optimization. Advances in neural information processing systems, 31, 2018.
[17] M. C. Kennedy and A. O'Hagan. Predicting the output from complex computer code when fast approximations are available. Biometrika, 87(1):1-13, 2000.
[18] V. Balabanov and G. Venter. Multi-fidelity optimization with figh-fideitty analysis and low- fidelity gradients. 10th AIAAVISSOO Multidisciplinary Analysis and Optimization Conference, page 4459, 2004.
[19] N. V. Queipo, R. T. Haftka, W. Shyy, T. Goel, Vaidyanathan, and P. K. Tucker. Surrogate- based analysis and optimization. Progress in aerospace sciences, 41(1))1-28, 2005.
[20] R. Lam, D. L. Allaire, and K. E. Willcox. Multifidelity optimization using statistical surrogate modeling for non-hieranee con-hicrarchical information sources. In 56th ALAA/ASCEAAHSSSCC Structures, Structural Dynamics, and Materials Conference, page 0143, 2015.

---

[21] L. Le Gratiet and C. Cannamela. Cokriging based sequential design strategies using fast coossvvlidattonn techniques for multi-fidelity computer codes. Technometrics, 57(3):418-422 2015.
[22] E. V. Bonilla, K. Chai, and C. Williams. Multi-task gaussian process prediction. Advances in neural information processing systems, 20, 2007.
[23] K. Swersky, J. Snoek, and R. P. Adams. Multi-task bayesian optimization. Advances in neural information processing systems, 26, 2013.
[24] A. Klein, S. Falkner, Bartels, P. Hennig, and F. Hutter. Fast bayesian optimization of machine learning hyperparameters on large datasets. In Artificial intelligence and statistics, pages 528-536. PMLR, 2017.
[25] O. Lübsen, C. Hespe, and A. Eichler. Safe multi-task bayesian optimization. arXiv preprint arXiv:23/2.0728// 2023.
[26] A. Marco, F. Berkenkamp, P. Hennig, A. P. Schoellig, A. Krause, S. Schaal, and S. Trimpe. Virtual vs. real: Trading off simulations and physical experiments in reinforcement learning with bayesian optimization. In 2017 IEEE International Conference Robotics and Automation (ICRA), pages 1557-1563, IEEE, 2017.
[27] M. Poloczek, J. Wang, and P. Frazier. Multi-informatiin source optimization. Advances in neural information processing systems, 30, 2017.
[28] M. E. O' Kelly, H. Abbas, S. Gao, Kato, Shiraishi, and R. Mangharam. Apex: Autonomous vehicle plan verification and execution. Technical report, SAE Technical Paper, 2016.
[29] M. Kwiatkowska, G. Norman, and D. Parker. Prism 4.0: Verification of probabilistic real-time systems. In Computer Aided Verification: 23rd International Conference, CAV 2011, Snowbird, UT, USA, July 14-20, 2011. Proceedings 23, pages 585 591.599..333 Springer, 2011.
[30] M. Althoff and M. Dolan. Online verification of automated road vehicles using reachability analysis. IEEE Transactions on Robotics, 30(4):903-918; 2014.
[31] S. A. Seshia, D. Sadigh, and S. S. Sastry. Formal methods for semi-autonommuus driving. In Proceedings of the 52nd Annual Design Automation Conference, pages 1-5, 2015.
[32] C. E. Tuncali, T. P. Pavlic, and G. Fainekos. Utilizing s-taliro as an automatic test generation framework for autonomous vehieles. In 2016 ieee 19th international conference on intelligent transportation systems (itsc), pages 1470-1475. IEEE, 2016.
[33] R. Lee, O. Mengshoel, A. Saksena, R. W. Gardner, D. Genin, J. Silbermann, M. Owen, and M. J. Kochenderfer. Adaptive stress testing: Finding likely failure events with reinforcement learning. Journal of Artificial Intelligence Research, 69: 69:1166 16651111,,, 2020.
[34] M. Koren, S. Alsaif, R. Lee, and M. J. Kochenderfer. Adaptive stress testing for autonomous vehicles. In 2018 IEEE Intelligent Vehicles Symposium (IV), pages 1-7. IEEE, 2018
[35] A. Sinha, M. O'Kelly, R. Tedrake, and J. C. Duchi. Neural bridge sampling for evaluating safety-critical autonomous systems. Advances in Neural Information Processing Systems, 33: 2020.
[36] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer. Reluplex: An efficient smt solver for verifying deep neural networks. In Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30, pages 97-117. Springer, 2017.
[37] Norden, M. O' Kelly, and Sinha. Efficient -ox assessment of autonomous vehicle safety. arXiv preprint /0/i::/9/.006/// 2019.
[38] A. Corso, R. Moss, M. Koren, R. Lee, and M. Kochenderfer. A survey of algorithms for black- box safety validation of cyber-physical systems. Journal of Artificial Intelligence Research, 72: 377-428, 2021.
10

---

[39] E. Bronstein, S. Srinivasan, S. Paul, A. Sinha, M. O'Kelly, P. Nikdel, and S. Whiteson. Embed- ding synthetic off pooicy experience for autonomous driving via zero-shot curricula. In K. Liu, D. Kulic, and J. Ichnowski, editors, Proceedings of The 6th Conference on Robol Learning, volume 205 of Proceedings of Machine Learning Research, pages 188-198. PMLR, 14-18 Dec 2023. URL https: http///oooceeenngeee/eeeeeee mlr, a.l..pess// wwwpsereeeaa w2205//boonteeinneeee.
[40] N. Montali, J. Lambert, P. Mougin, A. Kuefler, N. Rhinehart, M. Li, C. Gulino, T. Emrich,
Z. Yang. S. Whiteson, B. White, and D. Anguelov. The waymo open sim agents challenge. In Advances in Neural Information Processing Systems Track on Datasets and Benchmarks, 2023.
[41] R. Wickramarachchi, C. Henson, and A. Sheth. An evaluation of knowledge graph embeddings for autonomous driving data: Experience and practice. arXiv preprint arXiv:2003. 2000:0000300 2020.
[42] C. Choi, J. H. Choi, J. Li. and S. Malla. Shared cross-modal trajectory prediction for au- tonomous driving. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 244 253, 2021.
[43] S. Segal, E. Kee, W. Luo, A. Sadat, E. Yumer, and R. Urtasun. Universal embeddings for spatio-temporal tagging of self-driving logs. In Conference on Robot Learning, pages 973-983. PMLR, 2021.
[44] A. Hu, L. Russell, H. Yeo, Z. Murez, G. Fedoseev, A. Kendall, J. Shotton, and G. Corrado. Gaia-1: A generative world model for autonomous driving. arXiv preprint arXiv:2309.17000 17080, 2023.
[45] N. Nayakanti, R. AI-Rfou, A. Zhou, K. Goel, K. S. Refaat, and B. Sapp. Wayformer: Motion forecasting via simple & efficient attention networks. 2023 IEEE International Conference on Robotics and Automation (ICRA), pages 2980-2987. IEEE, 2023.
[46] C. Chevalier, J. Bect, D. Ginsbourger, E. Vazquez, V. Picheny, and Y. Richet. Fast paral- lel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set. T'ccnometrics, 56(4):455-46, 2014.
[47] R. Y. Rubinstein and D. P. Kroese, The cross-entropy method: unified approach to combi- natorial optimization, Monte-Carla simulation, and machine learning, volume 133. Springer, 2004.
[48] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980. 2014.
[49] A. Kulesza, B. Taskar, al. Determinantal point processes for machine learning. Foundations and Trends in Machine Learning, 5(2-3):123-288) 2012.
[50] J. Shi and J. Malik. Normalized cuts and image segmentation. IEEE Transactions on pattern analysis and machine intelligence, (188)8888905555 2000.
[51] A. Damle, V. Minden, and L. Ying. Simple, direct and efficient multi-way spectral clustering. Information and Inference: A Journal of the IMA, 8(1):181 --0880111000001111 06 2018. ISSN 2049-8772. 10i:100193//maaii//ee
URL https://doi.oooo org/10. 1et111999mmm/t htt///mmaaaiiay0......1
[52] C. D. Manning, P. Raghavan, and H. Schütze. Xml retrieval. Introduction Information Retrieval, 2008.
[53] Z. Drezner. Computation of the bivariate normal integral. Mathematics of Computation, 32 (141):277-279, 1978.
[54] Z. Drezner and G. O. Wesolowsky. On the computation of the bivariate normal integral. Journal of Statistical Computation and Simulation, 35(1-2): 1990.
[55] S. Feng, H. Sun, X. Yan, H. Zhu, Z. Zou, S. Shen, and H. X. Liu. Dense reinforcement learning for safety validation of autonomous vehicles. Nature, 615(7953):620 627,6 2023.
[56] M. Arief, Z. Huang. G. K. S. Kumar, Y. Bai, S. He, W. Ding. H. Lam, and D. Zhao. Deep probabilistic accelerated evaluation: A certifiable rare-event simulation methodology for black- box autonomy, arXiv preprint Xiv:2000611522. 1(4), 2020.

---

[57] W. Ding, C. Xu, M. Arief. H. Lin. B. Li, and D. Zhao. A survey on safety-critical driving scenario generation- a methodological perspective. IEEE Transactions on Intelligent Transportation Systems, 24(7)6991-69888 2023.
[58] Z. Zhong, Y. Tang, Y. Zhou, V. d. O. Neves, Y. Liu, and B. Ray. A survey on scenario-based test- ing for automated driving systems in high-fidelity simulation, arXiv preprint (...v:2112.009900 2021.
[59] B. Schütt, J. Ransiek, T. Braun, and E. Sax. 1001 ways of scenario generation for testing of self-driving cars: A survey. In 2023 IEEE Intelligent Vehicles Symposium (IV), pages 1-8. IEEE, 2023.
12

---

A Algorithmic details
A.1 BAMS algorithm
For notational brevity in Algorithm 1, we use AJj(m) = (((,,)))) ((,)))) (1(x,1))) 11
Algorithm 1 BAMS
data points {wr)) for all I E {0, T.), posterior On, number of clusters S, budget m, overbudget parameter
input:
Create S clusters of points (Algorithm 2)
for 1 to S
cost, -
Queue, <- [0]
while cost, < [ymnN./N]]
Ut < arg min of optimization problem (5) in cluster cost, f- cost, + c(yt)
Store (u,, (3n) in Queue,
cost,
Set, - {0) while cost, < m
Candidates + [0]
for = 1 to
(yr, (3r), c(w:) + next item in Queue,
if costs + c(vi) <
Candidates Candidates (ye, (J)( c(ui)
Candidates = [0 break
(yi, (3n), <(()))c(())))))))))))) Inea <- arg candicas A
cost, + cost, + (c(())) nest
Sel, <- Set, + )ass)
return the set of next points to evaluate Set,
A.2
Computation of Equation (4)
As noted in Proposition of Chevalier et al. [46], have
Pn (x; (z, (etr) == Es,, (0x+m(x)l((, = x, 1u]
=2 ( (~((z1) (1 (1-1(a) ((x)
t(x)
1 1- 1-t(x))
((x) )
where
s(x) = 7 (x)
Ontm (xr)
o (z)
cov, COv, (2, Xm cov, (Xm, Xmm) ( (Xm x.
t(x) = om (I)
0 o (x)
and Xm E Randd is the matrix of m new evaluation points. Departing from the form in Chevalier al. [46], further simplify this result form that eliminates ontm- In particular, note that we can rewrite
cov, (x, X,, )oov, (X,,,, Xm,) -ooo,, (X,,,
f(x) = t(x)
(²((r
We can then rewrite
22 ((5622)
s(x) (
t(x)
1- )
((=tt))
8(z)
f(x) (
1(x))
t(x)
$2
-
where s(z) = (=H.(2) By computing the Cholesky decomposition of cov, (Xm. Xm) once (an ((m3) operation), we can then compute f(z;) for all in O(m² N) time by solving the triangular linear system. Then computing the expectation Ex [8n (X; {xy Jntr)] O(N), resulting in an overall complexity of Om)
(n(x))
A.3 Recursive next point selection for expected point variance
We first outline the procedure for naive sequential selection followed by how we improve upon it.
13

---

Naive sequential selection goes as follows: given q < previously selected points, we select consider all N n q = ((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 ((N) iterations each of an Og computation for a total of O 22) cost. Since must do this for S S 1, the resulting naive sequential selection strategy is (m²²N))
We increase efficiency over the naive strategy by exploiting the fact that we are simply adding a single point the selected candidates, so we do not need to blindly perform matrix inversion on the entire cov,, (X ++1, x (+++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 En.tt1 = cov, (Xg: Xg). Then we can write (with the right hand side as shorthand),
Engg COVn ( Xg)
(X. ) (
Em =
covn (Zg+1
Defining shorthand cov, (X x)) = 1k)? we have via the block matrix inversion identity
2k
A22 (E) h h
cov.. Xx+1)
E n. cOV, (X 1) = k A + (k2 A-)
where h == cov,
A Note that the first term the right hand side is precisely Cov (z, x,)= cov,, n (Xg: 2), so we can define a recursive relationship as long as we can also write recursive relationship for k. A f. For an arbitrary vector (e]le)" matrix inversion identity again to show that
le)" we can utilize the block
T A'' 9
cov, (X)=1; Co = = -ii
A
(K. -
h
where the first term on the gight-hand-side is precisely formula for
eT _m covn cov, xx.
Then we have a recursive
co,,,.; Xx+1 ) ) n cov, (...v..(XX..... (x+1,
which allows us efficiently compute i(x) iteratively. particular, we can compute i(xr.) for all i in O(N2) time, after which we must also compute the acquisition function (4) for all N - n candidates and choose the argument of the minimum, an O(N2) computation. Thus, the overall time complexity is (m..3).
A.4 Clustering algorithm
First we divide P into S independent clusters. A natural choice for this task spectral clustering [50, 51], but this requires both computing and storing the entire posterior covariance matrix, which is expensive. Instead, input points v; and run K 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
use the trained Matérn lengthscales at 1 = 0 from the GP model to rescale the
points yt, each resultant cluster has N,, points at each of the L levels, with Σ N₂ 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 [Tm,,,/N Here, 17 2 1 is an "overbudget" parameter such that 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 5 equal-size clusters, the complexity of our sequential selection strategy is ([1m//5 N2 SS). We can trivially take advantage of parallel computation in this approach 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 > 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 be the original data point æ scaled by the Matérn lengthscales at / = 0 from the GP model. Let C, for the set of points in cluster with index i. Then, for points 2. in cluster i and 26 in cluster j, let the Euclidean distance be denoted as d(a, b). The
14

---

Algorithm 2 Clustering with Haussdorff merges
input: data points {(2x }e desired number of clusters S, initial number of clusters $
Create clusters of points via K-means
for i = 1 to $
arg m [c
arg ming (C, Ck)
C - uC.
Remove C.
return remaining clusters {C.)
Haussdorff distance between clusters and j is defined as:
du(C;,Cy) Cy) == max
max min d(a, b) max min d(a, b) aEC, bEC,
(6)
bEC, aEC,
The clustering procedure is shown in Algorithm 2.
B
Technical Results
B.1
Computation time for variance
The variance satisfies:
Var(p,, (0.)) Eon (0,(0m) Ea, (o.))))))
(a) Eo..
(p,(0n) Ex (x)))
Ea,. [Ex (x(, 00)) - Pm (X)))²
Foo, [Fxxx' [(g(X, ((((X,,)) 0n) - Pr (X)) (a(X', (0'())X')))))) (x()
Ex.x' [Es, ((9 b,) Px)))) (9(X', 0xx) PX)(X
Ex,x' (Cov(y(X,0,) 0, ,x),,
Σ
Covlg(æ 0n), (13, 0w)).
6-=1
where we assume Fubini's ' Theorem applies so that we may interchange expectations in equalities (a) and (b), and equality (c) is due to the fact that P is the empirical distribution function over N samples zi. The bivariate covariance above can be written as
Cov(g(X, (il,(X,,,),,te 9 (m)) Es. [9(X, 8,, (X X' 0,)] Pr, (X)p,,
P (0,, (X) 7, G (X') 7) Pm (X)),,(X))
42 7 on (X)
Un (X) 7 - (X') Cov, (X,
Pr. (X)),
(x'(') Tn (X)on(())
where a2 b,, is the cumulative distribution function for vector (a, b) under the bivariate normal distribution
N ())
We implement an efficient numerical procedure for this bivariate normal integral via the Gaus-
sian quadrature method outlined by Drezner [53], Drezner and Wesolowsky [54]. Computing Var(ū,(0..)) requires (+)) of these computations for i + j. When
j. can simplify
Cov(g(a,, 0w), g( 0n)) P.. (1:)(1 - Pv. (a:)), and p. (x) = + (,(()) where ((.) is the cumulative distribution function for standard (one-dimensional) normal. Overall, the computation time is dominated by the (?)) evaluations of 42, the complexity is O(N²).
15

---

B.2
Proof of Proposition 3.1
Proof.
Let the probability distribution function be p(r) such that Sx p(z)dz = Then the variance satisfies:
Var(ū)((..)) Egn. (0,) Eg,. (0,))))
(a) Eu,.. (10(00)) Ex 1pn (X)))²
Een.
(Exl9(X, 0n) - Pm (x)1
27
CoLo, 0n) - (=))))
Bon.
x (g(2,0n) On) -Pa(x))) p(z)da / p(uidy)
s Bon
I (g(2, P. (x)) p(o))d
Ee..
Fo..
[Ex [[o(X,00 0n)) (x)]
(b) Fxx [Fo,, [(g(X,0m) ) [X]
Ex (Px(() (x))]
Ex [4n. [x,),
where the inequality is the Sauchy-Schward inequality, and we assume Fubini's Theorem applies so that we may interchange expectations in equalities (a) and (b).
B.3 Proof of Corollary 3.2
Proof. The proof is immediate following Proposition 3.1. Namely,
Eg. (On = x, "
Eg. [Ex [hn+m (X)] = 2 [
Ex [Eo, [ (X) (X, 22 2,>]tr]
J ({x,) ((ry)nt?))
where inequality (a) follows from Proposition 3.1, and again we assume Fubini's Theorem applies so we may interchange expectations in equality (b).
C
Parameter-tuning studies
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.
C.1 Number of clusters
We explore the impact of varying the number of clusters on both computational time and recall achieved at retention budget of 5p,NN after Batch 3. Figure 6 illussrates the relationship between the number of clusters and the recall achieved at this retention level, averaged over random seeds.
Based on these findings, our clustering approach can significantly speed up the process while maintaining acceptable recall. We select S = clusters for BAMS, achieving substantial speedup of 17x compared to the non-clustering approach while maintaining competitive recall performance.
C.2 IS parameter
In BAMS. the importance sampler is defined as (x)x)" " The parameter a plays crucial role in aligning the IS with the objective of rate-informed discovery. We select minimize the relative variance of the rate estimate while maintaining acceptable recall for importance-sampled items. Figure showcases the impact of changing on both the relative variance of the rate estimate and the recall achieved by IS with budget N samples. Based on this analysis, we choose a = 2.5 as it minimizes the relative variance and achieves satisfactory recall.
16

---

0.85
1200
0.8
Time
1000
0.75
800
Recall
Time (8)
600
0.65
400
0.6
0.55
200
Number of clusters S
Figure 6: Exploring how the number of clusters affects both recall and computational time of BAMS's sequential selection stage N retention. Recall is shown in blue, time in red.
0.8
Recall
0.7
RV
0.6
3.5
0.5
Rocall
2.5
1.5
0.1
0.5
0.5
1
1.5
2.5
3
3.5
Figure Influence of IS parameter on relative variance of rate estimate and recall achieved using budget of Np.N samples. Recall shown blue, red.
C.3 Multifidelity cost
To determine the optimal level 1 - cost in the AV setting, we conducted a parameter-tuning study, We varied the number of AV planner rollouts used to calculate the Q0.25 (mTTCp mT and evaluated BAMS recall at retention budget of N after Batch 3. The cost factor for level I i/32, where i represents the number of rollouts used. We employed the same parameters as our main experiment, with m, = 15 and mn = 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 these findings, we selected = 5 (and c(y1) 5/32).
D Synthetic Experiments
Objective We consider a wwo-dimensional input X' R and is the empirical distribution of N = 20000 samples from Gaussian distribution N N/(0, I). Denoting zμ as the 7h dimension of E X' our objective function is f(z) [[(=)]]] c, Z[p - Il, where
1.95 determines the centers (#c. c) of two diamonds for the failure modes. These modes areas where f() < 7 with
7 0.56, resulting in a rate PA = 0.005. Figure 2b visualizes P and highlights failure modes.
Multifidelity setup Level 0 represents ((x without noise. Level I has a cost c(y) ) = 0.10 and adds Gaussian noise [0, 0.13) to fz). For the initial batch (Batch 1), employ budget of m, 10 evaluations. Subsequently, for batches 2 and 3, we allocate a budget of mn = 5 evaluations across both fidelity levels.
17

---

0.85
0.8
0.75
Recall
0.65
0.6
0.55
Loweliieeity cost c(jt))
Figure Impact of varying the cost of the lowest fidelity level (level 1) on recall N retention after Batch The data points represent the average recall across random seeds.
HAMS
0.8
BAS
U.S
0.5
MC
Recall
0.6
MI-CP
Recall
0.6
Recall
0.6
MCM-GP
0.2
0.2
0.2
0
0
10
10 Retention (P,N)
10P
10
10a1
10 Retention (P,N)
10'
10
10
10'
10'
10²
Retention (p,N)
Batch
(b) Batch
(c) Batch 3
Figure 9: retention-recall curves for synthetic experiments. The æ-axis represents the retention budget scaled by P. N while the g-axis shows the recall at each corresponding retention budget.
Results To illustrate the performance of our acquisition function, Figure 2c visualizes the samples selected by BAS. For clarity in the plot, we use budget my 50 (and we keep mi = 10). After the initial 10 random samples, the acquisition function guides exploration towards the top-left and top-right areas in Batch 2, revealing potential regions of interest. Finally, in Batch boundaries of the set of failure modes, precisely targeting those points with the highest uncertainty. This visualization emphasizes the effectiveness our acquisition function in progressively refining its focus, transitioning from initial exploration to targeted discovery.
samples the
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), highlight- ing 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 discover failure modes while accurately estimating the rate. This demonstrates BAMS's superior performance in rate-informed discovery. In Figure 10 utilize K N samples for IS rather than 5p.N as in the AV experiments.
E
Setup for the cross-entropy method
Our setup is follows: implemented the Cross-Entropy (CE) method using multivariate Gaussian sampling distribution with 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 most dangerous TTC values to update the mean and diagonal covariance the sampling distribution. After 3 batches (like all of our other methods), to generate importance scores for all run segments and performed rate estimation with a sample size of 5p,,N N like the other methods.
used the final multivariate distribution function in the third batch
18

---

0.03
Method
Recall +SE
100 SE)
edtinute
0.025
BAMS
1.00 + 0.000
2.00
0.008
BAS
0.917 0.003
3.85
0.015
MC
0.010 0.001
710. 214. 214.
0.01
MC-GP
0.394 0.004
39.4
0.005
MCM-GP
0.537
+ 0.005
29.76
BAMS
BAS
MC
MC-GPP
MCM-GP
Figure 10: Synthetic setting 90% confidence intervals for IS rate estimate with K 2 samples.
Table evaluation for synthetic experiments.
Recall and relative variance rate estimates are calculated after Batch 3 using IS with K = 2p, samples. Standard error (SE) is calculated over 10 seeds.
F
Connections to other rare-event simulation and AV safety-evaluation
techniques
Our work on rare event discovery and rate estimation complements several areas within the field of AV safety evaluation. For instance, while we focus selecting informative segments from 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 al. [56] rely on specific assumptions about the failure region's 's topology (e.g., orthogonal monotonicity), our approach remains agnostic to such constraints, enhancing its applicability to 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, 591, categorizing our approach under data-driven or tearning-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 potentially brittle synthetic scenario generation.
This work was motivated by the desire to develop Bayesian formulation for Sequential Monte Carlo (SMC). While non-Bayesian SMC approaches exist [35, 37J. a Bayesian perspective offers several benefits. Our primary contribution, B BAMS, embodies the novel Bayesian component of this broader research direction and is designed for seamless integration into SMC frameworks.
key advantage of Bayesian formulation lies in reducing the sample complexity of Markov Chain Monte Carlo (MCMC) steps within the sequential levels of SMC. Specifically, by leveraging 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 final simulation- based Metropolis- Cetropolis-Hasting 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 6). the total simulation cost is
N T8N K, and the resulting relative variance is 5) Using the standard value of 6 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 + 24127 > 2973, where the lower bound (T = 1) 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 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 2 1.10, which is worse than BAMS and BAS.
19

---

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 discrete distribution over run segments) creates difficulty for implementing vanilla SMC baseline approaches. In particular, efficient MCMC techniques such Hamiltonian Monte Carlo cannot be easily utilized. this sense, the theoretical lower bounds are quite loose because it is likely that we need perform more than one MH step using the real simulator. Furthermore, creating competitive SMC baseline approaches would require developing 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 to develop 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 sampling distribution (built upon the GP in our case) that updates BAMS, we utilize an optimization problem to perform sampling rather than simply drawing from parametric sampler.
iteratively sample. In
20


Reviews:
[{"review_id": "uiPMwTAu5b", "reviewer": "Reviewer_EQTe", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 3, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "Accept: The paper should be accepted for its novel approach, strong theoretical foundations, and impressive empirical results.", "technical_quality": 4, "clarity_of_presentation": 4, "potential_impact": 3, "robotics_focus": 2, "quality_of_limitations": 1, "questions_for_rebuttal": "None", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8hYy0wblLr", "reviewer": "Reviewer_GzcW", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": 2, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "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.", "summary_of_recommendation": "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.", "technical_quality": 3, "clarity_of_presentation": 3, "potential_impact": 2, "robotics_focus": 2, "quality_of_limitations": 2, "questions_for_rebuttal": "1.\tFigure 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. \n2.\tWhat does “reducing the end-to-end cost” mean in the last sentence of the conclusion section?\n3.\tIn 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?\n4.\tIn 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”?\n5.\tI 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. \n6.\tWhat is the connection between “simulators” and “ﬁdelities”?\n7.\tIn line 246, why use $\\gamma=0.43$ as the threshold?\n8.\tIn 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.", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: 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.
Decision: Accept