arXiv CS open access

Clustered Randomized Smoothing for Stochastic Prediction Functions

2026 · arxiv_cs

arXiv CS · Papers · 2026

deep-learning machine-learning neural-networks

Conceptio Open Archive

Analyze & explore related research in Conceptio

Connect this record directly into your research workflow or AI agent. Explore citations across open repositories, query full text, and export structured proof bundles.

Explore in Search →Connect via MCP / API ↗

Clustered Randomized Smoothing for Stochastic Prediction Functions Eduardo Figueiredo⋆,1 Frederik Mathiesen⋆,1 Julian Schumann2 Jens Kober3 Arkady Zgonnikov2 Luca Laurenti1,4 1 Delft Center for Systems and Control, Delft University of Technology 2 Cognitive Robotics, Delft University of Technology 3 University of Stuttgart 4 AI4I {e.figueiredo, f.b.mathiesen, j.f.schumann, a.zgonnikov, l.laurenti}@tudelft.nl

arXiv:2608.12037v1 [cs.LG] 12 Aug 2026

Abstract Modern stochastic predictors can model rich, multi-modal outcome distributions. However, this expressive power comes with challenges in ensuring robust predictions—a critical requirement in safety-critical domains. Randomized smoothing is a leading technique for improving robustness, particularly against adversarial perturbations. Yet, in stochastic multi-modal regression settings, randomized smoothing often fails due to mode collapse, yielding averaged predictions that do not reflect the underlying distribution. To address this limitation, we propose clustered α-smoothing, a framework that (1) partitions noisy samples using an arbitrary clustering algorithm, (2) applies α-smoothing locally within each cluster, and (3) combines the resulting predictions into a mixture distribution. By interpreting the smoothing distribution as a mixture of α-smoothers, we derive a lower bound on the probability that the smoothed prediction lies within a union of compact regions corresponding to distinct modes. We empirically evaluate our framework on two benchmarks, demonstrating substantial improvements over state-of-the-art methods. In stochastic trajectory prediction on a driving simulator dataset, our approach achieves, on average, a 27% lower Wasserstein distance to the groundtruth distribution compared to α-smoothing. In quadrotor control, where modes correspond to distinct feasible paths to a target, our method reduces the collision rate by 81% relative to the state-of-the-art randomized smoothing.

Introduction

Stochastic predictors, such as variational autoencoders [10], normalizing flows [28], or Bayesian neural networks [43], are capable of modeling rich, multi-modal outcome distributions [21, 12]. However, this expressive power comes with challenges in ensuring robust predictions – a fundamental requirement in safety-critical applications such as electrical grids, healthcare, robotics, or autonomous driving, where this type of brittle behavior can have catastrophic consequences [9, 24, 23, 37]. A particular concern is vulnerability to adversarial attacks [15, 22, 35, 29, 6], in which carefully crafted, small perturbations of the input cause significant deviations in the network’s output. To address this issue, a wide range of methods has been developed for formally verifying desired input–output properties of neural networks, including mixed-integer programming [16], abstract interpretation [13], linear relaxations with branch-and-bound [48, 45], satisfiability modulo theory solvers [44, 11], and reachability analysis [41, 1]. Despite their strong guarantees, these methods are often computationally expensive, which limits their practical applicability. Randomized smoothing is a scalable alternative for certifying robustness [20, 17, 8]. Given a base predictor, randomized smoothing injects noise into the input and averages the resulting predictions, yielding a smoothed predictor with certified robustness guarantees. While most prior work focuses on classification [20, 17, 8], recent efforts extend these ideas to regression [7, 31, 30, 36]. In this setting, Preprint.

R1 R

R = R1 ∪ R2 R2

(a) base predictor

(b) randomized smoothing

Figure 1: Conceptual illustration of mode collapse in randomized smoothing. In (a), an ego vehicle (red) is predicting the trajectory of another vehicle (green) at a cross junction. The stochastic predictor outputs a distribution with two modes: turning left before the ego vehicle passes the intersection, or waiting for it to pass before initiating the turn. In (b), current randomized smoothing methods (e.g., [30]) would produce an averaged prediction, failing to capture any of the two modes. Instead, in (c), our approach is to smooth each mode individually and combine as a mixture model, which preserves the modes corresponding to the two tactical behaviors. Sets R represent the sets that can be certified with each method.

given a base predictor h : Rd → Rq and a Gaussian noise ε ∼ N (0, σ 2 I) with σ 2 > 0, the smoothed predictor is commonly defined as (some variant of) H(x) = Eε [h(x + ε)]. However, this approach is fundamentally limited for stochastic neural network predictors. In these models, predictions are often inherently multi-modal [10, 28, 43], yet randomized smoothing collapses this structure by averaging over noisy outputs, yielding a unimodal and potentially uninformative estimate (Figure 1). This limitation is particularly problematic in safety-critical applications where multi-modality is essential, such as trajectory prediction [27, 26, 4, 37] or failure prediction [33, 49]. To address this limitation and capture the multi-modal behavior of the base predictor, we propose clustered α-smoothing. Given an input x and a set of noise samples {ε1 , . . . , εN }, our method (i) clusters the prediction samples {hw (x + εi ), . . . , hw (x + εN )}, (ii) trims α outlier samples per cluster, and (iii) smooths locally within each cluster. The clusters are then combined according to a mixture distribution. By combining the modes/clusters as mixture model rather than averaging, we avoid mode collapse, and thus, the limitation of the previous state-of-the-art. Further, the per-mode α-trimming allows greater control over outliers. We prove by classical probability theory results and statistical methods that the output of the smoothed predictor belongs to a coverage region R with high probability for all inputs in a radius r. The guaranteed coverage region R is provided as the union of the cluster-wise coverage regions Rm , which are regions containing the samples of each cluster with high probability (see Figure 1), and the probability of the prediction belonging to the coverage region is provided both individually and jointly for all clusters. Our framework is agnostic to the clustering algorithm and provides flexibility in how to determine the number of clusters and sample assignment. We demonstrate the framework on two benchmarks: (1) stochastic multi-modal trajectory prediction in traffic (Figure 1), a key prerequisite for safe and efficient automated driving, and (2) robustification of a multi-modal deep RL quadrotor controller, where the framework helps improve robustness while preserving multiple feasible navigation strategies around obstacles. In summary, our main contributions are: (i) we introduce a framework for randomized smoothing for multi-modal stochastic predictors via clustering and α-trimming, (ii) we prove that the output of the clustered α-smoothed predictor with high probability falls within a small coverage region per cluster for all inputs in a radius, and (iii) we demonstrate the effectiveness of our framework in benchmarks from stochastic trajectory prediction in traffic and a multi-modal RL controller for a quadrotor, showing improved performance compared to state-of-the-art.

Preliminaries on randomized smoothing

Notation. For a vector x ∈ Rd , we denote by x(i) its i-element. Given N ∈ N>0 , the set {1, . . . , N } is denoted by [N ]. I denotes the identity matrix. The Pontryagin difference between sets A, B is denoted by A ⊖ B. For a set X ⊆ Rd , the indicator function for X is denoted as 1X (x) := 1 if x ∈ X ; otherwise 0. Given a Borel measurable space X ⊆ Rd , we denote by B(X ) the Borel sigma algebra over X and by P(X ) the set of probability distributions on X . For a random variable x taking values in X , x ∼ Px ∈ P(X ) represents the probability measure associated to x. Let fn,k (p) = nk pk (1 − p)n−k denote throughout the paper the probability mass function of a binomial random variable with parameters n and k. For a vector x ∈ X , the k-th order statistic of x is denoted by x(k) . Further, we define a stochastic predictor as follows. 2

Definition 1 (Stochastic Predictor). Consider a deterministic function hw : X → Y with X ⊆ Rd and Y ⊂ Rq , where w ∈ Rm is a parameter vector. Then, for x ∈ X and a possibly input dependent random vector w ∼ Pw|x on Rm , the stochastic prediction induced by hw and w at x is the random variable hw (x). Our definition of stochastic predictor is sufficiently general to encompass important stochastic neural models such as Bayesian neural networks (BNNs) [43], variational autoencoders (VAEs) [10], and normalizing flows [28]. In what follows, for simplicity, we assume that the distribution of w does not depend on the input x. Accordingly, we write Pw instead of Pw|x . However, our methods in Sections 3 and 4 also apply to the case where w depends on the input x. Randomized smoothing. Although the literature on randomized smoothing has been extensively developed for classification problems [8, 19, 46, 18], recent work has recently sparked interest in regression settings [7, 31, 30]. Given a stochastic predictor hw and x ∈ X , randomized smoothing can be viewed as the transformation H(x) = Eε [hw (x + ε)], i.e., the original stochastic function is replaced by its expectation w.r.t. noisy perturbations around x. In practice, computing H(x) in closed-form is often infeasible [14], thus commonly replaced by its sampled approximation PN HN (x) = N1 i=1 hw (x + εi ), where {ε1 , . . . , εN } are i.i.d. samples from the distribution of ε. However, as discussed in [30], such smoother presents a major weakness: the smoothed function is vulnerable to outliers in the prediction, with a single sample being capable of pushing the smoothed output beyond the desired range R. Therefore, [30] proposes a α-trimming smoothing technique, which allows trimming these outliers. The α-smoothed predictor is defined as follows. Definition 2 (α-smoothing [30]). Let hw be a stochastic predictor with X , Y, being its input and output spaces and x ∈ X . Given N ∈ N>0 i.i.d. samples {ε1 , . . . , εN } from a distribution Pε ∈ P(X ⊖ {x}), define the random vector H(x) = hw (x + εi ) i∈[N ] . Then, for α ∈ [0, 12 ), the α-smoothed predictor is defined as 1 H̃N,α (x) = N − 2⌊αN ⌋

N −⌊αN ⌋

i=1+⌊αN ⌋

H(x)(i) ,

(1)

where the i-th order statistic in H(x)(i) is taken element-wise. α-smoothing allows tighter control over the probability that H̃N,α (x′ ) ∈ R for any convex set R ⊆ Y. However, when the predictor distribution is multi-modal, the above-mentioned technique inadvertently collapses its modes (Section 3). To solve this issue, we propose clustered α-smoothing. Before introducing clustered α-smoothing, we should note that clustering has previously been applied successfully in randomized smoothing to robustify large language models [39, 42] under the justification that adversarial examples are often located in small clusters adjacent to the dominant mode. However, unlike these previous works, the safety-critical settings considered in this manuscript require that we account for behavior of the predictor in those minor clusters and thus cannot simply prune them, as done in [39, 42].

