Neural Bridge Sampling for Evaluating Safety-Critical Autonomous Systems

Introduction

Data-driven and learning-based approaches have the potential to enable robots and autonomous systems that intelligently interact with unstructured environments. Unfortunately, evaluating the performance of the closed-loop system is challenging, limiting the success of such methods in safety-critical settings. Even if we produce a deep reinforcement learning agent better than a human at driving, flying a plane, or performing surgery, we have no tractable way to certify the system’s quality. Thus, currently deployed safety-critical autonomous systems are limited to structured environments that allow mechanisms such as PID control, simple verifiable protocols, or convex optimization to enable guarantees for properties like stability, consensus, or recursive feasibility (see e.g. ). The stylized settings of these problems and the limited expressivity of guaranteeable properties are barriers to solving unstructured, real-world tasks such as autonomous navigation, locomotion, and manipulation.

The goal of this paper is to efficiently evaluate complex systems that lack safety guarantees and/or operate in unstructured environments. We assume access to a simulator to test the system’s performance. Given a distribution $X\sim P_{0}$ of simulation parameters that describe typical environments for the system under test, our governing problem is to estimate the probability of an adverse event

Problem (1) is often solved in practice by naive Monte Carlo estimation methods, the simplest of which explore the search space via random samples from $P_{0}$. These methods are unbiased and easy to parallelize, but they exhibit poor sample complexity. Naive Monte Carlo can be improved by adding an adaptive component exploiting the most informative portions of random samples drawn from a sequence of approximating distributions $P_{0},P_{1},\dots,P_{K}$. However, standard adaptive Monte Carlo methods (e.g. ), though they may use first-order information on the distributions $P_{k}$ themselves, fail to use first-order information about $f$ to improve sampling; we explicitly leverage this to accelerate convergence of the estimate through optimization.

Naive applications of first-order optimization methods in the estimation problem (1)—for example biasing a sample in the direction $-\nabla f(x)$ to decrease $f(x)$—also require second-order information to correct for the distortion of measure that such transformations induce. Consider the change of variables formula for distributions $\rho(y)=\rho(g^{-1}(y))\cdot|\det J_{g^{-1}}(y)|$ where $y=g(x)$. When $g(x)$ is a function of the gradient $\nabla f(x)$, the volume distortion $|\det J_{g^{-1}}(y)|$ is a function of the Hessian $\nabla^{2}f(x)$. Hessian computation, if even defined, is unacceptably expensive for high-dimensional spaces $\mathcal{X}$ and/or simulations that involve the time-evolution of a dynamical system; our approach avoids any Hessian computation. In contrast, gradients $\nabla f(x)$ can be efficiently computed for many closed-loop systems or through the use of surrogate methods .

To that end, we propose neural bridge sampling, a technique that combines exploration, exploitation, and optimization to efficiently solve the estimation problem (1). Specifically, we consider a novel Markov-chain Monte Carlo (MCMC) scheme that moves along an adaptive ladder of intermediate distributions $P_{k}$ (with corresponding unnormalized densities $\rho_{k}(x)$ and normalizing constants $Z_{k}:=\int_{\mathcal{X}}\rho_{k}(x)dx$). This MCMC scheme iteratively transforms the base distribution $P_{0}$ to the distribution of interest $P_{0}I\{f(x)\leq\gamma\}$. Neural bridge sampling adaptively balances exploration in the search space (via $\nabla\log\rho_{0}$) against optimization (via $\nabla f$), while avoiding Hessian computations. Our final estimate $\hat{p}_{\gamma}$ is a function of the ratios $Z_{k}/Z_{k-1}$ of the intermediate distributions $P_{k}$, the so-called “bridges” . We accurately estimate these ratios by warping the space between the distributions $P_{k}$ using neural density estimation.

Section 2 presents our method, while Section 3 provides guarantees for its statistical performance and overall efficiency. A major focus of this work is empirical, and accordingly, Section 4 empirically demonstrates the superiority of neural bridge sampling over competing techniques in a variety of applications: (i) we evaluate the sensitivity of a formally-verified system to domain shift, (ii) we consider design optimization for high-precision rockets, and (iii) we perform model comparisons for two learning-based approaches to autonomous navigation.

1 Related Work

Several communities have attempted to evaluate the closed-loop performance of cyber-physical, robotic, and embodied agents both with and without learning-based components. Existing solutions are predicated on the definition of the evaluation problem: verification, falsification, or estimation. In this paper we consider a method that utilizes interactions with a gradient oracle in order to solve the estimation problem (1). In contrast to our approach, the verification community has developed tools (e.g. ) to investigate whether any adverse or unsafe executions of the system exist. Such methods can certify that failures are impossible, but they require that the model is written in a formal language (a barrier for realistic systems), and they require whitebox access to this formal model. Falsification approaches (e.g. ) attempt to find any failure cases for the system (but not the overall probability of failure). Similar to our approach, some falsification approaches (e.g. ) utilize gradient information, but their goal is to simply minimize $f(x)$ rather than solve problem (1). Adversarial machine learning is closely related to falsification; the key difference is the domain over which the search for falsifying evidence is conducted. Adversarial examples (e.g. ) are typically restricted to a $p$-norm ball around a point from a dataset, whereas falsification considers all possible in-distribution examples. Both verification and falsification methods provide less information about the system under test than estimation-based methods: they return only whether or not the system satisfies a specification. When the system operates in an unstructured environment (e.g. driving in an urban setting), the mere existence of failures is trivial to demonstrate . Several authors (e.g. ) have proposed that it is more important in such settings to understand the overall frequency of failures as well as the relative likelihoods of different failure modes, motivating our approach.

