Bayesian Quadrature for Probability Threshold Robustness of Partially Undeﬁned Functions

In engineering design, one often wishes to calculate the probability that the performance of a system is satisfactory under uncertainty. State of the art algorithms exist to solve this problem using active learning with Gaussian process models. However, these algorithms cannot be applied to problems which often occur in the autonomous vehicle domain where the performance of a system may be undeﬁned under certain circumstances. Naïve modiﬁcation of existing algorithms by simply masking undeﬁned values will introduce a discontinuous system performance function, and would be unsuccessful because these algorithms are known to fail for discontinuous performance functions. We solve this problem using a hierarchical model for the system performance, where undeﬁned performance is classiﬁed before the performance is regressed. This enables active learning Gaussian process methods to be applied to problems where the performance of the system is sometimes undeﬁned, and we demonstrate this by testing our methodology on synthetic numerical examples for the autonomous driving domain. The code to generate these experiments is available as an open source repository: https://github.com/fiveai/hGP_experiments/


Introduction
Estimating the probability that the performance of a system is satisfactory under uncertain or variable operating circumstances is an important step towards deploying such systems safely in the real world. This is especially important in safety critical application such as autonomous driving, where finding rare but catastrophic failures has been identified as a important challenge [1]. Powerful active learning approaches based on Gaussian processes (GP) have been proposed as a solution to this problem [2][3][4][5][6][7] and achieve state of the art performance (we describe related work in detail in Section 2). However, such approaches cannot be applied to problems where the performance of the system may be undefined under certain specific circumstances, a situation which often occurs in the autonomous vehicle domain [1,[8][9][10]. For example, consider the case of an autonomous vehicle waiting to join a road at a T-junction, where other cars are travelling on the road at constant velocity. Assuming that the autonomous vehicle behaves deterministically with respect to variables like the initial starting positions and velocities of the vehicles, one could specify a measure of safe system performance, for example the distance of closest approach between the autonomous vehicle and the other vehicles as a function of these variables. However, in some cases the closest distance function could be undefined because the autonomous vehicle never joined the road at all, for example if the autonomous vehicle's planner deems the situation to be unsafe. Naïvely masking away regions where the performance of the system is undefined would introduce discontinuities and leads to poor performance since Gaussian processes with stationary kernels are not well suited to the regression of discontinuous targets [11]. In this work, we extend these methods from first principles to the case where the performance function can be undefined by using a hierarchical model (termed hGP) for the system performance, where undefined performance is classified before the performance is regressed.

Related Work
The PTR measure has recently become of interest in the robust optimisation literature [12], for example Inatsu et al. [4] show how system designs can be adjusted to optimise the measure. The measure has a much longer history in reliability engineering [14]. Moustapha et al. [6] and Teixeira et al. [15] provide reviews of active learning methods for estimating this measure. The Adaptive Kriging Monte Carlo methodology is perhaps the most well known of these methods, and achieves close to state of the art results [2,3,6]. Efficient methods of estimating the PTR measure also exist in reinforcement learning [13]. Beglerovic et al. [16] use a Bayesian optimisation approach to identify failure cases for an autonomous vehicle but do not exhaustively search for all x such that f (x) < 0 and also do not calculate the PTR measure, p f . Similarly, Wang et al. [17] use a realistic LiDAR simulator to modify real-world LiDAR data which can then be used to test end-to-end autonomous driving systems while searching for adversarial traffic scenarios with Bayesian Optimisation. A related problem in the autonomous vehicle space is finding the most likely x, i.e. with largest p(x), leading to f (x) < 0 [18]. This is closely related to first order methods for estimating the PTR measure [14].
Although there exists literature related to Gaussian Process modelling for discontinuous targets [11], there is little literature on active learning specifically for the PTR measure for discontinuous targets.

Approach
We modify the AK-MCS active learning algorithm by using a different Gaussian Process and acquisition function to Echard et al. [2]. We use a hierarchical Gaussian process model for the rule function, consisting of separate regression and classification Gaussian processes, and a modified acquisition function which minimises the catastrophic event classification error to yield an optimal surrogate model of the rule function. Otherwise our proposed algorithm proceeds in the same way as the AK-MCS algorithm, i.e. an initial training set is chosen to train a Gaussian process, and then subsequent evaluations of the performance function, g, are chosen iteratively by maximising a function of the Gaussian process known as the acquisition function, which are then used to retrain the Gaussian process. The algorithm terminates when the coefficient of variation (CoV) of the failure probability computed using the Gaussian process is below some threshold, and the predicted misclassification probability is also below some threshold. This algorithm is shown in Algorithm 1.
Let y * be the predicted performance for the test input x * ∈ X where y ∈ R ∪ {NaN}, and let the dataset of training examples D = {(x i , y i )|i = 1, ..., n}. We model the predictive distribution p(y * |x * , D) hierarchically as where p(y * |x * , D, y * = NaN) is the predicted regression distribution for y * at the test input x * given that y * is defined, and p(y * = NaN|x * , D) is the predicted probability that y * is undefined for x * . We model these distributions with separate GPs; for p(y * |x * , D, y * = NaN) GP regression is used, and for p(y * = NaN|x * , D) GP classification is used. The conditional prediction of the failure event can easily be calculated as p(y * < 0, y * = NaN|x * , D), which can be used to define an acquisition function, p misclassification (x), based on probability of misclassification of y * < 0 ∩ y * = NaN, as in Echard et al. [2]. We give more details about our modelling approach in Appendix B. Finally, our hierarchical model Eq.
(2) can be used to compute the failure probability as The termination criteria for the AK-MCS algorithm will determine the error in the failure probability computed using the Gaussian process in Eq. (3), in addition to bounding the error of the hierarchical Gaussian process model. This ensures that the model is sufficiently accurate to be used by engineers to make predictions about the behaviour of the system.
Note that in this paper we only consider systems with a single performance function, however Yun et al. [3] demonstrate how acquisition functions for multiple performance functions can be combined when one is interested in a combined PTR measure for the performance functions. This can be applied to our hierarchical model directly. Finally, we note that our modifications to the AK-MCS algorithm are fairly general and only involve changing the performance function model and acquisition function, and therefore these changes could possibly also be used with different active learning algorithms. We do not explore these possible applications in this paper and instead leave this as a topic for future research.
Algorithm 1: Hierarchical Gaussian Process PTR quadrature Input: GP prior GP(0, k), termination threshold η, model g(x) Define proposal set S: sample n mc points from p(x * ). Define initial design of experiments: sample n E points uniformly from S and evaluate with model Choose most likely misclassified x: x * = arg max x∈S p misclassification (x) Observe y * = g(x * ) and Add (x * , y * ) toŜ end while Estimate p f using Monte Carlo simulation with Gaussian Process on S using Eq. (3) Sample n mc points from p(x * ) and evaluate with model g(x) to add to S end while Output: Fitted GP and p f computed using Eq. (3).

Experiments
Benchmark tasks: We evaluate our methodology on two benchmark problems where the system performance is partially undefined and for which p f can be calculated analytically: • Toy function The system performance (plotted in Fig. 2a) is given by where the uncertain variable x is distributed with • Autonomous driving (AD) A mathematical model of an autonomous vehicle joining a main road at a T-junction, as depicted in Fig. 1. The thresholded minimum lateral distance between the autonomous vehicle and an approaching vehicle is given by the performance function (plotted in Fig. 2b) where we have defined the closest approach distance between the vehicles as x a , x lim , d threshold , a ego and v a are physical variables parameterised by distributions shown in Fig. 1. A full derivation of the performance function and explanation of physical variables is provided in Appendix A. Undefined values here represent scenarios when the autonomous vehicle decides not to join the main road, and so the lateral distance along the road is undefined.
Baselines: The hierarchical GP (hGP) will be compared with the following baseline methodologies: • Masked GP: AK-MCS [2] with a regression GP where NaN values are masked with positive • Active GP Classification (GPC): similarly to Kapoor et al. [19], we apply a GP classifier to classify the event y * < 0 ∩ y * = NaN and use this to replace the Kriging model in the AK-MCS active learning loop. Metrics: We use the following metrics, some of which were used in the original AK-MCS paper [2], and some of which we introduce specifically to measure aspects of performance related to our problem.  • Maximum predicted misclassification (max x∈S p misclassification (x)) bounds the "risk" suffered over p(x), and hence can be used to measure the convergence of the internal state of the algorithm, i.e. it does not require ground truth data. Used in Echard et al. [2].
• Failure probability (p f ) can be used to measure if the identified failure region has the correct size, by comparing to the p f computed in some other way as ground truth (i.e. analytically). Used in Echard et al. [2].
• Coefficient of Variation (CoV) is used to assess the algorithms internal uncertainty in the estimated p f , and can be calculated as described in Algorithm 1. Used in Echard et al. [2].
• F 1 score can be used to check if the identified failure region is correctly located in the space of x. We introduce this metric because in reliability problems the failure probability is usually low and hence the class distribution is imbalanced, and F 1 score is known to perform well in such cases. The F 1 score is defined as F 1 = precision×recall precision+recall where precision and recall are measured with a test set of 10 5 points, treating failure as the positive class and the safe region as the negative class. Used in Iwazaki et al. [5].
• Average precision penalises an incorrectly located failure region in a similar way to F 1 score, however it has the additional advantage that the ranking of the predicted class scores is tested. Average precision is insensitive to correctly ranked but miscalibrated scores. It is also measured with a test set of 10 5 points.
Model Accuracy: Firstly we run each algorithm for 150 iterations, i.e. neglect the termination criteria by setting η = 0, in order to study the accuracy of the models independently of the termination criteria. Fig. 3 shows the convergence of the failure probability and maximum predicted misclassification probability for the models for the Toy and Autonomous experiments. In addition, in Appendix D.1 the convergence of the failure probability (p f ) and F 1 score for the models is shown. We observe that the convergence for the hGP is far faster than for the other models. The convergence of p f and predicted misclassification probability for the masked GP is erratic, which is indicative of the misspecification of the model. Visualising the fitted models in Fig. 4, reveals that the low F 1 score and average precision are caused by the length scale for the masked GP becoming extremely short due to the discontinuity in the masked performance function, resulting in a large predicted variance and hence erroneous class scores. hGP outperforming GPC is unsurprising, as the hGP utilises more information about the magnitude of the performance function in order to make a more educated selection of the next point to query.
Overall, although all models eventually accurately estimate p f , the hGP clearly provides a more accurate classification of the failure region, more stable training, and class probabilities which better represent the state of knowledge given the available data.
Termination Criteria: We analyse the ability of the models to terminate the active learning loop after an appropriate number of iterations by repeating the experiments in the previous section with the termination criteria enabled. Table 1 shows the number of samples, the p f , and F 1 score when the algorithm terminates. Appendix D.1 shows additional properties of the models upon termination. The hGP AK-MCS terminates with an accurate model in both experiments. For both experiments we see that GPC eventually learns an accurate model but does not terminate, because the maximum predicted misclassification probability remains high. In the toy function experiment the masked GP terminates far too early while the F 1 score is very low compared to hGP. The masked AK-MCS does not terminate for the AD experiment; the short length scale in the masked GP means that the predicted variance is large and a large amount of iterations will be required to significantly reduce this. For the AD experiment all models obtain the correct value of p f within estimated error, however the F 1 score is much lower for the masked GP, suggesting the identified failure region is incorrectly located.
It is clear that due to the ability of the hGP to predict appropriate classification probabilities it is the only model which can reliably terminate when the model is sufficiently accurate.

Conclusion
In this paper, we tackled the problem of efficiently estimating the performance of systems whose performance function can be undefined under certain circumstances. We presented a novel extension of the AK-MCS algorithm based on modelling the partially defined performance function using a hierarchical Gaussian process, and derived an acquisition function based on the probability of misclassification for undesirable system performance. We showed that our method performs better than existing approaches on two synthetic benchmark systems.
In future work the integration of our hierarchical model into other active learning algorithms could be explored, in addition to testing different termination criteria for the algorithm presented here. In addition, it would be of interest to study the convergence of the algorithm theoretically, and also to extend the proposed algorithm to simultaneously optimise over possible designs of the system itself.

A Mathematical Model
We consider a vehicle moving in the nearside lane at a constant velocity v a , which has starting longitudinal position x(t = 0) = x a < 0 relative to the ego vehicle. The ego vehicle will perceive the adversarial vehicle if the distance between the two vehicles is less than the limiting distance for the sensor (x lim ), otherwise a false negative detection will occur. The scenario will be considered safe provided the distance between the vehicles in the same lane is no less than a threshold d threshold , corresponding to the stopping distance of a typical vehicle. Ego will attempt to merge into the nearside lane and accelerate at its maximum velocity a ego until it reaches v a , however if ego perceives that this action will result in a collision due to the position and velocity of the adversarial vehicle then ego will not merge into the road and hence the longitudinal distance between the vehicles in the lane will be undefined. This is shown in Fig. 1.
Modelling each vehicle as a particle moving in a one dimensional space, we can write the rule numerically as where we have defined the closest approach distance between the vehicles as In order to ensure that the scale of the performance function is appropriate for the GP hyperparameters we have chosen, we rescale the performance function by dividing by 20 resulting in performance values of lower magnitude.

B Hierarchical model details
In this section we provide additional details which describe the model proposed in Section 3 in further detail. Let X be a matrix containing all the x in D and Y be a vector containing all the y in D.
We assume that p(y * |x * , D, y * = NaN) follows a Gaussian process, i.e. p(y * |x * , D, where k : X × X → R is a positive definite kernel with 0 < k(x, x) ≤ 1 for all x ∈ X , and K(X 1 , X 2 ) is matrix containing evaluations of the kernel at all points in X 1 and X 2 , and σ noise is a small positive constant which should be inversely proportional to how deterministic the evaluation of y is. We assume that p(y * = NaN|x * , D) is given by a classification Gaussian process, i.e. p(y * = NaN|x * , D) = probit(h)p(h|x * , D)dh = p nan (x * ), where p(h|x * , D) is the predictive distribution of a Gaussian process which can be calculated using the expectation propagation method as described in Rasmussen and Williams [20].

C Experimental Hyperparameters
For the classification part of the hierarchical GP we use a Matern52 kernel and fix the variance to 10 5 as we find that this ensures a quick convergence in practice, and corresponds to a prior belief that the classification GP should model a deterministic function. For all regression GPs we use a Matern52 kernel with variance and lengthscale determined by optimisation on the training data with some relatively weak constraints: the lengthscale falls in [0, 0.2] and variance falls in [0.5, 1]. For GPC we use a Matern52 kernel with variance 100 and length scale constrained within [0, 0.2]. For all GPs we set the likelihood variance to a small positive number (0.005 2 ), representing a belief that the system performance is deterministic. For AK-MCS we add n mc = 5 × 10 3 proposal samples when the maximum misclassification is below η = 0.02 and terminate the algorithm when the coefficient of variation is below 0.1, which are similar criteria to those used in Echard et al. [2]. In all experiments, the initial design of experiments is 12 samples. The Gaussian Process models were created using GPy [21], and Emukit was used for the active learning algorithms [22].

D Additional Experimental Results
D.1 Additional Figures   Fig. 5 and Fig. 6 show the convergence of the failure probability (p f ), maximum predicted misclassification probability, average precision and F 1 score for the models for the toy function and autonomous driving experiments. In Table 2 we show the number of samples required for each algorithm to terminate and the p f , average precision and F 1 score when the algorithm terminates.

Maximum predicted misclassif cation probability
Hierarchical GP Masked GP α = 1 Masked GP Classif cation Figure 6: Convergence of failure probability (p f ), average precision, maximum predicted misclassification probability and F 1 score for fitted Gaussian processes for the autonomous driving experiment. The shaded area represents the minimum and maximum of 5 repeats, and the dark line represents the mean. All models eventually obtain a p f which is correct within the Monte Carlo error calculated using the CoV (approximately 0.0027).   Fig. 7-8 we repeat the analysis from Section 4 to show the convergence of failure probability (p f ), average precision, maximum predicted misclassification probability and F 1 score for fitted masked Gaussian processes with different α. In comparison to Hierarchical GP, all masked GPs are similar. All models show slower convergence than Hierarchical GP. Table 3 shows the termination results for different values of α.  Figure 7: Convergence of failure probability (p f ), average precision, maximum predicted misclassification probability and F 1 score for fitted masked Gaussian processes with different α for the simple toy function. The shaded area represents the minimum and maximum of 5 repeats, and the dark line represents the mean. All models eventually obtain a p f which is correct within the Monte Carlo error calculated using the CoV (approximately 0.0027).

Maximum predicted misclassification probability
Masked GP α = 1 Masked GP α = 0.5 Masked GP α = 0.1 Figure 8: Convergence of failure probability (p f ), average precision, maximum predicted misclassification probability and F 1 score for fitted masked Gaussian processes with different α for the autonomous driving experiment. The shaded area represents the minimum and maximum of 5 repeats, and the dark line represents the mean. All models eventually obtain a p f which is correct within the Monte Carlo error calculated using the CoV (approximately 0.0027).