Clustered α-Smoothing

In this section, we introduce a novel randomized smoothing method tailored to stochastic predictors, which we call clustered α-smoothing. The key idea is to (i) cluster samples in the output space, (ii) apply smoothing locally within each cluster, and (iii) combine the resulting predictors into a mixture distribution, with weights proportional to the size of each cluster. By combining the locallysmoothed predictors as a mixture, as illustrated in Figure 2, we prevent the mode collapse observed in α-smoothing [30]. To formally define clustered α-smoothing, we let V = {V1 , . . . , VM } be a partition of the output space Y into M disjoint sets, each corresponding to a mode of the prediction. For now, we assume that V is given; in Section 4, we describe how to construct V via clustering so that it aligns with the modes of the predictor. Definition 3 (Clustered α-smoothing). Let hw be a stochastic predictor, x ∈ X a point in its domain, and V = {V1 , . . . , VM } a partition of Y. Given N ∈ N>0 i.i.d. samples {ε1n, . . . , εN }

from a distribution Pε ∈ P(X ⊖ {x}), we define the set of indexes Im ⊂ [N ] as Im = i ∈ [N ] : 3

1.0

(a)

Original model α-smoothing Robustness region

Density

0.8

(b)

Original (cluster 1) Original (cluster 2) Clustered α-smoothing Robustness region Voronoi cell split

0.6 0.4 0.2 0.0 −2

−2

Output

Figure 2: An example of α-smoothing and clustered α-smoothing at x = 2 of a stochastic base predictor P hw (x) = wx where w ∼ 3i=1 πi N (· | µi , σi2 ) is a mixture of three Gaussians with means µ = (1, 0, 2) and equal variance σi2 = 0.01, and weights π = (0.2, 0.2, 0.6). The histograms are obtained by Monte Carlo sampling the distributions of hw (x + ε), H̃N,α,{Y} (x) (in (a), equivalent to α-smoothing, see Remark 1), and our clustered α-smoothing technique, H̃N,α,V (x), in (b) constructed by Algorithm 2. The Voronoi cell split represents the partition V = {V1 , V2 } into two clusters, and the robustness regions in (b) denote sets that H̃N,α,V (x) belongs to with high probability (see Algorithm 1).

o hw (x + εi ) ∈ Vm . Then, define HVm (x) = hw (x + εi ) i∈Im and, for α ∈ [0, 12 ) call H̃N,α,Vm (x) =

1 |Im | − 2⌊α|Im |⌋

|Im |−⌊α|Im |⌋

i=1+⌊α|Im |⌋

HVm (x)(i) ,

(2)

where the i-th ordering statist in HVm (x)(i) is taken element-wise. Then, the clustered α-smoothing is defined as follows M X H̃N,α,V (x) = 1z=m H̃N,α,Vm (x) (3) m=1

where z is a categorical random variable with class weights π (m) = |INm | . Intuitively, the set of indexes Im identifies to which partition cell each prediction sample hw (x + εi ) belongs. Then H̃N,α,Vm applies the α-trimming average to each partition cell individually. Finally, the clustered α-smoother H̃N,α,V (x) is constructed by randomly selecting each of the H̃N,α,Vm with a probability proportional to how many samples are in each Vm . The parameter α in Definition 3 refers to trimming a 2α fraction of samples, i.e., treating them as P outliers, per cluster. When α = 0, 1 H̃N,α,Vm (x) reduces to the mean smoothing, H̃N,0,Vm (x) = |Im | i∈Im hw (x + εi ). When α → 12 , H̃N,α,Vm (x) is equivalent to the median-smoothing, H̃N,α,Vm (x) = mediani∈Im hw (x + εi ) . Remark 1 (Relation with α-smoothing in [30]). When V = {V} = {Y}, then Definition 3 is equivalent to the α-smoothing in (1). However, when distributions exhibit multi-modality or are not strongly symmetric, then the flexibility of our approach can lead to a substantial improvement of performance over [30], as we will show in Section 5.

Robustness Guarantees. A core advantage of randomized smoothing is the ability to guarantee that, for sufficiently small perturbations, the output remains within a desired set with high probability. In this subsection, we develop this type of guarantee for clustered α-smoothing. Our approach, which results in Theorem 3.1 below, is to apply the Neyman-Pearson Lemma [8] to each cluster/mode, which are then combined to obtain worst-case bounds for the smoothed predictor via a union-bound argument combined with correction to the fact that V is a partition of X .To apply the lemma, we follow the predominant convention in the randomized smoothing literature and adopt ε ∼ N (0, σ 2 I) for σ 2 > 0, which requires X = Rd . Theorem 3.1 (Robustness certification for clustered α-smoothing). Let V = {V1 , . . . , VM } be a partition of the output space Y. Assume that p̌Vm ≤ P(w,ε) (hw (x + ε) ∈ Vm ) ≤ p̂Vm for each Vm . Furthermore, let L ⊂ [M ] and define R̃ = ∪l∈L Rl for convex sets Rl ⊂ Vl . Assume that P(w,ε) (hw (x + ε) ∈ Rl ) ≥ p̌Rl for each Rl . For r ≥ 0, define for each m ∈ [M ] and l ∈ L r r r pV = Φ Φ−1 (p̌Vm ) − , p̄Vm = Φ Φ−1 (p̂Vm ) + , and pR = Φ Φ−1 (p̌Rl ) − . m l σ σ σ (4) 4

pV = 0.8

pV = 1.0

1.0 0.5 0.0

N = 10 N = 25 N = 50

α = 0.2

1.0 0.5 0.0

α = 0.4

Certified probability

pV = 0.5

α = 0.0

pV = 0.2

1.0 0.5 0.0

0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

Ratio pR /pV

P Ps pR s Figure 3: Probability bounds N fN,s (pV ) from (6) attributable to each s=1 N j=s−⌊αs⌋ fs,j pV mode for different combinations of number of smoothing samples N , trimming parameter α, partition set probability pV and coverage ratio pR /pV . The right-most column (where pV = 1) is equivalent to the bound in [30], by noting that only the term s = N in (6) is non-zero when pV = pV = p̄V = 1.

Then, for any δ ∈ Rd such that ∥δ∥2 ≤ r, it holds that P(w,ε1 ,...,εN ) H̃N,α,V (x + δ) ∈ R̃ ≥

N XX s

inf

s X

N pVm ∈[p ,p̄Vm ] l∈L s=1 j=s−⌊αs⌋ P M Vm p =1 V m m=1

fs,j

(5)

fN,s (pVl ) .

pVl

The proof of Theorem 3.1 can be found in Appendix A. For each partition region Vm , the interval [pV , p̄Vm ] bounds the probability that hw (x + δ + ε) ∈ Vm for any perturbation δ with norm m ∥δ∥2 ≤ r. Analogously, pR represents a probability lower bounds on the event hw (x + δ + ε) ∈ Rl . l Consequently, Theorem 3.1 combines these bounds to obtain a lower bound for the probability that H̃N,α,V (x + δ) ∈ R̃. Figure 3 shows how the various parameters affect the resulting bounds. It can be observed that a higher trimming parameter α leads to a significantly faster increase of the bound, p l as outlier output values are eliminated from the average predictor in (2), while higher ratios pR – i.e. Vl larger coverage in the partition set – naturally increases the bound. Before showing how to compute (5) in practice, it is instructive to examine the properties of the case |L| = 1 first. In particular, the objective function in this case is independent from all other indices [M ] \ L, and so, as we formally prove in Lemma A.1, the bound in (5) simply reduces to P(w,ε1 ,...,εN ) H̃N,α,V (x + δ) ∈ Rl p N s n o X X s Rl ≥ fs,j min fN,s pV , fN,s (p̄Vl ) . l N p̄Vl s=1

(6)

j=s−⌊αs⌋

Solving the general case. Unfortunately, the computational convenience of the case |L| = 1 does not extend to the general case |L| > 1, as in this setting (5) consists of a non-convex optimization problem. Instead, to achieve a computationally tractable result, we show how the optimization problem in (5) can be lower bounded by a linear program via anchor points and an error that can be controlled by the user. This is presented in Proposition 3.2. Proposition 3.2. Assume that pV , p̄Vm for every m ∈ [M ] and pR for every l ∈ L are given. m

(1)

(2)

(K)

(k)

Define for each component l ∈ L the anchor points pVl < pVl < · · · < pVl , pVl ∈ [pV , p̄Vl ]. l

Then, it holds that N XX s

inf

s X