When sampling rare events and estimating their probability, there are two main branches of related work: parametric adaptive importance sampling (AIS) and nonparametric sequential Monte Carlo (SMC) techniques . Both of these literatures are advanced forms of variance reduction techniques, and they are complementary to standard methods such as control variates . Parametric AIS techniques, such as the cross-entropy method , postulate a family of distributions for the optimal importance-sampling distribution. They iteratively perform heuristic optimization procedures to update the sampling distribution. SMC techniques perform sampling from a sequence of probability distributions defined nonparametrically by the samples themselves. The SMC formalism encompasses particle filters, birth-death processes, and smoothing filters . Our technique blends aspects of both of these communities: we include parametric warping distributions in the form of normalizing flows within the SMC setting.

Our method employs bridge sampling , which is closely related to other SMC techniques such as umbrella sampling , multilevel splitting , and path sampling . The operational difference between these methods is in the form of the intermediate distribution used to calculate the ratio of normalizing constants. Namely, the optimal umbrella sampling distribution is more brittle than that of bridge sampling . Multilevel splitting employs hard barriers through indicator functions, whereas our approach relaxes these hard barriers with smoother exponential barriers. Path sampling generalizes bridge sampling by taking discrete bridges to a continuous limit; this approach is difficult to implement in an adaptive fashion.

The accuracy of bridge sampling depends on the overlap between intermediate distributions $P_{k}$. Simply increasing the number of intermediate distributions is inefficient, because it requires running more simulations. Instead, we employ a technique known as warping, where we map intermediate distributions to a common reference distribution . Specifically, we use normalizing flows , which efficiently transform arbitrary distributions to standard Gaussians through a series of deterministic, invertible functions. Normalizing flows are typically used for probabilistic modeling, variational inference, and representation learning. Recently, Hoffman et al. explored the benefits of using normalizing flows for reparametrizing distributions within MCMC; our warping technique encompasses this benefit and extends it to the SMC setting.

This paper assumes that the generative model $P_{0}$ of the operating domain is given, so all failures are in the modeled domain by definition. When deploying systems in the real world, anomaly detection can discover distribution shifts and is complementary to our approach (see e.g. ). Alternatively, the problem of distribution shift can be addressed offline via distributional robustness , where we analyze the worst-case probability of failure under an uncertainty set composed of perturbations to $P_{0}$.

Proposed approach

As we note in Section 1, naive Monte Carlo measures probabilities of rare events inefficiently. Instead, we consider a sequential Monte Carlo (SMC) approach: we decompose the rare-event probability $p_{\gamma}$ into a chain of intermediate quantities, each of which is tractable to compute with standard Monte Carlo methods. Specifically, consider $K$ distributions $P_{k}$ with corresponding (unnormalized) probability densities $\rho_{k}$ and normalizing constants $Z_{k}:=\int_{\mathcal{X}}\rho_{k}(x)dx$. Let $\rho_{0}$ correspond to the density for $P_{0}$ and $\rho_{\infty}(x):=\rho_{0}(x)I\{f(x)\leq\gamma\}$ be the (unnormalized) conditional density for the region of interest. Then, we consider the following decomposition:

Although we are free to choose the intermediate distributions arbitrarily, we will show below that our estimate for each ratio $Z_{k}/Z_{k-1}$ and thus $p_{\gamma}$ is accurate insofar as the distributions sufficiently overlap (a concept we make rigorous in Section 3). Thus, the intermediate distributions act as bridges that iteratively steer samples from $P_{0}$ towards $P_{K}$. One special case is the multilevel splitting approach , where $\rho_{k}(x):=\rho_{0}(x)I\{f(x)\leq L_{k}\}$ for levels $\infty=:L_{0}>L_{1}\ldots>L_{K}:=\gamma$. In this paper, we introduce an exponential tilting barrier

which allows us to take advantage of gradients $\nabla f(x)$. Here we use the “negative ReLU” function defined as $[x]_{-}:=-[-x]_{+}=xI\{x<0\}$, and we assume that the measure of non-differentiable points, e.g. where $\nabla f(x)$ does not exist or $f(x)=\gamma$, is zero (see Appendix A for a detailed discussion of this assumption). We set $\beta_{0}:=0$ and adaptively choose $\beta_{k}>\beta_{k-1}$. The parameter $\beta_{k}$ tilts the distribution towards the distribution of interest: $\rho_{k}\to\rho_{\infty}$ as $\beta_{k}\to\infty$. In what follows, we describe an MCMC method that combines exploration, exploitation, and optimization to draw samples $X_{i}^{k}\sim P_{k}$. We then show how to compute the ratios $Z_{k}/Z_{k-1}$ given samples from both $P_{k-1}$ and $P_{k}$. Finally, we describe an adaptive way to choose the intermediate distributions $P_{k}$. Algorithm 1 summarizes the overall approach.

Gradient-based MCMC techniques such as the Metropolis-adjusted Langevin algorithm (MALA) or Hamiltonian Monte Carlo (HMC) use gradients $\nabla\log\rho_{0}(x)$ to efficiently explore the space $\mathcal{X}$ and avoid inefficient random-walk behavior . Classical mechanics inspires the HMC approach: HMC introduces an auxiliary random momentum variable $v\in\mathcal{V}$ and generates proposals by performing Hamiltonian dynamics in the augmented state-space $\mathcal{X}\times\mathcal{V}$. These dynamics conserve volume in the augmented state-space, even when performed with discrete time steps .