N ,p̄Vm ] pVm ∈[p l∈L s=1 j=s−⌊αs⌋ PM Vm m=1 pVm =1 ≥

inf

K XX

p, λ1 ,...,λ|L|

s.t.

l∈L k=1

M X m=1

pVl

X l∈L

(7)

fN,s (pVl ) "

(k) (k) λl gl (pVl ) −

pVm ∈ [pV , p̄Vm ], m

fs,j

N N + 1 − p̄Vl pV − p R

pVm = 1,

K X k=1

# (8)

hl,K l

(k) (k)

(9)

λl pVl = pVl , λl ∈ ∆K , (k)

where the variable hl,K and function gl are given by hl,K = maxp∈[p ,p̄Vl ] mink∈[K] |p − pVl | and Vl p PN Ps l gl (p) = s=1 Ns fN,s (p) j=s−⌊αs⌋ fs,j R . p The intuition behind Proposition 3.2 is the following. We seek a tractable lower bound on the infimum in (7). This is generally intractable. Consequently, instead, we replace pRl by pR to obtain a valid l lower bound, and then approximate the convex envelope of each inner component using a finite set of evaluation points – the anchor points – leading to a linear program in which the coupling constraint is handled natively. The bound is asymptotically exact as the number of anchor points K increases, at the cost of increasing the size of the linear program. We call the second term of (8) the residual, which represents an adjustment needed due to the fact that we only approximate the convex envelope. Note that, as we show in Appendix G, one can always achieve an arbitrarily small residual value by increasing K. Remark 2 (Extension to Lρ -norm attack). Given ρ ∈ [2, ∞) ∪ {∞}, the guarantee (5) also holds 1 1 for δ ∈ Rd such that ∥δ∥ρ ≤ d ρ − 2 r by a straight-forward application of Holder’s inequality [34].

Algorithm

In Theorem 3.1, given a point x ∈ X , we need to compute probability bounds for the events hw (x + ε) ∈ Vm and hw (x + ε) ∈ Rl , where V = {Vm }M m=1 is a partition of the output space Y into M sets, and Rl ⊆ Vl are given subsets. Given the (possibly) non-linearity of the underlying deterministic function hw on both the parameters w and input x, computing those quantities explicitly is generally intractable. Thus, one needs to resort to sample-based methods, for which probability bounds are estimated with high-confidence. Proposition 4.1 (High-confidence probability bounds for sets). Let V = {Vm }M m=1 be a partition of the space X , and Rl ⊆ Vl be given subsets. Given a confidence parameter β ∈ [0, 1], and a number PN̄ PN̄ of samples N̄ ∈ N>0 , let ZVm = i=1 1hwi (x+εi )∈Vm and ZRm = i=1 1hwi (x+εi )∈Rm , where (wi , εi ) are i.i.d. samples from P(w,ε) . If p̌Vm , p̂Vm , p̌Rm ∈ [0, 1] are such that1 ZVm N̄ N̄ X X X β β β = fN̄ ,j (p̌Vm ) , = fN̄ ,j (p̂Vm ) , and = fN̄ ,j (p̌Rm ) , 3M 3M 3M j=0 j=ZVm

(10)

j=ZRm

then, the assumptions in Theorem 3.1 hold with confidence 1 − β. Proposition 4.1 comes from the Clopper-Pearson Lemma [8] combined with a union-bound argument. It gives a natural way to construct the partition V and sets R: after observing the samples {hwi (x + εi )}, one may place the sets Vm and Rm so that the number of outputs inside them (approximately) induces the desired probability level. This reasoning is the core principle of Algorithm 1. Given the user’s desired coverage level p, Algorithm 1 generates the sets Rm to cover a fraction p of the samples within Vm (Line 9). This parameter can be used to calibrate the trade-off between the size of the coverage sets and their certified probability. Next, by choosing the partition of the input 1 If Z Vm = 0 (resp. ZRm = 0), then p̌Vm = 0 (resp. p̌Rm = 0). Alternatively, if ZVm = N̄ , then p̂Vm = 1.

Algorithm 1 Construction of partition V = V1 , . . . , VM and sets R = R1 , . . . , RM samples such that Ri ⊂ Vi and P(hw (x + ε) ∈ Ri )/P(hw (x + ε) ∈ Vi ) ≈ p.

from

1: function C ONSTRUCT C OVERAGE S ETS(x, hw , p, N ) 2: Outputs ← ∅ 3: for i ← 1, . . . , N do 4: Draw εi ∼ N (0, σ 2 I) and wi ∼ Pw 5: Compute hwi (x + εi ) and append to Outputs 6: Cluster the elements of Outputs ▷ e.g., using DBSCAN, DP-means, UMAP. 7: Denote the resulting clusters by C1 , . . . , CM 8: for m ← 1, . . . , M do 9:

10: 11: 12:

(p)

Cm ← Select approximately a fraction p of the points in Cm (p) Rm ← E NVELOPE(Cm ) ▷ e.g., convex hull, hyperrectangle. V ← generalized Voronoi partition associated with R return partition V and coverage regions R

space as a generalized Voronoi partition (closest region instead of point, Line 11; see Appendix E for details) [25], we enforce that the regions are contained in their partition, which is a requirement for Theorem 3.1. Although computing the (generalized) Voronoi partition associated with a set of regions R is hard (see Algorithm 2), we only need to check whether a point belongs to Vm , which is equivalent to checking if Rm is the closest region to this point. Thus, in practice, V is implicitly defined by R. With partition V, we can compute the robustness bounds introduced in Theorem 3.1 for the clustered α-smoother H̃N,α,V , as presented in Algorithm 2. Appendix C expands on the parameters of Algorithm 2 and their impact on the coverage regions and certified probabilities. Algorithm 2 Robustness certification of H̃N,α,V at x ∈ X 2 .

Require: Point x ∈ Rd , stochastic function hw , smoothing parameters α ∈ [0, 1] and N ∈ N>0 , coverage level p ∈ [0, 1], robustness radius r ≥ 0, number of samples N̄ ∈ N>0 , and confidence level β ∈ (0, 1). Ensure: With confidence 1 − β, for any m ∈ [M ], H̃N,α,V (x + δ) ∈ R̃ with probability given by (5) for any δ ∈ Rd such that ∥δ∥2 ≤ r. 1: V, R ← C ONSTRUCT C OVERAGE S ETS(x, hw , p, N̄ ) 2: Construct H̃N,α,V (x) 3: Compute p̌Vm , p̂Vm , p̌Rm ∈ [0, 1] with Proposition 4.1 for m ∈ [M ] (using β and N̄ ) , p̄Vm , pR using (4) for m ∈ [M ] 4: Compute p V m

5: return H̃N,α,V (x) and (5) certificate.

The time complexity of Algorithm 2 is generally dominated by the repeated inference cost of the base classifier. It is thus governed by the number of samples required, which is determined by (i) the data patterns and clustering algorithm, (ii) the confidence bound (binomial inference, see Proposition 4.1), and (iii) the probability bound in Theorem 3.1.

Experiments

We evaluated our framework on two benchmarks: stochastic trajectory prediction and multi-modal RL quadrotor control. For both experiments, we use the minimum bounding hyperrectangle for the envelope of each cluster (Appendix D provides further details on the experiment setup)3 . All experiments were run on an Intel Core i7-1365U CPU with 16GB of RAM. Trajectory Prediction. We consider the stochastic TrajFlow prediction model [26], a normalizing flow-based model (details in Appendix D.1), and apply it to the L-GAP driving simulator dataset [47]. 2 The quality of the smoothed predictor H̃ N,α,V (x) depends critically on the choice of partition V. In Algorithm 1, we construct V locally around x to ensure good performance in its neighborhood. However, to preserve the validity of Theorem 3.1, the same partition must be used for all points x′ within a δ-neighborhood of x. 3 The code is publicly available at https://github.com/EduardoFMDCosta/ClusteredRandomizedSmoothing.

RS-Reg

α-smoothing

Clustered α-smoothing

Risk rate

1.0 16.5%

P (vehicle stays)

0.8 0.6 0.4

7.5%

0.2 2.5%

0.0 Vehicle Stays

Vehicle Goes

Vehicle Stays

Vehicle Goes

Vehicle Stays

Vehicle Goes

RS-Reg

α-smooth.

Clustered α-smooth.

Figure 4: Safety analysis for 200 randomly selected inputs xi from the L-GAP dataset.

In this dataset, human drivers face a decision whether to turn left in front of an oncoming vehicle (go) or wait for it to pass before turning (yield). The recorded trajectories are converted into control inputs (acceleration and curvature) using the unicycle model [37]. The resulting stochastic predictor is then defined as hw : R18 → R40 , where the input space R18 encodes the two control actions that underlie the 9 transitions between the 10 past recorded positions of the target agent. The output space R40 corresponds to the position of the target vehicle over 20 time-steps, here corresponding to a horizon of 2 s. We start with a safety analysis comparing our approach to α-smoothing and RS-Reg [31] (which is equivalent to smoothing without any trimming, i.e. α = 0 in (1)) using the L-GAP dataset. In particular, we evaluate the predicted probability that the vehicle will adopt the behavior stay4 and we make the ego vehicle go if it predicts a stay probability for the other higher than a decision threshold of 95%. We then use the actual trajectory of the predicted vehicle to assess the risk of this decision: if both vehicles decided to go, we classify it as a risky decision. We then report the risk rate for 200 inputs in our database. We show that RS-Reg has a risk rate of 7.5%, α-smoothing of 16.5%, while ours only 2.5%. The main mechanism behind this result is that while α-smoothing and RS-Reg are often highly sure of the oncoming vehicle’s behavior (due to the mode collapse for both RS-Reg and α-smoothing), our clustered smoother allows for behavioral uncertainty that makes it more cautious about adopting the go decision, as illustrated in Figure 4. Second, for a quantitative comparison of the discrepancy between the output distributions from each approach, we select a random subset of 200 input trajectories xi from L-GAP and compute a sample approximation of the 2-Wasserstein distance W2 between the noisy predictor hw (xi + ε) and α-smoothing, RS-Reg, and ours − note that we are thus comparing distributions embedded in a 40-dimensional space. We found that our clustered α-smoothing reduces the mean 2-Wasserstein distance to hw (x + ε) in this subset of data in 27% compared to α-smoothing (3.80 vs. 5.18; Table 1). This confirms the intuition that the distributional flexibility of our method generally yields predictions that are closer to the original predictor, thereby better capturing its uncertainty (as also illustrated in Figure 7 in Appendix B). Wasserstein W2

α=0 α = 0.2 α = 0.4

Method

Trajectory (R40 )

Last step (R2 )

Ours RS-Reg α-smoothing

3.80 ± 1.59 4.70 ± 2.09 5.18 ± 2.19

1.49 ± 0.63 1.84 ± 0.84 2.04 ± 0.84

Certified lower-bound

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

Robustness radius r

Table 1: Comparison between smoothers and hw (xi + ε) in terms of W2 for 200 inputs xi .

Figure 5: Robustness lower-bounds (with one std. dev. shaded band) for 20 distinct inputs xi .

Finally, we perform a parameter analysis of the robustness lower-bounds in Theorem 3.1 for various trimming parameters α ∈ {0.0, 0.2, 0.4} and robustness radii r ≥ 0, where the sets R are constructed according to Algorithm 1 enveloping a fraction p = 0.85, for different inputs in our database. In Figure 5, we observe that certified lower bounds are higher when the robustness radius is smaller, as this allows smaller perturbations from the reference input, and when the trimming parameter α is higher (so that outliers are filtered out, so we can be more certain about where the smoother 4 The behaviors stay and go are defined in the Appendix D (Behavioral clustering).

(a) nominal policy

(b) α-smoothing [30]

Figure 6: Top row: rollouts of the quadrotor from the starting state (blue marker) where the controller predicts a bi-modal action distribution, corresponding to navigating through two “holes” in the obstacles towards the goal region (green). Bottom row: action distributions at the initial state x = (−3.5, 0, −0.35, 0, 0.7, 0). The smoothed predictors are constructed using 30 samples. Insets: the certified probability of robustness of the smoothed policies with respect to their coverage regions (boxes) as a function of perturbation radius. The quantization in (b) arises from cross-mode smoothing with a limited number of samples (30 samples, α-trimmed to 6) (analyzed further in Appendix F).

probability mass is concentrated). Despite the solid certification suggested by Figure 5, decreasing the conservativeness of the lower bounds in Theorem 3.1 constitutes an interesting research direction. Multi-modal quadrotor control. Quadrotors are safety-critical systems where robustness directly affects deployability. We consider a navigation task in which the quadrotor must reach a goal while avoiding numerous obstacles (see Figure 6), under stochastic dynamics with full state observability [3]. Details of benchmark can be found in Appendix D.2. The task is inherently multi-modal: obstacles between the quadrotor and its goal admit passage on either side, and wind can shift probability mass between such modes. The controller is an RL policy that outputs a Gaussian mixture over actions with predicted component means and mixture weights, and fixed variances [32]. The prediction of a mixture enables navigation around obstacles through either path and rejection of wind-induced disturbances. This is precisely the setting where randomized and α-smoothing fails: perturbations are averaged such that different modes are collapsed to their (weighted) mean, which often is to navigate directly into an obstacle. Clustered smoothing recovers meaningful guarantees by treating each cluster individually. Visual inspection of policies (Figure 6) shows that the nominal policy is bi-modal, corresponding to navigating through one of two paths through the obstacles to the goal. However, the α-smoothing destroys the multi-modality and ultimately predicts the (trimmed) average, which corresponds to navigating directly at the obstacles. Instead, our framework preserves the modes and only tightens the clusters. In the simulation shown in Figure 6, the nominal policy reaches the goal in 72 of the 100 trajectories sampled and crashes in 28 trajectories. The α-smoothed policy crashes in 48 trajectories, while our clustered α-smoothed policy only crashes in 9 trajectories. The certified probability as a function of perturbation radius in Figure 6 suggests that our method is penalized more for larger input radii r, by the correction of pV , p̄Vm in (4). m

Note that the guarantees of the smoothed predictors are with respect to the controller output at each time step, while the benchmark operates in a receding horizon fashion. Nonetheless, the empirical results show significant benefit in applying our method compared to both the base policy and regular α-smoothing. Moreover, our approach is agnostic to the underlying predictor and is thus applicable to controllers such as MPC and motion planners such as RRT, both of which are important methods in robotic applications. 9

Limitations

Randomized smoothing relies on repeated inference. In certain contexts (e.g. large neural networks), this may be prohibitive, in particular when the base predictor has a non-negligible computational cost. That our method is agnostic to the clustering algorithm is both a strength and a weakness: the quality of the result depends strongly on the choice, and thus, one will need to assess the clustering to choose a suitable algorithm. Theory-wise, we identify four primary limitations. First, defining R as a union of disjoint and convex sets is a key component to enabling the guarantees, but this is restrictive, e.g., if the clusters cannot be cleanly separated as convex sets such as the famous half-moon toy example. The second limitation – that certification is with respect to the smoothed predictor – is key, as one cannot apply our method to robustify the prediction if parity with the base regressor is crucial, such as physics and dynamical systems simulations. Third, in (4), if one wants to reduce the noise but maintain the probability bound, it is necessary to fix r/σ, and thus, r must be proportional to the reciprocal of σ. Fourth, as we have shown empirically, the lower bounds from Theorem 3.1 can become more conservative with the inclusion of new clusters. Consequently, either a lower certified radius or a larger noise level may be needed to guarantee robustness in highly multi-modal settings.

Conclusion

In this work, we addressed a key limitation of randomized smoothing for regression: the collapse of multi-modal predictions under averaging. We proposed clustered α-smoothing, which preserves multi-modality by applying α-trimmed smoothing within clusters and combining the results as a mixture. We established probabilistic robustness guarantees and demonstrated significant empirical gains, reducing the Wasserstein distance by 27% in trajectory prediction and the collision rate by 81% in quadrotor control. Several directions remain for future work. In particular, our current approach assumes a fixed partition of the output space V. However, recent work in randomized smoothing [40, 2] suggests that input-dependent smoothing distributions can improve performance. Extending our framework to allow input-dependent smoothing and partitions V(x) is therefore a promising and challenging direction.

References [1] Steven Adams, Andrea Patane, Morteza Lahijanian, and Luca Laurenti. Bnn-dp: robustness certification of bayesian neural networks via dynamic programming. In International Conference on Machine Learning, pages 133–151. PMLR, 2023. [2] Motasem Alfarra, Adel Bibi, Philip HS Torr, and Bernard Ghanem. Data dependent randomized smoothing. In Uncertainty in Artificial Intelligence, pages 64–74. PMLR, 2022. [3] Thom Badings, Licio Romao, Alessandro Abate, David Parker, Hasan A Poonawala, Marielle Stoelinga, and Nils Jansen. Robust control for dynamical systems with non-gaussian noise via formal abstractions. Journal of Artificial Intelligence Research, 76:341–391, 2023. [4] Inhwan Bae, Young-Jae Park, and Hae-Gon Jeon. Singulartrajectory: Universal trajectory predictor using diffusion model. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 17890–17901, 2024. [5] Holger Caesar, Varun Bankiti, Alex H. Lang, Sourabh Vora, Venice Erin Liong, Qiang Xu, Anush Krishnan, Yu Pan, Giancarlo Baldan, and Oscar Beijbom. nuScenes: A Multimodal Dataset for Autonomous Driving. In IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11618–11628. IEEE Computer Society, June 2020. [6] Yipu Chen, Haotian Xue, and Yongxin Chen. Diffusion policy attacker: Crafting adversarial attacks for diffusion-based policies. Advances in Neural Information Processing Systems, 37:119614–119637, 2024. [7] Ping-yeh Chiang, Michael Curry, Ahmed Abdelkader, Aounon Kumar, John Dickerson, and Tom Goldstein. Detection as regression: Certified object detection with median smoothing. Advances in Neural Information Processing Systems, 33:1275–1286, 2020. [8] Jeremy Cohen, Elan Rosenfeld, and Zico Kolter. Certified adversarial robustness via randomized smoothing. In international conference on machine learning, pages 1310–1320. PMLR, 2019. [9] Ruilong Deng, Gaoxi Xiao, Rongxing Lu, Hao Liang, and Athanasios V Vasilakos. False data injection on state estimation in power systems—attacks, impacts, and defense: A survey. IEEE transactions on industrial informatics, 13(2):411–423, 2016.

[10] P Kingma Diederik and Welling Max. An introduction to variational autoencoders. Foundations and Trends® in Machine Learning, 12(4):307–392, 2019. [11] Hai Duong, ThanhVu Nguyen, and Matthew B Dwyer. Neuralsat: A high-performance verification tool for deep neural networks. In International Conference on Computer Aided Verification, pages 409–423. Springer, 2025. [12] Andrew Foong, David Burt, Yingzhen Li, and Richard Turner. On the expressiveness of approximate inference in bayesian neural networks. Advances in Neural Information Processing Systems, 33:15897–15908, 2020. [13] Timon Gehr, Matthew Mirman, Dana Drachsler-Cohen, Petar Tsankov, Swarat Chaudhuri, and Martin Vechev. Ai2: Safety and robustness certification of neural networks with abstract interpretation. In 2018 IEEE symposium on security and privacy (SP), pages 3–18. IEEE, 2018. [14] Agathe Girard, Carl Rasmussen, Joaquin Q Candela, and Roderick Murray-Smith. Gaussian process priors with uncertain inputs application to multiple-step ahead time series forecasting. Advances in neural information processing systems, 15, 2002. [15] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy. Explaining and harnessing adversarial examples. arXiv preprint arXiv:1412.6572, 2014. [16] Guy Katz, Clark Barrett, David L Dill, Kyle Julian, and Mykel J Kochenderfer. Reluplex: An efficient smt solver for verifying deep neural networks. In International conference on computer aided verification, pages 97–117. Springer, 2017. [17] Mathias Lecuyer, Vaggelis Atlidakis, Roxana Geambasu, Daniel Hsu, and Suman Jana. Certified robustness to adversarial examples with differential privacy. In 2019 IEEE symposium on security and privacy (SP), pages 656–672. IEEE, 2019. [18] Alexander Levine and Soheil Feizi. Wasserstein smoothing: Certified robustness against wasserstein adversarial attacks. In International conference on artificial intelligence and statistics, pages 3938–3947. PMLR, 2020. [19] Bai Li, Changyou Chen, Wenlin Wang, and Lawrence Carin. Certified adversarial robustness with additive noise. Advances in neural information processing systems, 32, 2019. [20] Xuanqing Liu, Minhao Cheng, Huan Zhang, and Cho-Jui Hsieh. Towards robust neural networks via random selfensemble. In Proceedings of the european conference on computer vision (ECCV), pages 369–385, 2018. [21] Yulong Lu and Jianfeng Lu. A universal approximation theorem of deep neural networks for expressing probability distributions. Advances in neural information processing systems, 33:3094–3105, 2020. [22] Aleksander Madry, Aleksandar Makelov, Ludwig Schmidt, Dimitris Tsipras, and Adrian Vladu. Towards deep learning models resistant to adversarial attacks. In 6th International Conference on Learning Representations, ICLR, 2018. [23] Charles Meyers, Tommy Löfstedt, and Erik Elmroth. Safety-critical computer vision: an empirical survey of adversarial evasion attacks and defenses on computer vision systems. Artificial Intelligence Review, 56(Suppl 1):217–251, 2023. [24] Rhiannon Michelmore, Matthew Wicker, Luca Laurenti, Luca Cardelli, Yarin Gal, and Marta Kwiatkowska. Uncertainty quantification with statistical guarantees in end-to-end autonomous driving control. In 2020 IEEE international conference on robotics and automation (ICRA), pages 7344–7350. IEEE, 2020. [25] Victor Milenkovic. Robust construction of the voronoi diagram of a polyhedron. In Canadian Conference on Computational Geometry, volume 93, pages 473–478, 1993. [26] Anna Mészáros, Julian F. Schumann, Javier Alonso-Mora, Arkady Zgonnikov, and Jens Kober. TrajFlow: Learning Distributions over Trajectories for Human Behavior Prediction. In 2024 IEEE Intelligent Vehicles Symposium (IV), Jeju, June 2024. [27] Nigamaa Nayakanti, Rami Al-Rfou, Aurick Zhou, Kratarth Goel, Khaled S. Refaat, and Benjamin Sapp. Wayformer: Motion Forecasting via Simple & Efficient Attention Networks. In 2023 IEEE International Conference on Robotics and Automation (ICRA), pages 2980–2987, May 2023. [28] George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1–64, 2021. [29] Zeyu Qin, Yanbo Fan, Yi Liu, Li Shen, Yong Zhang, Jue Wang, and Baoyuan Wu. Boosting the transferability of adversarial attacks with reverse adversarial perturbation. Advances in neural information processing systems, 35:29845– 29858, 2022. [30] Aref Rekavandi, Farhad Farokhi, Olga Ohrimenko, and Benjamin Rubinstein. Certified adversarial robustness via randomized α-smoothing for regression models. Advances in Neural Information Processing Systems, 37:134127– 134150, 2024.

[31] Aref Miri Rekavandi, Olga Ohrimenko, and Benjamin IP Rubinstein. Rs-reg: Probabilistic and robust certified regression through randomized smoothing. Transactions on Machine Learning Research, 2025. [32] Jie Ren, Yewen Li, Zihan Ding, Wei Pan, and Hao Dong. Probabilistic mixture-of-experts for efficient deep reinforcement learning. arXiv preprint arXiv:2104.09122, 2021. [33] Andy Rivas, Gregory Kyriakos Delipei, and Jason Hou. Predictions of component remaining useful lifetime using bayesian neural network. Progress in Nuclear Energy, 146:104143, 2022. [34] Walter Rudin. Real and complex analysis, 1987. [35] Hadi Salman, Jerry Li, Ilya Razenshteyn, Pengchuan Zhang, Huan Zhang, Sebastien Bubeck, and Greg Yang. Provably robust deep learning via adversarially trained smoothed classifiers. Advances in neural information processing systems, 32, 2019. [36] Julian F Schumann, Eduardo Figueiredo, Frederik Baymler Mathiesen, Luca Laurenti, Jens Kober, and Arkady Zgonnikov. Evaluating randomized smoothing as a defense against adversarial attacks in trajectory prediction. arXiv preprint arXiv:2603.10821, 2026. [37] Julian F. Schumann, Jeroen Hagenus, Frederik Baymler Mathiesen, and Arkady Zgonnikov. Realistic Adversarial Attacks for Robustness Evaluation of Trajectory Prediction Models via Future State Perturbation. ACM Journal on Autonomous Transportation Systems, 2026. [38] Julian F Schumann, Anna Mészáros, Jens Kober, and Arkady Zgonnikov. STEP: Structured training and evaluation platform for benchmarking trajectory prediction models. arXiv preprint arXiv:2509.14801, 2025. [39] Guangzhi Su, Shuchang Huang, Yutong Ke, Zhuohang Liu, Long Qian, and Kaizhu Huang. Smoothguard: Defending multimodal large language models with noise perturbation and clustering aggregation. arXiv preprint arXiv:2510.26830, 2025. [40] Peter Sukenik, Aleksei Kuvshinov, and Stephan Günnemann. Intriguing properties of input-dependent randomized smoothing. In International Conference on Machine Learning, 2021. [41] Hoang-Dung Tran, Xiaodong Yang, Diego Manzanas Lopez, Patrick Musau, Luan Viet Nguyen, Weiming Xiang, Stanley Bak, and Taylor T Johnson. Nnv: the neural network verification tool for deep neural networks and learning-enabled cyber-physical systems. In International conference on computer aided verification, pages 3–17. Springer, 2020. [42] Zixia Wang, Gaojie Jin, Jia Hu, and Ronghui Mu. Clucert: Certifying llm robustness via clustering-guided denoising smoothing. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, pages 37998–38006, 2026. [43] Anqi Wu, Sebastian Nowozin, Edward Meeds, Richard E. Turner, José Miguel Hernández-Lobato, and Alexander L. Gaunt. Deterministic variational inference for robust bayesian neural networks. In 7th International Conference on Learning Representations, ICLR, 2019. [44] Haoze Wu, Omri Isac, Aleksandar Zeljić, Teruhiro Tagomori, Matthew Daggitt, Wen Kokke, Idan Refaeli, Guy Amir, Kyle Julian, Shahaf Bassan, et al. Marabou 2.0: a versatile formal analyzer of neural networks. In International Conference on Computer Aided Verification, pages 249–264. Springer, 2024. [45] Kaidi Xu, Zhouxing Shi, Huan Zhang, Yihan Wang, Kai-Wei Chang, Minlie Huang, Bhavya Kailkhura, Xue Lin, and Cho-Jui Hsieh. Automatic perturbation analysis for scalable certified robustness and beyond. Advances in Neural Information Processing Systems, 33:1129–1141, 2020. [46] Greg Yang, Tony Duan, J Edward Hu, Hadi Salman, Ilya Razenshteyn, and Jerry Li. Randomized smoothing of all shapes and sizes. In International conference on machine learning, pages 10693–10705. PMLR, 2020. [47] Arkady Zgonnikov, David Abbink, and Gustav Markkula. Should I Stay or Should I Go? Cognitive Modeling of Left-Turn Gap Acceptance Decisions in Human Drivers. Human Factors: The Journal of the Human Factors and Ergonomics Society, 66(5):1399–1413, May 2024. [48] Huan Zhang, Tsui-Wei Weng, Pin-Yu Chen, Cho-Jui Hsieh, and Luca Daniel. Efficient neural network robustness certification with general activation functions. Advances in neural information processing systems, 31, 2018. [49] Liangliang Zhuang, Ancha Xu, and Xiao-Lin Wang. A prognostic driven predictive maintenance framework based on bayesian deep learning. Reliability Engineering & System Safety, 234:109181, 2023.

Technical Proofs

A.1

Proof of Theorem 3.1

As Theorem 3.1 is built from a union-bound argument over Rm , we start by introducing the following Lemma A.1, which bounds the probability P(w,ε1 ,...,εN ) H̃N,α,V (x + δ) ∈ Rm . For notational convenience, we write ε to mean ε1 , . . . , εN for the remainder of this appendix. Lemma A.1. Let V = {V1 , . . . , VM } be a partition of the space X . Assume that p̌Vm ≤ P(w,ε) (hw (x + ε) ∈ Vm ) ≤ p̂Vm and P(w,ε) (hw (x + ε) ∈ Rm ) ≥ p̌Rm for a convex set Rm ⊂ Vm . For r ≥ 0, define r r r , p̄Vm = Φ Φ−1 (p̂Vm ) + , and pR = Φ Φ−1 (p̌Rm ) − . pV = Φ Φ−1 (p̌Vm ) − m m σ σ σ (11) Then, for any δ ∈ Rd such that ∥δ∥2 ≤ r, it holds that P(w,ε) H̃N,α,V (x + δ) ∈ Rm p s N n o X X s Rm fs,j ≥ min fN,s pV , fN,s (p̄Vm ) . (12) m N p̄Vm s=1 j=s−⌊αs⌋

Proof. First, we observe that using the Neyman-Pearson’s Lemma [8], for any δ ∈ Rd such that ∥δ∥2 ≤ r, ∥δ∥2 r −1 P(w,ε) (hw (x + δ + ε) ∈ Vm ) ≥ Φ Φ (p̌Vm ) − = pV , (13) ≥ Φ Φ−1 (p̌Vm ) + m σ σ ∥δ∥2 r P(w,ε) (hw (x + δ + ε) ∈ Vm ) ≤ Φ Φ−1 (p̂Vm ) + ≤ Φ Φ−1 (p̂Vm ) + = p̄Vm , (14) σ σ and ∥δ∥2 r P(w,ε) (hw (x + δ + ε) ∈ Rm ) ≥ Φ Φ−1 (p̌Rm ) − = pR , ≥ Φ Φ−1 (p̌Rm ) − m σ σ (15) due to the monotonicity of the function Φ. Further, because by construction Rm ⊂ Vm , we have that p (16) P(w,ε) (hw (x + δ + ε) ∈ Rm | hw (x + δ + ε) ∈ Vm ) ≥ Rm . p̄Vm Having these probability bounds for hw (x + δ + ε), we can now proceed to study the probability of the event H̃N,α,V (x + δ) ∈ Rm . For any x′ ∈ X , by applying the Law of Total Probability twice, it holds that P(w,ε) H̃N,α,V (x′ ) ∈ Rm = P(w,ε) H̃N,α,V (x′ ) ∈ Rm | H̃N,α,Vm (x′ ) ∈ Rm P(w,ε) H̃N,α,Vm (x′ ) ∈ Rm + P(w,ε) H̃N,α,V (x′ ) ∈ Rm | H̃N,α,Vm (x′ ) ∈ / Rm P(w,ε) H̃N,α,Vm (x′ ) ∈ / Rm = | {z } =0

N X s=0

P(w,ε) H̃N,α,V (x′ ) ∈ Rm | H̃N,α,Vm (x′ ) ∈ Rm , |Im | = s · P(w,ε) H̃N,α,Vm (x′ ) ∈ Rm | |Im | = s P(w,ε) (|Im | = s).

(17)

In particular, let x′ = x + δ. We will examine each of these terms separately. Recall that, by construction, the random set of indexes Im is given by Im = {i ∈ [N ] : hw (x + δ + εi ) ∈ Vm }, so that |Ij | ∼ Binomial N, P(w,ε) (hw (x + δ + ε) ∈ Vm ) . Therefore, n o (18) P(w,ε) (|Im | = s) ≥ min fN,s pV , fN,s (p̄Vm ) , m

h i as P(w,ε) (hw (x + δ + ε) ∈ Vm ) ∈ pV , p̄Vm by (13) and (14). That the minimum is attained at m pV or p̄Vm follows from the fact that the continuous function fN,s (p) = Ns ps (1 − p)N −s has m a positive derivative f ′ (p) > 0 for p < Ns (resp. f ′ (p) < 0 for p > Ns ). Hence, p = Ns is the maximizer within [0, 1] and the minimum of f (p) for p ∈ [p, p̄] must be attained at either p = p or p = p̄. We then move on to the event H̃N,α,V (x + δ) ∈ Rm given H̃N,α,Vm (x + δ) ∈ Rm , |Im | = s. Trivially, conditioned on H̃N,α,Vm (x + δ) ∈ Rm , the event H̃N,α,V (x + δ) ∈ Rm takes place if the component H̃N,α,Vm is sampled from the mixture, which by construction is equal to |INm | . Thus, s P(w,ε) H̃N,α,V (x + δ) ∈ Rm | H̃N,α,Vm (x + δ) ∈ Rm , |Im | = s = . (19) N Finally, to conclude, we need to lower-bound the probability of the event H̃N,α,Vm (x+δ) ∈ Rm given |Im | = s. We observe that conditioned on the knowledge |Im | = s, the random variable H̃N,α,Vm (x+ δ) conditioned on |Im | = s behaves as the α-smoother Ps introduced in [30] with number of samples s. Therefore, we follow a similar proof. Let Zm = i=1 1hw (x+δ+εi )∈Rm |hw (x+δ+εi )∈Vm , thus Zm ∼ Binomial s, P(w,ε) (hw (x + δ + ε) ∈ Rm | hw (x + δ + ε) ∈ Vm ) .

(20)

From the Law of Total Probability, it holds that P(w,ε) H̃N,α,Vm (x + δ) ∈ Rm | |Im | = s = s X z=0

P(w,ε) H̃N,α,Vm (x + δ) ∈ Rm | |Im | = s, Zm = z P(Zm = z),

(21)

Now, note that the event z ≥ s − ⌊αs⌋ means that the α-trimming only discards points out of Rm (i.e., the average in (2) only contains points in Rm ), thus it must hold that H̃N,α,Vm (x + δ) ∈ Rm conditioned on the given events by the convexity of the set Rm . On the other hand, when z < N − ⌊αs⌋, the average in (2) contains points out of Rm . In the absence of more information about how far those points are, we have no guarantees that the average over them will belong to Rm (as only one point suffices to displace the average arbitrarily far from Rm [31]). Therefore, conditioned on these given events, we may conservatively consider that H̃N,α,Vm (x + δ) ∈ / Rm . Thus, from (21), P(w,ε) H̃N,α,Vm (x + δ) ∈ Rm | |Im | = s = s−⌊αs⌋−1

P(w,ε) H̃N,α,Vm (x + δ) ∈ Rm | |Im | = s, Zm = z PZm (Zm = z)+ | {z }

X z=0

≥0

s X

P(w,ε) H̃N,α,Vm (x + δ) ∈ Rm | |Im | = s, Zm = z PZm (Zm = z) {z } z=s−⌊αs⌋ | =1

≥ ≥

s X

PZm (Zm = z)

z=s−⌊αs⌋ s X j=s−⌊αs⌋

fs,j

p̄Vm

(22)

from the fact that Zm is distributed as the binomial in (20) and its success probability is lower bounded by (16). The proof is concluded by replacing (18), (19), and (22) in (17). We are now ready to prove Theorem 3.1. 14

Proof. By assumption, we have Rl ∩ Rl′ = ∅ for l, l′ ∈ L where l ̸= l′ . Therefore, we can employ a union-bound argument over Lemma A.1 for the following bound P(w,ε) H̃N,α,V (x + δ) ∈ R̃ X = P(w,ε) H̃N,α,V (x + δ) ∈ Rl l∈L

≥ =

X l∈L

inf

pVl ∈[p

inf

pVm ∈[p

N X s ,p̄Vl ]

s=1

s X

l∈L s=1

s X

fN,s (pVl )

pVl

j=s−⌊αs⌋

N XX s ,p̄Vm ]

fs,j

j=s−⌊αs⌋

pVl

fN,s (pVl ) ,

(23)

where the last infimum is for all m ∈ [M ]. Naturally, one would compute this infimum independently PM per l ∈ L. However, extracting it reveals that we can add the constraint m=1 pVm = 1, i.e., ensure that the probability over all partition cells Vm sum to one. Thus, we arrive at P(w,ε) H̃N,α,V (x + δ) ∈ R̃ p N s XX X s Rl ≥ inf fs,j fN,s (pVl ) , (24) N p pVm ∈[p ,p̄Vm ] Vl Vm s=1 PM

m=1 pVm =1

l∈L

j=s−⌊αs⌋

which concludes the proof. A.2

Proof of Proposition 3.2

Before proving Proposition 3.2, we first prove an auxiliary Lemma. Lemma Lipschitz constant Lg . Let (1) A.2.(K)Let g : [ℓ, u] → R be a Lipschitz continuous function with be a set of K ∈ N>0 anchor points such that ℓ = p(1) < p(2) < · · · < p(K) = u. p ,...,p Let convK (g)(p) be defined as (K ) K X X (k) (k) (k) (k) convK (g)(p) = min λ g(p ) s.t. λ p =p . (25) λ∈∆K

k=1

Then, convK (g)(p) ≤ g(p) + Lg hK , ′

where hK = maxp′ ∈[ℓ,u] mink∈[K] |p − p

(k)

(26)

Proof. For any p ∈ [ℓ, u], there exist consecutive anchor points p(k) , p(k+1) such that p(k) ≤ p ≤ p(k+1) . We can write p as a convex combination of those anchor points, i.e. p = λp(k) +(1−λ)p(k+1) p(k+1) −p by taking λ = p(k+1) ∈ [0, 1]. Thus, the vector λ = [0, . . . , 0, λ, 1 − λ, 0, . . . , 0] (i.e., non-zero −p(k) at positions k and k + 1) belongs to ∆K and satisfies the constraint in (25), thus convK (g)(p) ≤ λg(p(k) ) + (1 − λ)g(p(k+1) ).

(27)

where Lg is the Lipschitz constant of g in [ℓ, u]. Thus, from (27) and (28), we have that, for any p ∈ [ℓ, u], convK (g)(p) ≤ g(p) + Lg hK .

We are now ready to prove Proposition 3.2. Proof. The goal of Proposition 3.2 is to bound the following optimization problem. p N s XX X s Rl fs,j fN,s (pVl ) . inf N p ,p̄Vm ] pVm ∈[p Vl Vm s=1 PM

m=1 pVm =1

l∈L

(29)

j=s−⌊αs⌋

Due to the non-convexity of (29), we will search for a tractable convex lower-bound. First, note that the objective function is separable in pVl as a sum of functions gl , which allows the objective in P (29) to be written as G(p) = l∈L gl (pVl ), where p = {pVl }l∈L . For each component l, select K anchor points: (1) (2) (K) (k) pVl < pVl < · · · < pVl , pVl ∈ [pV , p̄Vl ], (30) l

(k)

and evaluate gl

(k)

:= gl (pVl ) at each anchor. Then, let (K ) K X X (k) (k) (k) (k) convK (gl )(p) = min λ gl (p ) s.t. λ p =p . λ∈∆K

k=1

(31)

k=1

From Lemma A.2, it holds that (32)

convK (gl )(pVl ) ≤ gl (pVl ) + Lgl hl,K ,

where Lgl is the Lipschitz constant of the function gl in the interval [pV , p̄Vl ] and the maximum l

(k)

interval size hl,K = maxp∈[p ,p̄Vl ] mink∈[K] |p − pVl |. Therefore, Vl X X X convK (gl )(pVl ) ≤ inf gl (pVl ) + Lgl hl,K . inf pVm ∈[p ,p̄Vm ] l∈L P M Vm m=1 pVm =1

,p̄Vm ] pVm ∈[p l∈L P M Vm m=1 pVm =1

(33)

l∈L

As the λ variables within the definition of each convK (gl )(pVl ) are independent from each other, we can write the left-hand side term as the following LP problem inf

pVm ∈[p ,p̄Vm ] l∈L PM Vm p m=1 Vm =1

convK (gl )(pVl ) =

p, λ1 ,...,λ|L|

s.t.

To conclude, we apply Lemma A.3 to obtain Lgl ≤

K XX

inf

l∈L k=1

(k)

λl gl (pVl )

 pVm ∈ [pV , p̄Vm ],  m     M  X    pVm = 1,   m=1

K  X   (k) (k)   λl pVl = pVl ,     k=1   λl ∈ ∆K .

(34)

(35)

N N + . 1 − p̄Vl pV − pR l

Lemma A.3. Let g : (0, 1] → R be the function defined by j s−j N s X X s p̃ p̃ N s s 1− p (1 − p)N −s , g(p) = N j p p s s=1 j=s−⌊αs⌋

where α ∈ [0, 21 ), N ∈ N>0 , and 0 < p̃ < p < p̄ < 1. Then g is Lipschitz continuous on [p, p̄] with

Lipschitz constant

Lg ≤

N N + . 1 − p̄ p − p̃ 16

(36)

Proof. Because g is a sum of continuous and differentiable functions on [p, p̄], its Lipschitz constant is given by Lg = supp∈[p,p̄] |g ′ (p)|. For notational simplicity, we write the summand as j s−j p̃ p̃ s N s N −s Ts,j (p) = Cs,j p (1 − p) 1− , Cs,j = . p p j s s j N d 1− First, note that, for ps (1 − p)N −s , dp p (1 − p)N −s ≤ ps (1 − p)N −s 1− p̄ . Then, for p̃/p s−j j s−j N , the absolute value of the derivative is at most p̃/p 1− p̃/p p̃/p p−p̃ . Since all coefficients are nonnegative and sum to at most 1, it holds that N N |g ′ (p)| ≤ + , p ∈ [p, p̄], 1 − p̄ p − p̃

which concludes the proof. A.3

Proof of Proposition 4.1

Before proving Proposition 4.1, we introduce the well-known Clopper-Pearson Lemma. Lemma A.4 (Clopper-Pearson lower bound). Let {z1 , . . . , zM } be i.i.d. samples from P, and let PM h : X → Y and R ⊆ Y be given. Define Z = i=1 1h(zi )∈R . Given a confidence level β ∈ [0, 1], PM if p̌ ∈ [0, 1] is such that β = j=Z fM,j (p̌), then P(z1 ,...,zM ) (Pz (h(z) ∈ R) ≥ p̌) ≥ 1 − β

(37)

We are now ready to prove Proposition 4.1. Proof. From the Fréchet bound and Lemma A.4, P{(wi ,εi )}N ∩M m=1 p̌Vm ≤ P(w,ε) (hw (x + ε) ∈ Vm ) ≤ p̂Vm , p̌Rm ≤ P(w,ε) (hw (x + ε) ∈ Rm ) i=1 3M X β ≥ 1− − (3M − 1) = 1 − β. 3M i=1

Additional Experimental Results

B.1

Trajectory Prediction

We examine a specific input instance from the dataset and demonstrate that the introduction of a small input perturbation ε ∼ N (0, I) increases the behavioral variability of the target agent by raising the probability of the tactical behavior go. While α-smoothing tends to suppress the go mode – thereby suggesting that the target vehicle is likely to yield – our clustered α-trimming preserves the predicted multi-modality of hw (x + ε), as illustrated in Figure 7, where we show the last time step of the sampled trajectories.

Summary of method parameters

The clustered α-smoothing certificate is parameterized by three groups. The clustering group depends on the strategy. For the purposes of the experiments in Section 5, we use DBSCAN. The parameters are listed in Table 2.

Experiment Details

D.1

Trajectory Prediction

TrajFlow architecture and training procedure. TrajFlow [26] is built with a combination of a normalizing flow network F with a recurrent autoencoder, consisting of encoder ERNN and decoder 17

hw (x)

hw (x + ε)

α-smoothing

Clustered α-smoothing

−5

−10 −10

−5

10−10

−5

10−10

−5

10−10

−5

Figure 7: Given a particular trajectory input x from the L-GAP dataset [47], we show the last time step of trajectory samples from hw (x), hw (x + ε), where ε ∼ N (0, I), the α-smoothing in (1), and our clustered α-smoothing. We also show the projection to the last time step subspace of the sets Rm produced with Algorithm 1 for coverage p = 0.8. Table 2: Clustered α-smoothing parameters. Parameter

Symbol

Impact

Clustering (DBSCAN strategy) Maximum number of clusters

Mmax

Coverage Neighborhood radius Min. samples

p εdb –

If the output of DBSCAN has more than Mmax clusters, then we merge them in order of fewest to most samples. See Line 9 of Algorithm 1 DBSCAN neighbourhood radius in output space. DBSCAN core-point threshold.

Smoothing core Std.dev. of Gaussian input-perturbation

Number of samples for smoothing Trimming fraction

N α

Confidence Number of samples for confidence

N̄

Confidence parameter

Larger std.dev. allows larger perturbation radius r, at the cost of larger deviation from hw (x). See (4). More samples tightens the modes of the smoothed distribution. The certificate discards the most extreme 2α tails of each per-cluster output coordinate, 0 ≤ α < 1/2. Number of samples used for Proposition 4.1. More samples leads to better confidence, at an increased computational cost. The joint confidence for all estimated probability bounds is 1 − β.

ERNN (see Figure 8). Furthermore, the model is completed with the networks ϕRNN (a gated recurring unit encoding the predicted agents past trajectory), ψGNN (a graph neural network encoding the other agents), and ϕCNN (a convolutional network encoding the environment, which is given as a bird’s-eye image). This model is trained in a two-stage process, where first the recurrent autoencoder is trained to encode the future trajectories of all involved agents. Once done, the normalizing flow is trained on the encoded future sampled as well as the past and environment information. An detailed description of the model components and the training loss can be found in the original work [26]. For model training, we employ the STEP framework [38] (MIT License) to train TrajFlow for our chosen scenario. In particular, the model is trained on a dataset combining NuScenes [5] (CC BY-NC-SA) and L-GAP [47] (CC-BY). Following previous works [37, 36], we do not use the standard input for NuScenes (i.e., 4 past time steps at a frequency of 2 Hz), but instead use STEP to extract past trajectories with 10 time steps at 10 Hz. This choice is mostly motivated by the fact that this results in a smoother trajectory, for which we can extract the underlying control states much easier. Importantly, while these datasets provide trajectories with states including positions, velocities and accelerations, for TrajFlow, we only consider the recorded positions. The code is available at https://github.com/DAI-Lab-HEReALD/General-Framework/blob/ main/Framework/simulations_randomized_smoothing.py. Behavioral clustering. For our experiments, we run Algorithm 1 with a behavior-based two-cluster heuristic in endpoint space, using prior knowledge of two tactical modes (stay and go). For each sample, let ∆ = yend − yinit . A sample is tagged as stay when ∥∆∥2 ≤ rstay , and as go when ∆(2) ≤ −δ↓ and ∆(1) ≤ −δ← , excluding already-tagged stay samples. Endpoint prototypes are then 18

Normalizing direction y

Past

Generative direction p0 (z0 )

xtar x E

ERNN

b y

DRNN Normal. Flow F

ϕRNN ψGNN ϕCNN

p0 (z0 )

b x

p0 (z0 ) det J −1 (z0 ) F

Past

xtar x E

ϕRNN ψGNN ϕCNN

b y

Normal. Flow F −1 b x

ypred

Future

p0 (z0 ) det J −1 (z0 ) F

Concatenate

Figure 8: An overview of TrajFlow’s architecture, taken from Figure 2 in [26]. During training, future b Using trajectories y are encoded with the encoder ERNN and transformed to the abstracted features y. b which follows a standard the normalizing direction, they are then transformed into a sample z0 = F (y), normal distribution p0 . For inference we then use the generative direction, in which a sample z0 ∼ p0 is b = F −1 (z0 ) inversely transformed by the Normalizing Flow to generate the abstracted future trajectories y b is that are decoded with DRNN into the actual trajectories ypred . The likelihood of the encoded trajectory y obtained with p0 (z0 )|detJF −1 (z0 )|−1 . The encoding ϕCNN of map E, and the encoding ψGNN of social interactions are optional blocks, which can provide richer context information.

built from these tagged sets (with fallback rules if one set is empty), and all samples are assigned to the nearest prototype in 2D endpoint space. Finally, for each cluster j ∈ 1, 2, we construct a tight axis-aligned hyperrectangle Rj in the full trajectory space (40D) that covers a fraction p of the samples assigned to that cluster. D.2

Quadrotor Control

The quadrotor control benchmark was first described in [3]. We use it solely to obtain a safety-critical, stochastic, multi-modal closed-loop system on which to perform clustered randomized smoothing. Relative to [3], we augment the system with a reward structure and extended termination criteria. Environment. The state is sk = [xk , ẋk , yk , ẏk , zk , żk ] ∈ R6 and the action uk = [uxk , uyk , uzk ] ∈ [−4, 4]3 . The dynamics are a per-axis double integrator with unit time step, sk+1 = Ask + Buk + Cwk ,

wk ∼ N (0, σ 2 I)

(38)

with velocities clipped to [−7, 7] after each update. In particular, per axis the dynamics are xk+1 = xk + τ ẋk + 12 τ 2 uxk and ẋk+1 = ẋk + τ uxk + wk where wk ∼ N (0, σ 2 ). The corridor is [−15, 15] × [−9, 9] × [−7, 7], the goal region is [11, 15] × [1, 5] × [−7, −3], and 14 axis-aligned hyperrectangular obstacles match [3]. The termination criteria are in order: 1. if [xk , yk , zk ] is outside the corridor, 2. if any point on the line segment between [xk , yk , zk ] and [xk+1 , yk+1 , zk+1 ] touches an obstacle, 3. if [xk , yk , zk ] is inside the goal region, or 4. if the max episode length is reached. The parameters of the environment are summarized in Table 3. Policy and value networks. The policy is a tanh-squashed Gaussian mixture: at state s, sample a mode z ∼ Cat(softmax(gϕ (s))) where Cat is a categorical distribution, sample an unscaled action from this mode r ∼ N (µzθ (s), diag e2 log σz ), and rescale the action a = tanh(r) · ascale . The mixture is results in a multi-modal policy, which indeed requires clustered α-smoothing. The architecture is summarized in Table 4. Training. The policy is trained using standard clipped PPO with GAE and a curriculum on the start state. The curriculum using a linear schedule where a scalar f ∈ [0, 1] linearly anneals from 0 to 1 over the first 40,000 episodes. At each rollout we sample a valid (collision-free, outside the goal) start position uniformly inside an axis-aligned box around the goal center with half-width r = 3 + 25f , and initial velocities N (0, (0.1 + 0.4f )2 ). Without curriculum learning, the policy 19

Table 3: Environment parameters. ∆dist denotes the change in distance ∥pos − goal∥ to the goal center goal = (13, 3, −5) between subsequent steps.

Quantity

Value

Notes

State / action dim Action bounds Process noise σ Velocity clip Max steps / episode Reward (goal) Reward (collision / exit) Reward (timeout) Reward (alive) State normalization

6/3 [−4, 4]3 0.05 [−7, 7] 64 +10 −5 −0.05 · ∥pos − goal∥ 0.5 · ∆dist − 0.01 s/[15, 7, 9, 7, 7, 7]

matches JAIR benchmark on velocity coords; small enough that obstacles dominate prevents PPO blow-up under exploration matches rollout horizon terminal terminal distance-shaped: far timeouts penalised more dense progress shaping per-dim rescale for the policy net

Table 4: Network details. A denotes the action dimension. Component Trunk Mean head Mixture head (K > 1) Log-std Critic K (mixture components) ascale

Specification

√ 2 × Linear(128) + tanh, orthogonal init (gain 2) Linear(K · A), orthogonal init gain 0.01 Linear(K), orthogonal init gain 0.01 global per-component [K, A], init −1.2, clamped to [−5, 2] identical 2 × Linear(128) + tanh trunk, scalar head 2 4.0 (matches action bound; tanh saturation gives the true control limit)

frequently terminates in collisions before it ever sees the goal. The remaining training parameters are listed in Table 5. Table 5: PPO hyperparameters. Hyperparameter Optimizer Total episodes Rollout horizon T Episodes per update PPO epochs / batch Minibatch size Discount γ GAE λ Clip range Value loss coefficient Entropy coefficient Max grad norm Advantage normalization

Value

Notes −4

Adam, lr 3 × 10 200,000 64 steps 32 10 512 0.99 0.95 0.2 0.5 0.02 0.5 per-batch (mean / std)

shared across actor and critic matches env timeout ≈ 2k transitions / batch

mild; the mixture itself adds exploration global clip across actor + critic

Randomized smoothed policy parameters. Figure 6 is produced with the clustered α-smoothing parameters in Table 6 below. α-smoothing uses the same parameters, except for no clustering.

Generalized Voronoi partition

Given a collection of points C = {c1 , . . . , cM } ⊂ X , the Voronoi partition Ṽ induced by C is defined5 as Ṽi = {x ∈ X : ∥x − ci ∥2 ≤ ∥x − cj ∥2 , i ̸= j}, i.e., the point x belongs to the cell Ṽi if it is closer to ci than to any other point cj in C. The Generalized Voronoi partition extend this notion 5 Although we introduce the Voronoi partition for the Euclidean distance, we must stress that it can be defined for any other

metric.

Table 6: Clustered α-smoothing parameters for the quadrotor benchmarks. Parameter

Value

Notes

σ N α Mmax Coverage p Clustering alg. εdb dbscan_min_samples N̄ β

0.05 30 0.4 3 0.9 DBSCAN 0.45 50 4000 10−2

Output-space DBSCAN radius. Core-point threshold relative. Samples for Proposition 4.1.

(a) d(p, R) = minr∈R kp − rk2

(b) Generalized–Voronoi partition

d(p, R3 )

−2

p R4

−4 −6

−6

−4

−2

−6

−4

−2

Figure 9: Generalized Voronoi partition of four disjoint axis-aligned rectangles R1 , . . . , R4 ⊂ R2 under the Euclidean metric. (a) The distance from point p to the closest point in each region Ri is measured and p is assigned V3 as R3 is the closest. (b) Each cell Vi collects every point closer to Ri than to any other site. Cell boundaries in the case of rectangles in 2D consist of straight segments and parabolic arcs and meet at Voronoi vertices where three cells agree on the distance to a common point.

to sets, provided a shortest-distance measure from x to a set R ⊆ X . We write d(x, R) = inf ∥x − r∥2 r∈R

(39)

for the (Euclidean) point-to-set distance, which equals 0 if and only if x ∈ R for closed R. Definition 4 (Generalized Voronoi partition). For a collection R = {R1 , . . . , RM } of pairwise disjoint closed sets in X , the generalized Voronoi cell associated with Ri is given by Vi = x ∈ X : d(x, Ri ) ≤ d(x, Rj ) for all j ̸= i , and V = {Vi }M i=1 is the generalized Voronoi partition induced by R.

Unlike Ṽ, the cells Vi in the generalized Voronoi partition V are generally not convex, even when the sets Ri are convex or simple shapes such as rectangles, as shown in Figure 9. This fact, together with the combinatorial size of the partition that grows rapidly with both the number of sets M and the dimension of the ambient space X , implies that computing the partition is almost certainly intractable. d Indeed, even for the point-induced partition Ṽ, the worst-case complexity in X ⊆ R is already ⌈d/2⌉ Θ M . When the sets Ri are convex, the cell boundaries are (d − 1)-dimensional algebraic surfaces, and for arbitrary compact sets the topology can be unbounded. Explicit construction of the partition is, therefore, intractable in Rd beyond toy examples. Crucially, however, our procedure does not require the explicit construction of the partition, but only identifying to which cell Vi a given point x ∈ X belongs, which reduces to computing M point-to-set distances (39). This admits a trivial implementation when, for instance, the sets Ri are axis-aligned hyperrectangles.

Quantization artifact in α-smoothing

As shown in Figure 6, α-smoothing on bi-modal distributions unintuitively outputs a “quantized” smoothed prediction. The prediction concentrates on a small, near-uniformly spaced set of quantization levels whose positions and weights are fully determined by the mixture and trimming parameters. 21

We describe this artifact with the simplest non-trivial case, a 1-D two-component Gaussian mixture, which is enough to expose the mechanism. The same reasoning extends mode-by-mode to mixtures with more components or to higher-dimensional outputs (where each axis is quantized independently after coordinate-wise sorting). F.1

Setup

Let the output distribution be the bi-modal Gaussian mixture p(y) = πA N (y | µA , σ 2 ) + πB N (y | µB , σ 2 ),

πA + πB = 1,

(40)

with µA < µB and noise scale σ small relative to the mode gap µB − µA . At a fixed smoothed input, i.i.d. the predictor draws N samples Y1 , . . . , YN ∼ p, sorts them as Y(1) ≤ · · · ≤ Y(N ) , and returns the symmetric α-trimmed mean µ̂α =

N −L 1 X Y(i) , M i=L+1

L = ⌊αN ⌋,

M = N − 2L.

(41)

Each Monte Carlo realization of the smoothed predictor is one independent draw of µ̂α . We call the set {Y(L+1) , . . . , Y(N −L) } the trimmed core. F.2

Quantization in the noiseless limit

Let K be the (random) number of samples drawn from mode A. By construction K ∼ Binom(N, πA ). In the noiseless limit σ → 0 the modes are perfectly separated, so after sorting the first K values come from mode A and the remaining N − K from mode B. The trimmed core then contains cA (K) = clip (K − L, 0, M ) values from mode A and cB (K) = M − cA (K) from mode B. The trimmed mean then reduces to cA (K)µA + cB (K)µB . (42) M Since cA takes only integer values 0, 1, . . . , M , the predictor collapses onto M + 1 discrete quantization levels jµA + (M − j)µB ℓj = , j = 0, 1, . . . , M, (43) M spaced uniformly between µA and µB at increments of (µA −µB )/M . Each level carries the binomial probability that K takes a value mapping to it,  j = 0,  Pr[K ≤ L] πj = Pr[K = L + j] j = 1, . . . , M − 1, (44)   Pr[K ≥ N − L] j = M. µ̂σ=0 α (K) =

The two boundary levels absorb the tails of Binom(N, πA ) because any K ≤ L saturates cA (K) = 0 and any K ≥ N − L saturates cA (K) = M . Addressing the failure mode of α smoothing. Increasing N and decreasing α help alleviate the behavior by enlarging the trimmed core, and as N → ∞ or α → 0, the distribution approaches a continuous distribution. However, both are problematic. Decreasing α means including more outliers, and increasing N collapses the distribution to a tight peak around the mean, which is even more disjoint from the individual modes. Our clustered α-smoothing avoids the artifact by summarizing each mode with its own trimmed mean and encoding the binomial split in the partition weights π (m) instead of in the level positions. F.3

Order-statistic correction at finite σ

The noiseless scenario highlights the origin of the quantization behavior, but neglects a bias in the level positions for σ > 0. The mode-A samples that survive trimming are not arbitrary mode-A draws but specific order statistics of the K mode-A samples (those closer to mode-B after sorting), and analogously for mode B. Both contributions are biased toward the gap between the modes. 22

Conditioning on K = k, the i-th overall order statistic is A(i:k) , d Y(i) | K = k = B(i−k: N −k) ,

i ≤ k, i > k,

(45)

where A(j:n) and B(j:n) denote the j-th order statistic among n i.i.d. draws from the corresponding mode. Writing νj:n = E[Z(j:n) ] for the expected j-th order statistic of n standard normals, the conditional expected trimmed mean is 1 [cA (k)µA + cB (k)µB + σ (SA (k) + SB (k))] , M with order-statistic correction terms E [µ̂α | K = k] =

SA (k) =

min(k,N −L)

(0 if k ≤ L),

νi:k

i=L+1

SB (k) =

N −L X i=max(L+1,k+1)

νi−k:N −k

(46)

(47)

(0 if k ≥ N − L).

(48)

Because the expected order statistic is monotone in its rank, SA (k) is non-negative (we are averaging the upper order statistics of the mode-A sub-sample whenever k > L) and SB (k) is non-positive (lower order statistics of mode B). Both shifts pull E[µ̂α | K = k] inward, i.e., toward the middle of the gap. The shifts are equal and opposite when cA (k) = cB (k) = M/2, so the central level ℓM/2 remains at (µA + µB )/2. The marginal expected level for the j-th quantization bin is the binomial weighted average of (46) over the k values mapping to that bin (a single k for interior bins, a tail of values for the two boundary bins). F.4

Worked example: N = 30, α = 0.4, πA = 0.4

With these parameters, the trimmed core is M = 6 samples, giving M + 1 = 7 quantization bins. The bins are nearly uniformly spaced; the order-statistic correction compresses the outer bins inward by roughly 0.22 ≈ 0.7σ and leaves the center bin invariant. Most of the probability mass concentrates at the boundary level closest to µB because E[K] = N πA = 12 = L, placing the median of the binomial split exactly at the lower trimming boundary. Figure 10 shows the underlying mixture (top) alongside the empirical histogram of µ̂α over 20,000 Monte Carlo realizations (bottom). The dashed lines mark the corrected levels ℓ̄j and labels report the binomial weights πj from (44). Both line up tightly with the empirical peaks, confirming (43)-(46).

Control of the residual in Proposition 3.2

One can achieve a desired residual level ϵ > 0 by considering a uniform partition of the interval [pV , p̄Vl ] into K > maxl∈L 1 + |L| ϵ (p̄Vl − pVl ) l P N that l∈L 1−Np̄V + p −p hl,K ≤ ϵ. l

N 1−p̄Vl + p

N −p

anchor points. Then, it holds

density

Underlying 2-component Gaussian mixture 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

πA = 0.40

πB = 0.60

Empirical α-trimmed mean (N = 30, α = 0.4, core 0.578 size = 6)

empirical (20000 simulations) quantization level

density

2.0 1.5 1.0 0.5 0.0

0.021 0.027 4

0.049

0.078

0.110

0.136

Figure 10: Quantization of the single-cluster α-trimmed predictor on a bimodal output. Top: the underlying mixture 0.4N (−3, 0.32 )+0.6 N (+3, 0.32 ). Bottom: histogram of µ̂α for N = 30, α = 0.4 over 20,000 realizations. Dashed lines are the predicted quantization levels ℓ̄j from (46); labels are the binomial weights πj from (44).

Clustered Randomized Smoothing for Stochastic Prediction Functions

Analyze & explore related research in Conceptio

Related documents