By including the barrier $\exp\left(\beta_{k}\left[\gamma-f(x)\right]_{-}\right)$, we combine exploration with optimization; the magnitude of $\beta_{k}$ in the barrier modulates the importance of $\nabla f$ (optimization) over $\nabla\log\rho_{0}$ (exploration), two elements of the HMC proposal (see Appendix A for details). We discuss the adaptive choice for $\beta_{k}$ below. Most importantly, we avoid any need for Hessian computation because the dynamics conserve volume. As Algorithm 1 shows, we perform MCMC as follows: given $N$ samples $x_{i}^{k-1}\sim P_{k-1}$ and a threshold $\beta_{k}$, we first resample using their importance weights (exploiting the performance of samples that have lower function value than others) and then perform $T$ HMC steps. In this paper, we implement split HMC which is convenient for dealing with the decomposition of $\log\rho_{k}(x)$ into $\log\rho_{0}(x)+\beta_{k}[\gamma-f(x)]_{-}$ (see Appendix A for details).

Bridge sampling allows estimating the ratio of normalizing constants of two distributions by rewriting

where $\rho_{k}^{B}$ is the density for a bridge distribution between $P_{k-1}$ and $P_{k}$, and $Z_{k}^{B}$ is its associated normalizing constant. We employ the geometric bridge $\rho_{k}^{B}(x):=\sqrt{\rho_{k-1}(x)\rho_{k}(x)}$. In addition to being simple to compute, bridge sampling with a geometric bridge enjoys the asymptotic performance guarantee that the relative mean-square error scales inversely with the Bhattacharyya coefficient, $G(P_{k-1},P_{k})=\int_{\mathcal{X}}\sqrt{\frac{\rho_{k-1}(x)}{Z_{k-1}}\frac{\rho_{k}(x)}{Z_{k}}}dx\in$ (see Appendix B for a proof). This value is closely related to the Hellinger distance, $H(P_{k-1},P_{k})=\sqrt{2-2G(P_{k-1},P_{k})}$. In Section 3, we analyze the ramifications of this fact on the overall convergence of our method.

Both HMC and bridge sampling benefit from warping samples $x_{i}$ into a different space. As Betancourt notes, HMC mixes poorly in spaces with ill-conditioned geometries. Girolami and Calderhead and Hoffman et al. explore techniques to improve mixing efficiency by minimizing shear in the corresponding Hamiltonian dynamics. One way to do so is to transform to a space that resembles a standard isotropic Gaussian .

Conveniently, transforming $P_{k}$ to a common distribution (e.g. a standard Gaussian) also benefits the bridge-sampling estimator (4). As noted above, the error of the bridge estimator grows with the Hellinger distance between the distributions $H(P_{k-1},P_{k})$. However, normalizing constants $Z_{k}$ are invariant to (invertible) transformations. Thus, transformations that warp the space between distributions reduce the error of the bridge-sampling estimator (4). Concretely, we consider invertible transformations $W_{k}$ such that $y^{k}_{i}=W_{k}(x_{i}^{k})$. For clarity of notation, we write probability densities over the space $\mathcal{Y}$ as $\phi$, the corresponding distributions for $Y^{k}$ as $Q_{k}$, and the inverse transformations $W_{k}^{-1}(y)$ as $V_{k}(y)$. Then we can write the bridge-sampling estimate (4) in terms of the transformed variables $y$. The numerator and denominator are as follows:

By transforming all $P_{k}$ into $Q_{k}$ to resemble standard Gaussians, we reduce the Hellinger distance $H(Q_{k-1},Q_{k})\leq H(P_{k-1},P_{k})$. Note that the volume distortions in the expression (5) are functions of the transformation $V_{k}$, so they do not require computation of the Hessian $\nabla^{2}f$. However, computing $\rho_{k}(V_{k}(y))$ requires evaluations of $f$ (e.g. calls of the simulator). We consider the cost-benefit analysis of warping in Section 3.

The KL divergence is an upper bound to the Hellinger distance; we found minimizing the former to be more stable than minimizing the latter. Furthermore, to improve training efficiency, we exploit the iterated nature of the problem and warm-start the weights $\theta_{k}$ with the trained values $\theta_{k-1}$ when solving problem (6) via stochastic gradient descent (SGD). As a side benefit, the trained flows can be repurposed as importance-samplers for the ladder of distributions from nominal behavior to failure.

Because we assume no prior knowledge of the system under test, we exploit previous progress to choose the intermediate $\beta_{k}$ online; this is a key difference to our approach compared to other forms of sequential Monte Carlo (e.g. ) which require a predetermined schedule for $\beta_{k}$. We define the quantities

The first is the fraction of samples that have achieved the threshold. The second is an importance-sampling estimate of $E_{k+1}$ given samples $x_{i}^{k}\sim P_{k}$, written as a function of $\beta$. For fixed fractions $\alpha,s\in(0,1)$ with $\alpha<s$, $\beta_{k+1}$ solves the following optimization problem:

Since $b_{k}(\beta)$ is monotonically decreasing and $b_{k}(\beta)\geq a_{k}$, this problem can be solved efficiently via binary search. The constant $\alpha$ tunes how quickly we enter the tails of $P_{0}$ (smaller $\alpha$ means fewer iterations), whereas $s$ is a stop condition for the last iteration. Choosing $\beta_{k+1}$ via (8) yields a crude estimate for the ratio $Z_{k+1}/Z_{k}$ as $\alpha$ (or $a_{K-1}/s$ for the last iteration). The bridge-sampling estimate $\widehat{E}_{k+1}$ corrects this crude estimate once we have samples from the next distribution $P_{k+1}$.

Performance analysis

We can write the empirical estimator of the function (2) as

where $\widehat{E}_{k}$ is given by the expression (4) without warping, or similarly, as a Monte Carlo estimate of the expression (5) with warping. We provide guarantees for both the time complexity of running Algorithm 1 (i.e. the iterations $K$) as well as the overall mean-square error of $\hat{p}_{\gamma}$. For simplicity, we provide results for the asymptotic (large $N$) and well-mixed MCMC (large $T$) limits. Assuming these conditions, we have the following:

See Appendix B for the proof. We provide some remarks about the above result. Intuitively, the first term in the bound (10) accounts for the variance of $\widehat{E}_{k}$. The denominator of $\widehat{E}_{k-1}$ and numerator of $\widehat{E}_{k}$ both depend on $x_{i}^{k}$; the second sum in (10) accounts for the covariance between those terms. Furthermore, the quantities in the bound (10) are all empirically estimable, so we can compute the mean-square error from a single pass of Algorithm 1. In particular,

Our method can exploit two further sources of efficiency. First, we can employ surrogate models for gradient computation and/or function evaluation during the $T$ MCMC steps. For example, using a surrogate model for a fraction $d\leq 1-1/T$ of the MCMC iterations reduces the factor $T$ to $T_{s}:=(1-d)T$ in the overall cost. Surrogate models have an added benefit of making our approach amenable for simulators that do not provide gradients. The second source of efficiency is parallel computation. Given $C$ processors, the factor $N$ in the cost drops to $N_{c}:=\left\lceil{N/C}\right\rceil$.

The overall efficiency of the estimator (9)—relative error multiplied by cost —depends on $p_{\gamma}$ as $\log(p_{\gamma})^{2}$. In contrast, the standard Monte Carlo estimator has cost $N$ to produce an estimate with relative error $\frac{1-p_{\gamma}}{p_{\gamma}N}$. Thus, the relative efficiency gain for our estimator (9) over naive Monte Carlo is $O(1/({p_{\gamma}\log(p_{\gamma})^{2}}))$: the efficiency gains over naive Monte Carlo increase as $p_{\gamma}$ decreases.

Experiments

We evaluate our approach in a variety of scenarios, showcasing its use in efficiently evaluating the safety of autonomous systems. We begin with a synthetic problem to illustrate the methodology concretely as well as highlight the pitfalls of using gradients naively. Then, we evaluate a formally-verified neural network controller on the OpenAI Gym continuous MountainCar environment under a domain perturbation. Finally, we consider two examples of using neural bridge sampling as a tool for engineering design in high-dimensional settings: (a) comparing thruster sizes to safely land a rocket in the presence of wind, and (b) comparing two algorithms on the OpenAI Gym CarRacing environment (which requires a surrogate model for gradients) .

Figure 2(b) shows contours of $f(x)$. Notably, the failure region (dark blue) is an extremely irregular geometry with pathological curvature, which renders MCMC difficult for AMS and B . Quantitatively, poor mixing adversely affects the performance of AMS and B, and they perform even worse than MC (Figure 2(c)). Whereas gradients help B slightly over AMS, gradients and neural warping together help NB outperform all other methods. We next move to higher-dimensional systems.

We now consider the problem of autonomous, high-precision vertical landing of an orbital-class rocket (Figure 3(a)), a technology first demonstrated by SpaceX in 2015. Rigorous system-evaluation techniques such as our risk-based framework are powerful tools for quickly exploring design tradeoffs. In this experiment, the amount of thrust which the rocket is capable of deploying to land safely must be balanced against the payload it is able to carry to space; stronger thrust increases safety but decreases payloads. We consider two rocket designs and we evaluate their respective probabilities of failure (not landing safely on the landing pad) for landing pad sizes up to $15$ meters in radius. That is, $-f(x)$ is the distance from the landing pad’s center at touchdown and $\gamma=-15$. We evaluate whether the rockets perform better than a threshold failure rate of $10^{-5}$.

We let $P_{0}$ be the 100-dimensional search space parametrizing the sequence of wind-gusts during the rocket’s flight. Appendix C contains details for this parametrization and the closed-loop simulation of the rocket’s control law (based on industry-standard approaches ). Figure 3(b) shows the estimated performance of the two rockets. We show only MC and NB for clarity; comparisons with other methods are in Table LABEL:tab:results (with ground-truth values calculated using 50 million naive Monte Carlo simulations). Whereas both NB and MC confidently estimate Rocket2’s failure rate as higher than $10^{-4}$, only NB confidently estimates Rocket1’s failure rate as higher than $10^{-5}$, letting engineers quickly judge whether to increase the size of the landing pad or build a better rocket.

We can also distinguish between the modes of failure for the rockets. Namely, Figure 3(c) shows a PCA projection of failures (with $\gamma_{\rm test}=-15$) onto 2 dimensions. Analysis of the PCA modes indicates that failures are dominated by high altitude and medium altitude gusts. Even though Rocket2 has a higher probability of failure, its failure mode is more concentrated than Rocket1’s failures.

The CarRacing environment (Figure 4(a)) is a challenging reinforcement-learning task with a continuous action space and pixel observations. Similar observation spaces have been proposed for real autonomous vehicles (e.g. ). We compare two recent approaches, AttentionAgentRacer and WorldModelRacer that have similar average performance: they achieve average rewards of $903\pm 49$ and $899\pm 46$ respectively (mean $\pm$ standard deviation over 2 million trials). Both systems utilize one or more deep neural networks to plan in image-space, so neither has performance guarantees. We evaluate the probability of getting small rewards ($\gamma=150$).

The 24-dimensional search space $P_{0}$ parametrizes the generation of the racing track (details are in Appendix C). This environment does not easily provide gradients due to presence of a rendering engine in the simulation loop. Instead, we fit a Gaussian process surrogate model to compute $\nabla f(x)$ (see Appendix C). As these experiments are extremely expensive (taking up to 1 minute per simulation), we only use 2 million naive Monte Carlo samples to compute the ground-truth failure rates. Figure 4(b) shows that, even though the two models have very similar average performance, their catastrophic failure curves are distinct. Furthermore, MC is unable to distinguish between the policies below rewards of 160 due to its high uncertainty, whereas NB clearly shows that WorldModelRacer is superior. Note that, because even the ground-truth has non-negligible uncertainty with 2 million samples, we only report the variance component of relative mean-square error in Table LABEL:tab:results.

As with the rocket design experiments, we visualize the modes of failure (defined by $\gamma_{\rm test}=225$) via PCA in Figure 4(c). The dominant eigenvectors involve large differentials between radii and angles of consecutive checkpoints that are used to generate the racing tracks. AttentionAgentRacer has two distinct modes of failure, whereas WorldModelRacer has a single mode.

Conclusion

There is a growing need for rigorous evaluation of safety-critical systems which contain components without formal guarantees (e.g. deep neural networks). Scalably evaluating the safety of such systems in the presence of rare, catastrophic events is a necessary component in enabling the development of trustworthy high-performance systems. Our proposed method, neural bridge sampling, employs three concepts—exploration, exploitation, and optimization—in order to evaluate system safety with provable statistical and computational efficiency. We demonstrate the performance of our method on a variety of reinforcement-learning and robotic systems, highlighting its use as a tool for continuous integration and rapid engineering design. In future work, we intend to investigate how efficiently sampling rare failures—like we propose here for evaluation—could also enable the automated repair of safety-critical reinforcement-learning agents.

Broader Impact

This paper presents both foundational theory and methods for efficiently evaluating the performance of safety-critical autonomous systems. By definition, such systems can cause injury or death if they malfunction . Thus, improving the tools that practitioners have to perform risk-estimation has the potential to provide a strong positive impact. On the other hand, the improved scalability of our method could be used to more efficiently find (zero-day) exploits and failure modes in $P_{0}$ (the model of the operational design domain). However, we note that adversarial examples or exploits can also be found via a variety of purely optimization-based methods . The nuances of our method are primarily concerned with the frequency of adverse events, an extra burden; thus, we anticipate they will be of little interest to malicious actors who can manipulate the observations and sensor measurements of complex systems. Another potential concern about the use of our method is with respect to the identification of $P_{0}$, which we specifically assume to be known in this paper. The gap between $P_{0}$ in simulation and the real distribution of the environment could lead to overconfidence in the capabilities of the system under test. In Section 1.1 we outline complementary work in anomaly detection and distributionally robust optimization which could mitigate such risks. Still, more work needs to be done to standardize the operational domain of specific tasks by regulators and technology-stakeholders. Nevertheless, we believe that our method will enable the comparison of autonomous systems in a common language—risk—across the spectrum from engineers to regulators and the public.

The applications of our technology are diverse (cf. Corso et al. ), ranging from testing autonomous vehicles and medical devices to evaluating deep neural networks and reinforcement-learning agents . In the case of autonomous vehicles, Sparrow and Howard argue that it will be morally wrong not to deploy self-driving technology once performance exceeds human capabilities. Our work is an important tool for determining when this performance threshold is achieved due to the rare nature of serious accidents . While the widespread availability of autonomy-enabled devices could narrowly benefit public health, there are many external risks associated with their development. First, many learning-based components of these systems will require massive and potentially invasive data collection ; preserving privacy of the public via federated learning and differential privacy-based mechanisms should remain important initiatives within the machine-learning community. A second potential negative consequence of the applications like autonomous vehicles is the use of the real-world as a “simulator” within a reinforcement-learning scheme by releasing “beta” autonomy features (e.g. Tesla Autopilot ). Unlike established industries such as aerospace , many potential applications currently lack regulation and standards; it is important to ensure that industry works with policy makers to develop safety standards in a way that avoids regulatory capture. If widely adopted in regulatory frameworks, our tool would enable rational decisions about the impact, positive or negative, of safety-critical autonomous systems before real lives are affected.

More broadly, the advent of autonomy could spark significant societal changes. For example, the autonomous applications described previously could become core components of weapons systems and military technology that are incompatible with (modern interpretations of) just war theory . Similarly, the automation of the transportation industry has the potential to rapidly destroy the economics of public infrastructure and cost millions of jobs . Thus, Benkler highlights that there is a growing need for the academic community to take action on defining the broader performance criteria to which we will hold AI applications. Brundage et al. and Wing outline broad research agendas which are necessarily interdisciplinary. Still, much more work needs to be done to empower researchers to influence policy. These efforts will require systemic initiatives by research institutions and organizations to engage with local, national, and international governing bodies.

Acknowledgements

AS and JD were partially supported by the DAWN Consortium, NSF CAREER CCF-1553086, NSF HDR 1934578 (Stanford Data Science Collaboratory), ONR YIP N00014-19-2288, and the Sloan Foundation. MOK was supported by an NSF GRFP Fellowship. RT was supported by Lincoln Laboratory/Air Force Award No. PO# 7000470769 and Amazon Robotics Award No. CC MISC 00272683 2020 TR.

References

Appendix A Warped Hamiltonian Monte Carlo (HMC)

In this section, we provide a brief overview of HMC as well as the specific rendition, split HMC . Given “position” variables $x$ and “momentum” variables $v$, we define the Hamiltonian for a dynamical system as $H(x,v)$ which can usually be written as $U(x)+K(v)$, where $U(x)$ is the potential energy and $K(v)$ is the kinetic energy. For MCMC applications, $U(x)=-\log(\rho_{0}(x))$ and we take $v\sim\mathcal{N}(0,I)$ so that $K(v)=\|v\|^{2}/2$. In HMC, we start at state $x_{i}$ and sample $v_{i}\sim\mathcal{N}(0,I)$. We then simulate the Hamiltonian, which is given by the partial differential equations:

Of course, this must be done in discrete time for most Hamiltonians that are not perfectly integrable. One notable exception is when $x$ is Gaussian, in which case the dynamical system corresponds to the evolution of a simple harmonic oscillator (i.e. a spring-mass system). When done in discrete time, a symplectic integrator must be used to ensure high accuracy. After performing some discrete steps of the system (resulting in the state $(x_{f},v_{f})$), we negate the resulting momentum (to make the resulting proposal reversible), and then accept the state $(x_{f},-v_{f})$ using the standard Metropolis-Hastings criterion: $\min(1,\exp(-H(x_{f},-v_{f})+H(x_{i},v_{i})))$ .

The standard symplectic integrator—the leap-frog integrator—can be derived using the following symmetric decomposition of the Hamiltonian (performing a symmetric decomposition retains the reversibility of the dynamics): $H(x,v)=U(x)/2+K(v)+U(x)/2$. Using simple Euler integration for each term individually results in the following leap-frog step of step-size $\epsilon$:

where each step simply simulates the individual Hamiltonian $H_{1}(x,v)=U(x)/2$, $H_{2}(x,v)=K(v)$, or $H_{3}(x,v)=U(x)/2$ in sequence. As presented by Shahbaba et al. , this same decomposition can be done in the presence of more complicated Hamiltonians. In particular, consider the Hamiltonian $H(x,v)=U_{1}(x)+U_{0}(x)+K(v)$. We can decompose this in the following manner: $H_{1}(x,v)=U_{1}(x)/2$, $H_{2}(x,v)=U_{0}(x)+K(v)$, and $H_{3}(x,v)=U_{1}(x)/2$. We can apply Euler integration to the momentum $v$ for the first and third Hamiltonians and the standard leap-frog step to the second Hamiltonian (or even analytic integration if possible). For this paper, we have $U_{0}(x)=-\log\rho_{0}(x)$ and $U_{1}(x)=-\beta[\gamma-f(x)]_{-}$.

To account for warping, the modifications needed to the HMC steps above are simple. When performing warping, we simply perform HMC for a Hamiltonian $\hat{H}(y,v)$ that is defined with respect to the warped position variable $y$, where $x=V_{\theta}(y)$ for given parameters $\theta$. By construction of the normalizing flows, we assume $y\sim\mathcal{N}(0,I)$, so that we can perform the dynamics for $\hat{H}_{2}(y,v)$ analytically. Furthermore, the Jacobian $J_{V_{\theta}}(y)$ is necessary for performing the Euler integration of $H_{1}(y,v)$ and $H_{3}(y,v)$. This is summarized in Algorithm 2. Note that we always perform the Metropolis-Hastings acceptance with respect to the true Hamiltonian $H$, rather than the Hamiltonian $\hat{H}$ that assumes perfect training of the normalizing flows.

In Section 2, we assumed that the measure of non-differentiable points is zero for the energy potentials considered by HMC. As discussed by Afshar and Domke , the inclusion of the Metropolis-Hastings acceptance criterion as well as the above assumption ensures that HMC asymptotically samples from the correct distribution even for non-smooth potentials. An equivalent intuitive explanation for this can be seen by viewing the ReLU function $[x]_{+}$ as the limit of softplus functions $g_{k}(x):=\log(1+\exp(kx))/k$ as the sharpness parameter $k\to\infty$. We can freely choose $k$ such that, up to numerical precision, Algorithm 2 is the same whether we consider using a ReLU or sufficiently sharp (e.g. large $k$) softplus potential, because, with probability one, we will not encounter the points where the potentials differ. When further knowledge about the structure of the non-differentiability is known, the acceptance rate of HMC proposals can be improved .

Appendix B Performance analysis

We begin with showing the convergence of the number of iterations. To do this, we first show almost sure convergence of $\beta_{k}$ in the limit $N\to\infty$. We note that in the optimization problem (8), $\beta_{k}$ is a feasible point, yielding $b_{k}(\beta)=1$. Thus, $\beta_{k+1}\geq\beta_{k}\geq\beta_{0}:=0$. Due to this growth of $\beta_{k}$ with $k$, we have

Now, we consider the convergence of the solutions to the finite $N$ versions of problem (8), denoted $\beta_{k}^{N}$, to the “true” optimizers $\beta_{k}$ in the limit as $N\to\infty$. Leaving the dependence on $\beta_{k}$ implicit for the moment, we consider the random variable $Y:=g(X;\beta):=\exp\left((\beta-\beta_{k})[\gamma-f(X)]_{-}\right)$. Then, since $\beta\in\mathcal{B}$ is bounded and $g$ is continuous in $\beta$, we can state the Glivenko-Cantelli convergence of the empirical measure uniformly over $\mathcal{B}$: $\sup_{\beta\in\mathcal{B}}\|F^{N}(Y)-F(Y)\|_{\infty}\to 0$ almost surely, where $F$ is the cumulative distribution function for $Y$. Note that the constraints in the problem (8) can be rewritten as expectations of this random variable $Y$. Furthermore, the function $g$ is strictly monotonic in $\beta$ (and therefore invertible) for non-degenerate $f(X)$ (i.e. $f(x)>\gamma$ for some non-negligible measure under $P_{0}$). Thus, we have almost sure convergence of the argmin $\beta_{k+1}^{N}$ to $\beta_{k+1}$.

Until now, we have taken dependence on $\beta_{k}$ implicitly. Now we make the dependence explicit to show the final step of convergence. In particular, we can write $\beta_{k+1}$ as a function of $\beta_{k}$ (along with their empirical counterparts), For concreteness, we consider the following decomposition for two iterations:

We have already shown above that the first term on the right hand side vanishes almost surely. By the same reasoning, we know that $\beta_{1}^{N}\to\beta_{1}$ almost surely. The second term also vanishes almost surely since $\beta_{k+1}(\beta)$ is a continuous mapping. This is due to the fact that the constraint functions in problem (8) are continuous functions of both $\beta$ and $\beta_{k}$ along with the invertibility properties discussed previously. Then, we simply extend the telescoping series above for any $k$ and similarly show that all terms vanish almost surely. This shows the almost sure convergence for all $\beta_{k}$ up to some $K$.

Now we move to the relative mean-square error of $\hat{p}_{\gamma}$. We employ the delta method, whereby, for large $N$, this is equivalent to $\text{Var}(\log(\hat{p}_{\gamma}))$ (up to terms $o(1/N)$). For notational convenience, we decompose $\widehat{E}_{k}$ into its numerator and denominator:

By construction (and assumption of large $T$), Algorithm 1 has a Markov property that each iteration’s samples $x_{i}^{k}$ are independent of the previous iterations’ samples $x_{i}^{k-1}$ given $\beta_{k}$. For shorthand, let $\beta_{0:k}$ denote all $\beta_{0},\ldots,\beta_{k}$. Conditioning on $\beta_{0:k}$, we have

Since $\beta_{0:k}$ approaches constants almost surely as $N\to\infty$, the first term vanishes and the second term is the expectation of a constant. In particular, the second term is as follows:

Similarly, $\text{Var}(B_{k}|\beta_{0:k})=Z_{k-1}/Z_{k}-(Z_{k}^{B}/Z_{k})^{2}$. Next we look at the covariance terms:

Again, the first term vanishes since $\beta_{0:k}$ approach constants as $N\to\infty$. By construction, the second term is also 0 since the quantities are conditionally independent. Similarly, $\text{Cov}(B_{k-1},B_{k})=0$ and $\text{Cov}(A_{i},B_{j})=0$ for $j\neq i-1$. However, there is a nonzero covariance for the quantities that depend on the same distribution:

By the large $T$ assumption, the samples $x_{i}^{k}$ and $x_{j}^{k}$ are independent for all $i\neq j$ given $\beta_{k}$. Then we have

The last term in $\hat{p}_{\gamma}$, $\frac{1}{N}\sum_{i=1}^{N}\frac{\rho_{\infty}(x^{K}_{i})}{\rho_{K}(x^{K}_{i})}$, reduces to a simple Monte Carlo estimate since $\frac{\rho_{\infty}(X)}{\rho_{K}(X)}=I\{f(X)\leq\gamma\}$. Furthermore, this quantity is independent of all other quantities given $\beta_{0:K}$ and, as noted above, approaches $s$ almost surely as $N\to\infty$.

Putting this all together, the delta method gives (as $N\to\infty$ so that $\beta_{0:K}$ approach constants almost surely),

The Bhattacharrya coefficient can be written as

We remark that a special case of this formula is for $K=1$ and $s=1$ (so only the first term survives), which is the relative mean-square error for a single bridge-sampling estimate $\widehat{E}_{k}$.

Now, since $G(P,Q)\geq 0$, the terms in the second sum are $\geq-1$ so that the second sum is $\leq 2(K-1)/N$. Furthermore, since $s\geq 1/3$, the last term is also $\leq 2/N$. Thus, if we have $\frac{1}{G(P_{k-1},P_{k})^{2}}\leq D$ (with $D\geq 1$), then the asymptotic relative mean-square error (12) is $\leq 2KD/N$ (up to terms $o\left(\frac{1}{N}\right)$).

When performing warping, we follow the exact same pattern as the above results, conditioning on both $\beta_{0:k}$ and $W_{0:k}$, where $W_{0}$ is defined as the identity mapping. We follow the same almost-sure convergence proof for $W_{k}$ as above for $\beta_{k}$, which requires compactness of $\theta\in\Theta$, continuity of $W$ with respect to $\theta$ and $x$, and that we actually achieve the minimum in problem (6). Although the first two conditions are immediate in most applications, the last condition can be difficult to satisfy for deep neural networks due to the nonconvexity of the optimization problem.

Appendix C Experimental setups

The number of samples $N$ affects the absolute performance of all of the methods tested, but not their relative performance with respect to each other. For all experiments, we use $N=1000$ for B and NB to have adequate absolute performance given our computational budget (see below for the computing architecture used). Other hyperparameters were tuned on the synthetic problem and fixed for the rest of the experiments (with the exception of the MAF architecture for the rocket experiments). The hyperparameters were chosen as follows.

When performing Hamiltonian dynamics for a Gaussian variable, a time step of $2\pi$ results in no motion and time step of $\pi$ results in a mode reversal, where both the velocity and position are negated. The $\pi$ time step is in this sense the farthest exploration that can occur in phase space (which can be intuitively understood by recognizing that the phase diagram of a simple spring-mass system is a unit circle). Thus, we considered $T=4,8,12,$ and $16$ with time steps $\pi/T$. We found that $T=8$ provided reasonable exploration (as measured by autocorrelations and by the bias of the final estimator $\hat{p}_{\gamma}$) and higher values of $T$ did not provide much more benefit. For B, we allowed 2 more steps $T=10$ to keep the computational cost the same across B and NB. Similarly, for AMS, we set $T=10$. We also performed tuning online for the time step to keep the accepatance ratio between 0.4 and 0.8. This was done by setting the time step to $\sin^{-1}(\min(1,\sin(t)\exp((p-C)/2))$, where $t$ is the current time step, $p$ is the running acceptance probability for a single chain and $C=0.4$ if $p<0.4$ or $0.8$ if $p>0.8$. This was done after every $T$ HMC steps.

The MAF architectures for the synthetic, MountainCar, and CarRacing experiments were set at 5 MADE units, each with 1 hidden layer of 100 neurons. Because the rocket search space is very high dimensional, we decreased the MAF size for computational efficiency: we set it at 2 MADE units, each with hidden size 400 units. We used 100 epochs for training, a batch size of 100, a learning rate of 0.01 and an exponential learning-rate decay with parameter 0.95.

For the surrogate Gaussian process regression model for CarRacing, we retrained the model on the most recent $N$ simulations after every $NT$ simulations (e.g. after every $T$ HMC iterations). This made the amortized cost of training the surrogate model negligible compared to performing the simulations themselves. We used a Matern kernel with parameter $\nu=2.5$. We optimized the kernel hyperparameters using an L-BFGS quasi-Newton solver.

C.2 Environment details

The MountainCar environment considers a simple car driving on a mountain road. The car can sense horizontal distance $s$ as well as its velocity $v$, and may send control inputs $u$ (the amount of power applied in either the forward or backward direction). The height of the road is given by: $h(s)=0.45\sin(3s)+0.55$. The speed of the car, $v$, is a function of $s$ and $u$ only. Thus, the discrete time dynamics are: $s_{k+1}=s_{k}+v_{k+1}$ and $v_{k+1}=v_{k}+0.0015u_{k}-0.0025\cos(3s_{k})$. For a given episode the agent operating the car receives a reward of $-0.1u_{k}^{2}$ for each control input and $100$ for reaching the goal state.

In this experiment we explore the effect of domain shift on a formally verified neural network. We utilize the neural network designed by Ivanov et al. ; it contains two hidden layers, each of 16 neurons, for a total of 337 parameters. For our experiments we use the trained network parameters available at: https://github.com/Verisig/verisig. Ivanov et al. describe a layer-by-layer approach to verification which over-approximates the reachable set of the combined dynamics of the environment and the neural network. An encoding of this system (network and environment) is developed for the tool Flow∗ which constructs the (overapproximate) reachable set via a Taylor approximation of the combined dynamics.

The MountainCar environment is considered solved if a policy achieves an average reward of $90$ over $100$ trials. The authors instead seek to prove that the policy will achieve a reward of at least $90$ for any initial condition. By overapproximating the reachable states of the system, they show that the car always receives a total reward greater than $90$ and achieves the goal in less than $115$ steps for a subset of the intial conditions $\hat{p}_{0}\in\left[-0.59,-0.4\right]$.

C.2.2 Rocket design

The system under test is a rocket spacecraft with dynamics $m\ddot{p}=f-mge_{3}$ , where $m>0$ is the mass, $p(t)\in\mathbf{R}^{3}$ is the position, and $e_{3}$ is the unit vector in the z-direction. While it is possible to synthesize optimal trajectories for an idealized model of the system, significant factors such as wind and engine performance (best modeled as random variables) are unaccounted for . Without feedback control, even small uncorrected tracking errors result in loss of the vehicle. In the case of disturbances the authors suggest two approaches: (1) a feedback control law which tracks the optimal trajectory (2) receding horizon model predictive control. The system we consider tracks an optimal trajectory using a feedback control law. Namely, the optimal trajectory is given by the minimum fuel solution to a linearized mode of the dynamics. Specifically, we consider the thrust force discretized in time with a zero-order hold, such that $f_{k}$ applied for time $t\in[(k-1)h,kh]$ for a time step $h=0.2$. Then, the reference thrust policy solves the following convex optimization problem

where $B\sim\text{Bernoulli}(1/3)$ and $\hat{B}\sim\text{Bernoulli}(1/2)$. This results in 5 random variables for each second, or a total of 100 random variables since we have a 20 second simulation. The wind intensity experienced by the rocket is a linear function of height (implying a simplistic laminar boundary layer): $f_{w}=CWp_{3}$ for a constant $C$. Finally, the rocket has a proportional feedback control law for the booster thrusters to the errors in both the position $p_{k}$ and velocity $v_{k}$:

The maximum norm for clip-by-norm is $aF_{max}$, where $a=1.15$ for Rocket1 and $a=1.1$ for Rocket2, indicating that the boosters are capable of providing $15\%$ or $10\%$ of the thrust of the main engine.

C.2.3 Car Racing

We compare the failure rate of agents solving the car-racing task utilizing the two distinct approaches ( and ). The car racing task differs from the other experiments due to the inclusion of a (simple) renderer in the system dynamics. At each the step the agent recieves a reward of $-0.1+\mathcal{I}_{newtile}(1000/N)-\mathcal{I}_{offtrack}(100)$ where N is the total number of tiles visited in the track. The environment is considered solved if the agent returns an average reward of $900$ over $100$ trials. The search space $P_{0}$ is the inherent randomness involved with generating a track. The track is generated by selecting 12 checkpoints in polar coordinates, each with radian value uniformly in the interval $[2\pi i/12,2\pi(i+1)/12)$ for $i=0,\ldots 11$, and with radius uniformly in the interval $[R/3,R]$, for a given constant value $R$. This results in 24 parameters in the search space. The policies used for testing are described below (with training scripts in the code supplement).

Tang et al. utilize a simple self-attention module to select patches from a 96x96 pixel observation. First the input image is normalized then a sliding window approach is used to extract $N$ patches of size $M\times M\times 3$ which are flattened and arranged into a matrix of size $3M^{2}\times N$. The self-attention module is used to compute the attention matrix $A$ and importance vector (summation of each column of $A$). A feature extraction operation is applied to the top K elements of the sorted importance vector and the selected features are input to a neural network controller. Both the attention module and the controller are trained together via CMA-ES. Together, the two modules contain approximately 4000 learnable parameters. We use the pre-trained model available here: https://github.com/google/brain-tokyo-workshop/tree/master/AttentionAgent.

The agent of Ha and Schmidhuber first maps a top-down image of the car on track via a variational autoencoder to a latent vector $z$. Given $z$, the world model $M$ utilizes a recurrent-mixture density network to model the distribution of future possible states $P(z_{t+1}\mid a_{t},z_{t},h_{t})$. Note that $h_{t}$, the hidden state of the RNN. Finally, a simple linear controller $C$ maps the concatenation of $z_{t}$ and $h_{t}$ to the action, $a_{t}$. We use the pre-trained model available here: https://github.com/hardmaru/WorldModelsExperiments/tree/master/carracing.