Year: 2021

Abstract:
Randomized Smoothing (RS) is a promising method for obtaining robustness certificates by evaluating a base model under noise. In this work, we: (i) theoretically motivate why ensembles are a particularly suitable choice as base models for RS, and (ii) empirically confirm this choice, obtaining state-of-the-art results in multiple settings. The key insight of our work is that the reduced variance of ensembles over the perturbations introduced in RS leads to significantly more consistent classifications for a given input. This, in turn, leads to substantially increased certifiable radii for samples close to the decision boundary. Additionally, we introduce key optimizations which enable an up to 55-fold decrease in sample complexity of RS for predetermined radii, thus drastically reducing its computational overhead. Experimentally, we show that ensembles of only 3 to 10 classifiers consistently improve on their strongest constituting model with respect to their average certified radius (ACR) by 5% to 21% on both CIFAR10 and ImageNet, achieving a new state-of-the-art ACR of 0.86 and 1.11, respectively. We release all code and models required to reproduce our results at ANONYMIZED.

heading: 1 INTRODUCTION
text: Modern deep neural networks are successfully applied to an ever-increasing range of applications, but while they often achieve excellent accuracy on the data distribution they were trained on, they have been shown to be very sensitive to slightly perturbed inputs, called adversarial examples (Biggio et al., 2013; Szegedy et al., 2014). This limits their applicability to safety-critical domains. Heuristic defenses against this vulnerability have been shown to be breakable (Tramèr et al., 2020; Carlini & Wagner, 2017), highlighting the need for provable robustness guarantees.
A promising method providing such guarantees for large networks is randomized smoothing (RS) (Cohen et al., 2019). The core idea is to provide probabilistic robustness guarantees with arbitrarily high confidence by adding noise to the input of a base classifier and computing the expected classification over the perturbed inputs using Monte Carlo sampling. The key to obtaining high robust accuracies is a base classifier with low variance with respect to these perturbations that remains consistently accurate even under high levels of noise. Existing works use different regularization and loss terms to encourage such behavior (Salman et al., 2019; Zhai et al., 2020; Jeong & Shin, 2020) but are ultimately all limited by the bias-variance trade-off of individual models. We show first theoretically and then empirically how ensembles can be constructed to significantly reduce this variance component and thereby increase certified radii for individual samples, and consequently, certified accuracy across the whole dataset. We illustrate this in Fig. 1.
Ensembles are a well-known tool for reducing classifier variance at the cost of increased computational cost (Hansen & Salamon, 1990). However, in the face of modern architectures (He et al., 2016; Huang et al., 2017) allowing the stable training of large models, they have been considered computationally inefficient. Yet, recent work (Wasay & Idreos, 2021) shows ensembles to be more efficient than single monolithic networks in many regimes. In light of this, we develop a theoretical framework analyzing this variance reducing property of ensembles under the perturbations introduced by RS. Further, we show how this reduced variance can significantly increase the majority class’s prediction probability, leading to much larger certified radii than evaluating more perturbations with an individual model.
Certification with RS is computationally costly as the base classifier has to be evaluated many thousand times. To avoid exacerbating these costs by using ensembles as base models, we develop
two techniques: (i) an adaptive sampling scheme for RS, which certifies samples for predetermined certification radii in stages, reducing the mean certification time up to 55-fold, and (ii) a special aggregation mechanism for ensembles which only evaluates the full ensemble on challenging samples, for which there is no consensus between a predefined subset of the constituting models.
Main Contributions Our key contributions are:
• A novel, theoretically motivated, and statistically sound soft-ensemble scheme for randomized smoothing, reducing perturbation variance and increasing certified radii (§4 and §5).
• A data-dependent adaptive sampling scheme for RS that reduces the sample complexity for predetermined certification radii in a statistically sound manner (§6).
• An extensive evaluation, examining the effects and interactions of ensemble size, training method, and perturbation size. We obtain state-of-the-art results on ImageNet and CIFAR10 for a wide range of settings, including denoised smoothing (§7).

heading: 2 RELATED WORK
text: Adversarial Robustness Following the discovery of adversarial examples (Biggio et al., 2013; Szegedy et al., 2014), defenses aiming to robustify networks were proposed (Madry et al., 2018). Particularly relevant to this work are approaches that certify or enforce robustness properties. We consider probabilistic and deterministic approaches.
Deterministic certification methods compute the reachable set for given input specifications using convex relaxations (Singh et al., 2019; Xu et al., 2020), mixed-integer linear programming (Tjeng et al., 2019), semidefinite programming (Dathathri et al., 2020), or SMT (Ehlers, 2017), to reason about properties of the output. To obtain networks amenable to such approaches, specialized training methods have been proposed (Mirman et al., 2018; Balunovic & Vechev, 2020; Xu et al., 2020).
Probabilistic certification (Li et al., 2019; Lécuyer et al., 2019) introduces noise to the classification process to obtain probabilistic robustness guarantees, allowing the certification of larger models than deterministic methods. We review Randomized Smoothing (RS) (Cohen et al., 2019) in §3 and associated training methods (Jeong & Shin, 2020; Zhai et al., 2020; Salman et al., 2019) in App. G.3. Orthogonally to training, RS has been extneded in numerous ways (Yang et al., 2020; Zhang et al., 2020; Lee et al., 2019; Dvijotham et al., 2020; Fischer et al., 2020), which we review in App. B.
Ensembles have been extensively analyzed with respect to different aggregation methods (Kittler et al., 1998; Inoue, 2019), diversification (Dietterich, 2000), and the reduction of generalization errors (Tumer & Ghosh, 1996b;a). Randomized Smoothing and ensembles were first combined in Liu et al. (2020) as ensembles of smoothed classifiers. However, the method does not retain strong certificates for individual inputs; thus, we consider the work to be in a different setting from ours. We discuss this in App. A. While similar at first glance, Qin et al. (2021) randomly sample models to an ensemble, evaluating them under noise to obtain an empirical defense against adversarial attacks. They, however, do not provide robustness guarantees.

heading: 3 RANDOMIZED SMOOTHING
text: Here we review the relevant background on Randomized Smoothing (RS) as introduced in Cohen et al. (2019). We let f : Rd 7→ Rm denote a base classifier that takes a d-dimensional input and produces m numerical scores (pre-softmax logits), one for each class. Further, we let F (x) := arg maxq fq(x) denote a function Rd 7→ [1, . . . ,m] that directly outputs the class with the highest score. For a random variable ∼ N (0, σ2 I), we define a smoothed classifier G : Rd 7→ [1, . . . ,m] as
G(x) := arg max c
P ∼N (0,σ2 I)(F (x+ ) = c). (1)
This classifier G is then robust to adversarial perturbations as follows:
Theorem 3.1 (From Cohen et al. (2019)). Let cA ∈ [1, . . . ,m], pA, pB ∈ [0, 1]. If
P (F (x+ ) = cA) ≥ pA ≥ pB ≥ max c6=cA P (F (x+ ) = c),
then G(x+ δ) = cA for all δ satisfying ‖δ‖2 < R with R := σ 2 (Φ−1(pA)− Φ−1(pB)).
Algorithm 1 Certify from (Cohen et al., 2019) 1: function CERTIFY(F, σ , x, n0, n, α) 2: cnts0← SAMPLEUNDERNOISE(F, x, n0, σ ) 3: ĉA ← top index in cnts0 4: cnts← SAMPLEUNDERNOISE(F, x, n, σ ) 5: pA ← LOWERCONFBND(cnts[ĉA], n, 1− α) 6: if pA > 12 then 7: return prediction ĉA and radius σ Φ−1(pA) 8: return For radii R as per Theorem 3.1 with inverse Gaussian CDF Φ−1 and any x, showingR ≥ δ certifies G(x) = G(x′) for any x′ ∈ {x+ η | ‖η‖2 < δ}. Computing the exact probabilities P (F (x+ ) = c) is generally intractable. Thus, to allow practical application, CERTIFY (Cohen et al., 2019) (see Algorithm 1) utilizes sampling: First, n0 samples to determine the majority class, then n samples to compute a lower bound pA to the success probability with confidence 1− α via the Clopper-Pearson lemma (Clopper & Pearson, 1934). If pA > 0.5, we set pB = 1− pA and obtain radius R = σ Φ−1(pA) via Theorem 3.1 else we abstain (return ). See App. F for exact definitions of the sampling and lower bounding procedures.
To obtain high certified radii, the base model F must be trained to cope with the added Gaussian noise . To achieve this, several training methods, discussed in App. G.3, have been introduced.
We also see this in Fig. 1, where various models obtain different pA, and thereby different radii R.

heading: 4 RANDOMIZED SMOOTHING FOR ENSEMBLE CLASSIFIERS
text: In this section, we extend the methods discussed in §3 from single models to ensembles.
For a set of k classifiers {f l : Rd 7→ Rm}kl=1, we construct an ensemble f̄ via weighted aggregation, f̄(x) = ∑k l=1 w
l γ(f l(x)), where wl are the weights, γ : Rm 7→ Rm is a post-processing function, and f l(x) are the pre-softmax outputs of an individual model. Soft-voting (where γ denotes identity) and equal weights wl = 1k perform experimentally well (see App. H.3.1) while being mathematically simple. Thus, we consider soft-ensembling via the averaging of the logits:
f̄(x) = 1
k k∑ l=1 f l(x) (2)
The ensemble f̄ and its corresponding hard-classifier F̄ (x) := arg maxq f̄q(x), can be used without further modification as base classifiers for RS. We find that classifiers of identical architecture and trained with the same method but different random seeds are sufficiently diverse too, when ensembled with k ∈ [3, 50], exhibit a notably reduced variance with respect to the perturbations . As we will show in the next section, this increases both the true majority class probability pA and its lower confidence bound pA, raising the certified radius as per Theorem 3.1.

heading: 5 VARIANCE REDUCTION VIA ENSEMBLES FOR RANDOMIZED SMOOTHING
text: We now show how ensembling even similar classifiers f l (cf. Eq. (2)) reduces their variance over the perturbations introduced in RS and thereby increases their majority class probability pA, making these ensembles f̄ particularly well-suited as base classifiers for RS. To this end, we first model network outputs with general distributions before investigating our theory empirically for Gaussian distributions. We defer algebra to App. C and justify our modeling assumptions empirically in App. D.
Our modeling approach is mathematically closely related to Tumer & Ghosh (1996a), which focuses on analyzing ensemble over the whole dataset. However, we condition on a single input to study the interplay between the stochasticity in training and random perturbations encountered in RS.
Individual Classifier We consider individual classifiers f l : Rd 7→ Rm and perturbed inputs x+ for a single arbitrary but fixed x and Gaussian perturbations ∼ N (0, σ2 I). We model the presoftmax logits f l(x) =: yl ∈ Rm as the sum of two random variables yl = ylp + ylc. Here, ylc corresponds to the classifier’s behavior on the unperturbed sample and models the stochasticity in weight initialization and training with random noise augmentation. ylp describes the effect of the random perturbations applied during RS. Note that this split will become essential when analyzing the ensembles. We drop the superscript l when discussing an individual classifier to avoid clutter.
We model the distribution of the clean component yc over classifiers with mean c = El[f l(x)], the expectation for a fixed sample x over the randomness in the training process, and covariance Σc characterizing this randomness. We assume the distribution of the perturbation effect yp to be zero-mean (following from local linearization and zero mean perturbations) and to have a covariance Σp. While Σp might depend on the noise level σ , it is distinct from it. We do not restrict the structure of either covariance matrix and denote Σii = σ2i and Σij = σiσjρij . As yc models the global training effects and yp models the local behavior under small perturbations, we assume them to be independent. We thus obtain logits y with mean c, the expectation over the training process for the clean sample, and covariance matrix Σ = Σc + Σp, where Σc and Σp encode the stochasticity of the training process and over perturbations respectively.
The classifier prediction F l(x) = arg maxq yq is not determined by the absolute values of logits but rather by the differences between them. We thus call the difference between the majority class logit and others the classification margin. During certification with RS, the first step is to determine the majority class. Without loss of generality, let the class with index 1 be the majority class, leading to the classification margin zi = y1 − yi. Note that if zi > 0 for all i 6= 1, then the majority class logit y1 is larger than the logits of all other classes yi. Under the above assumptions, the statistics of the classification margin for a single classifier are:
E[zi] = c1 − ci Var[zi] = σ 2 p,1 + σ 2 p,i + σ 2 c,1 + σ 2 c,i − 2ρp,1iσp,1σp,i − 2ρc,1iσc,1σc,i
Ensemble Now, we construct an ensemble of k of such classifiers. We use soft-voting (cf. Eq. (2)) to compute the ensemble output ȳ = 1k ∑k l=1 y
l and then the corresponding classification margins z̄i = ȳ1 − ȳi. By the linearity of expectation we have E[z̄i] = E[zi] = c1 − ci. We consider similar classifiers, differing only in the random seed used for training. Hence, we assume that the correlation between the logits of different classifiers has a similar structure but smaller magnitude than the correlation between logits of one classifier. Correspondingly, we parametrize the covariance between yic and y j c for classifiers i 6= j with ζcΣc and similarly between yip and yjp with ζpΣp for ζc, ζp ∈ [0, 1]. Depending on these correlation coefficients ζc and ζp, this model captures the range from no correlation (ζ = 0) to perfect correlation (ζ = 1). We obtain the following variance:
Var[z̄i] = k + 2
( k 2 ) ζp
k2 (σ2p,1 + σ 2 p,i − 2ρp,1iσp,1σp,i)︸ ︷︷ ︸
σ2p(k)
+ k + 2
( k 2 ) ζc
k2 (σ2c,1 + σ 2 c,i − 2ρc,1iσc,1σc,i)︸ ︷︷ ︸
σ2c(k)
.
Variance Reduction We can split Var[z̄i] into the components associated with the clean prediction σ2c (k) and the perturbation effect σ 2 p(k), both as functions of the ensemble size k. We now compare
these variance terms independently to the corresponding terms of an individual classifier dropping the subscripts p and c from σ2c (k)/σ 2 c (1) and σ 2 p(k)/σ 2 p(1) as they follow the same structure:
σ2(k) σ2(1) = (1 + ζ(k − 1))(σ21 + σ2i − 2ρ1iσ1σi) k(σ21 + σ 2 i − 2ρ1iσ1σi) = 1 + ζ(k − 1) k k→∞−−−−→ ζ (3)
We observe that both variance components go towards their corresponding correlation coefficients ζc and ζp as ensemble size grows, highlighting the importance of non-identical classifiers and illustrated under a Gaussian assumption (explained later) in Fig. 2a.
Effect on Success Probability We can compute the probability of an ensemble of k classifiers predicting the majority class by integrating the probability distribution of the classification margin over the orthant1 z̄i > 0 ∀i ≥ 2 where majority class 1 is predicted:
p1 := P(F̄ (x+ ) = 1) = P (z̄i > 0 : ∀ 2 ≤ i ≤ m) = ∫
z̄s.t.z̄i>0, ∀ 2 ≤i≤m
P(z̄) dz̄. (4)
Assuming Gaussian distributions, we observe the increase in success probability shown in Fig. 2b. Without assuming a specific distribution, we can still lower-bound the success probability using Chebyshev’s inequality and the union bound. Given that a mean classification yields the majority class and hence c1 − ci > 0, we let ti = c1−ciσi(k) and have:
p1 ≥ 1− m∑ i=2 P ( |z̄i − c1 + ci| ≥ ti σi(k) ) ≥ 1− m∑ i=2 σi(k) 2 (c1 − ci) 2 (5)
where σi(k)2 = σc,i(k)2 + σp,i(k)2 is the variance of classification margin z̄i. We observe that as σi(k) decreases with increasing ensemble size k (see Eq. (3)), the lower bound to the success probability increases quadratically. The further we are away from the decision boundary, i.e., the larger ci, the smaller the absolute increase in success probability for the same variance reduction.
Given a concrete success probability p1, we compute the probability distribution over the certifiable radii reported by CERTIFY for a given confidence α, sample number n, and perturbation variance σ2 (up to choosing an incorrect majority class ĉA) as:
P ( R = σ Φ −1(p1(n1, n, α))) = B(n1, n, p1), for R > 0 (6) where B(s, r, p) is the probability of drawing s successes in r trials from a Binomial distribution with success probability p, and p(s, r, α) is the lower confidence bound to the success probability of a Bernoulli experiment given s successes in r trials with confidence α according to the Clopper-Pearson interval (Clopper & Pearson, 1934). We illustrate the resulting effect assuming Gaussian distributions in Fig. 2c.
1n-dimensional equivalent of a quadrant
Empirical Analysis via Gaussian Assumption To investigate our theory, we now assume ylc and ylp and hence also z̄i to be multivariate Gaussians. This choice is empirically well-fitting (see App. D) and follows from the central limit theorem for ylc and from Gaussian perturbations and local linearization for ylp. To estimate the free parameters c, Σc, Σp, ζc, and ζp, we evaluate ensembles of up to k = 50 GAUSSIAN trained ResNet20 (for details, see §7) at σ = 0.25. We obtain c and Σc as the mean and covariance of the output on a randomly chosen sample x. Subtracting the clean outputs from those for the perturbed samples, we estimate the covariance matrix Σp. We determine ζc ≈ 0 and ζp ≈ 0.82 as the median ratio of the inter- and intra-classifier covariance. ζc ≈ 0 implies that our models can be treated as independent conditioned on a single fixed and unperturbed input.
Plugging these estimates into our model under the Gaussian assumption, we observe a significant decrease in the variance of the classification margin to the runner-up class as ensemble size k increases (see Fig. 2a). This generally leads to more polarized success probabilities: the majority class’ probability is increased, while the other classes’ probabilities are decreased because the probability mass concentrates around the mean, and consequently on one side of the decision threshold z̄ = 0, which determines the success probability via Eq. (4). An increase, as in our case (see Fig. 2b), leads to much larger expected certified radii (see Fig. 2c) for a given number of sampled perturbations (here n = 103) via Eq. (6). In contrast, sampling more perturbations will, in the limit, only recover the true success probability. In our example, going from one classifier to an ensemble of 50 increases the expected certified radius by 191%, while drawing 50 times more perturbations for a single
model only yields a 28% increase (see Fig. 2c). As illustrated in Fig. 3, these effects are strongest close to decision boundaries (pA 1), where small increases in pA impact the certified radius much more than the number of samples n, which dominated at pA ≈ 1 .̧

heading: 6 COMPUTATIONAL OVERHEAD REDUCTION
text: To avoid exacerbating the high computational cost of certification via RS with the need to evaluate multiple models constituting an ensemble, we propose two orthogonal and synergizing approaches to reduce its sample complexity significantly. The first is an adaptive sampling scheme reducing the average sample requirement for certification. The second is a more efficient aggregation mechanism for ensembles and is also applicable to inference.
Adaptive Sampling We assume a setting where a target certification radius r is known a priori, as is the standard in deterministic certification, and aim to show R ≥ r via Theorem 3.1. After choosing the majority class ĉA as in CERTIFY, CERTIFYADP (see Algorithm 2) leverages the insight that relatively few samples are often sufficient to certify robustness at radius r. We perform an s ≥ 1 step procedure: At step i, we evaluate ni fresh samples and, depending on these, either certify radius r, abort certification, or continue with ni+1 ni. After s steps, we abstain. To obtain the same strong statistical guarantees as with CERTIFY, it is essential to correct for multiple testing by applying Bonferroni correction (Bonferroni, 1936). The following theorem summarizes the key properties of the algorithm:
Theorem 6.1. For α, β ∈ [0, 1], s ∈ N+, n1 < · · · < ns, CERTIFYADP:
1. returns ĉA if at least 1− α confident that G is robust with a radius of at least r. 2. returns before stage s only if at least 1− β confident that robustness of G at radius r can
not be shown.
3. for ns≥dn(1−logα(s))e has maximum certifiable radii at least as large as CERTIFY for n.
We provide a proof in App. E. For well-chosen parameters, this allows certifying radius r with, in expectation, significantly fewer evaluations of F than CERTIFY as often ni n samples are sufficient to certify or to show that certification will not be successful. At the same time, we choose
Algorithm 2 Adaptive Sampling 1: function CERTIFYADP(F, σ ,x, n0, {nj}j , s, α, β, r) 2: cnts0← SAMPLEUNDERNOISE(F,x, n0, σ ) 3: ĉA ← top index in cnts0 4: for i← 1 to s do 5: cnts← SAMPLEUNDERNOISE(F,x, ni, σ ) 6: pA ← LOWERCONFBND(cnts[ĉA], ni, 1− αs ) 7: If σ Φ−1(pA) ≥ r then return ĉA 8: pA ← UPPERCONFBND(cnts[ĉA], ni, 1− βs−1 ) 9: If σ Φ−1(pA) < r then return 10: return
sample with n1 = 1000
sample with n2 = 10000
sample with n3 = 125000 ĉA
else
else
cnts[ĉA] < 795cnts[ĉA] ≥ 880
cnts[ĉA] < 8270cnts[ĉA] ≥ 8538
elsecnts[ĉA] ≥ 105607
Figure 4: CERTIFYADP with β = 0.0001 obtains similar guarantee as CERTIFY at n = 100′000 while being more sample efficient for most inputs. Both use α = 0.001, σ = r = 0.25.
ns in accordance with (3) to retain the ability to certify large radii (compared to CERTIFY) despite the increased confidence level required due to Bonferroni correction. For an example, see Fig. 4.
K-consensus Aggregation To reduce both inference and certification times, we can adapt the ensemble aggregation process to return an output early when there is consensus. Concretely, first, we order the classifiers in an ensemble by their accuracy (under noisy inputs) on a holdout dataset. Then, when evaluating f̄ , we query classifiers in this order. If the first K individual classifiers agree on the predicted class, we perform soft-voting on these and return the result without evaluating the remaining classifiers. Especially for large ensembles, this approach can significantly reduce inference time without hampering performance, as the ratio kK can be large.
We note that (i) similar approaches for partial ensemble evaluation have been proposed (Inoue, 2019; an, 2018; Soto et al., 2016), and (ii) that this does not affect the mathematical guarantees of RS.

heading: 7 EXPERIMENTAL EVALUATION
text: In line with prior work, we evaluate the proposed methods on the CIFAR10 (Krizhevsky et al., 2009) and ImageNet (Russakovsky et al., 2015) datasets with respect to two key metrics: (i) the certified accuracy at predetermined radii r and (ii) the average certified radius (ACR). We show that on both datasets, all ensembles outperform their strongest constituting model and obtain a new state-of-the-art. On CIFAR10, we demonstrate that ensembles outperform individual models even when correcting for computational cost via model size or the number of used perturbations. Further, we find that using adaptive sampling and K-Consensus aggregation speeds up certification up to 55-fold for ensembles and up to 33-fold for individual models.
Experimental Setup We implement our approach in PyTorch Paszke et al. (2019) and evaluate on CIFAR10 with ensembles of ResNet20 and ResNet110 and on ImageNet with ensembles of ResNet50 (He et al., 2016), using 1 and 2 GeForce RTX 2080 Ti, respectively. Due to the high computational cost of CERTIFY, we follow previous work (Cohen et al., 2019) and evaluate every 20th image of the CIFAR10 test set and every 100th of the ImageNet test set (500 samples total). We train models with GAUSSIAN (Cohen et al., 2019), CONSISTENCY (Jeong & Shin, 2020), and SMOOTHADV (Salman et al., 2019) training and utilize pre-trained SMOOTHADV and MACER (Zhai et al., 2020) models. More details are provided in App. G. If not declared differently, we use n0 = 100, n = 100′000, α = 0.001, no CERTIFYADP or K-consensus and evaluate and train with the same σ .

heading: 7.1 MAIN RESULTS
text: Results on CIFAR10 In Table 1, we compare ensembles of 10 ResNet110 against individual networks at σ = 0.50 w.r.t. average certified radius (ACR) and certified accuracy at various radii. We consistently observe that ensembles outperform their constituting models by up to 21% in ACR and 45% in certified accuracy, yielding a new state-of-the-art at every radius. Improvements are more pronounced for larger ensembles (see Fig. 5) and at larger radii, where larger pA are required, which agrees well with our theoretical observations in §5. We present more extensive results in App. H.1, including experiments for σ = 0.25 and σ = 1.0 yielding new state-of-the-art ACR in all settings.
Results on ImageNet In Table 2, we compare ensembles of 3 ResNet50 with individual models at σ = 1.0, w.r.t. ACR and certified accuracy. We observe similar trends to CIFAR10, with ensembles outperforming their constituting models and the ensemble of 3 CONSISTENCY trained ResNet50 obtaining a new state-of-the-art ACR of 1.11, implying that ensembles are effective on a wide range of datasets. We visualize this in Fig. 6 and present more detailed results in App. H.2, where we also provide experiments for σ = 0.25 and σ = 0.5, yielding new state-of-the-art ACR in both cases.
Computational Overhead Reduction We evaluate CERTIFYADP in conjunction with KConsensus aggregation using β = 0.001, {nj} = {100, 1′000, 10′000, 120′000}, and K = 5 in Table 3 and report KCR, the percentage of inputs for which only K classifiers were evaluated and SampleRF and TimeRF, the factors by which sample complexity and certification time are reduced, respectively. We observe up to 67-fold certification speed ups and up to 55-fold sample complexity
reductions for individual predetermined radii, without incurring any accuracy penalties. This way, certifying an ensemble of k = 10 networks with CERTIFYADP can be significantly faster than certifying an individual model with CERTIFY while yielding notably higher certified accuracies. We observe that adaptive sampling and K-Consensus aggregation complement each other well: At larger radii r, more samples pass on to the later stages of CERTIFYADP, requiring more model evaluations. However, only samples with high success probabilities reach these later stages, making it more likely for K-Consensus aggregation to terminate an evaluation early. We provide extensive additional experiments in App. H.4.

heading: 7.2 ABLATION STUDY
text: Ensemble Size and Model Size Ablation In Fig. 5, we illustrate the effect of ensemble size on ACR for various training methods and model architectures, providing more extensive results in App. H.3.2. We observe that in all considered settings, even relatively small ensembles of just 3-5 models realize most of the improvements to be made. In fact, as few as 3 and 5 ResNet20 are enough to outperform a single ResNet110 under GAUSSIAN and CONSISTENCY training, respectively, while only having 47% and 79% of the parameters.
Equal Number of Inferences In Fig. 7, we compare the effect of evaluating a larger number of perturbations n for an individual model compared to an ensemble at the same computational cost. We observe that increasing the sample count for a single model only enables marginal improvements at most radii and only outperforms the ensemble at radii, which can mathematically not be certified using the smaller sample count. See also Fig. 3 for a theoretical explanation and App. H.3.4 for more experiments.
Other Base Classifiers While our approach works exceptionally well when ensembling similar classifiers, it can also be successfully extended to base models trained with different training methods (see App. H.3.5) and denoised classifiers (see App. H.3.6), yielding ACR improvements of up to 5% and 17% over their strongest constituting models, respectively. This way existing models can be reused to avoid the high cost of training an ensemble.
Adaptive Sampling Ablation In Table 4, we show the effectiveness of CERTIFYADP for certifying robustness at predetermined radii using β = 0.001, different target radii, training methods, and sampling count sets {nj}. We compare ASRj , the portion of samples returned in stage j, SampleRF, and TimeRF. Generally, we observe larger speed-ups for smaller radii (up to a factor of 34 at r = 0.25). There, certification only requires small pA, which can be obtained even with few samples. For large radii (w.r.t σ ), higher pA must be shown, requiring more samples even at true success rates of pA = 1. There, the first phase only yields early abstentions. The more phases we use, the more often we certify or abstain early, but also the bigger the costs due to Bonferroni correction. We observe that exponentially increasing phase sizes and in particular {nj} = {100, 1′000, 10′000, 120′000} perform very well across different settings. See App. H.4.1 for additional results.
K-Consensus Aggregation Ablation In Table 5, we compare ACR and certification time reduction for ensembles of 10 ResNet110 and 50 ResNet20 using K-Consensus aggregation in isolation across different K. Even when using only K = 2, we already obtain 81% and 70% of the ACR improvement obtainable by always evaluating the full ensembles (k = 10 and k = 50), respectively. The more conservative K = 10 for ResNet20 still reduces certification times by a factor of 2 without losing accuracy. See App. H.4.2 for more detailed experiments.

heading: 8 CONCLUSION
text: We propose a theoretically motivated and statistically sound approach to construct low variance base classifiers for Randomized Smoothing by ensembling. We show theoretically and empirically why this approach significantly increases certified accuracy resulting in state-of-the-art results. To offset the computational overhead of ensembles, we develop a generally applicable adaptive sampling mechanism, reducing certification costs up to 55-fold for predetermined radii and an ensemble aggregation mechanism, complementing it, and reducing evaluation costs on its own up to 6-fold.

heading: 9 ETHICS STATEMENT
text: Most machine learning techniques can be applied both in ethical and unethical ways. Techniques to increase the certified robustness of models do not change this but only make the underlying models more robust, again allowing beneficial and malicious applications, e.g., more robust medical models versus more robust weaponized AI. As our contributions improve certified accuracy, certification radii, and inference speed, both of these facets of robust AI are amplified. Furthermore, while we achieve state-of-the-art results, these do not yet extend to realistic perturbations. Malicious actors may aim to convince regulators that methods such as the proposed approach are sufficient to provide guarantees for general real-world threat scenarios, leading to insufficient safeguards.

heading: 10 REPRODUCIBILITY STATEMENT
text: We will make all code and pre-trained models required to reproduce our results publicly available at ANONYMIZED upon publication. Additionally, ANONYMIZED also contains detailed instructions and scripts facilitating the reproduction of our results. We explain the basic experimental setup in §7 and provide exact experimental details in App. G, where we extensively describe the datasets (App. G.1), model architectures (App. G.2), training methods (App. G.3), hyper-parameter choices (App. G.4), and experiment timings (App. G.5). We remark that the datasets, architectures, and training methods we use are all publicly available. Section §7 and App. H contain numerous experiments for various datasets, model architectures, training methods, noise levels, and ensemble compositions. The consistent trends we observe over this wide range of settings highlight the general applicability of our approach. In App. H.3.3, we analyze the variability of our results, reporting standard deviations for a range of metrics, find that our ensembling approach reduces variability, and note that our results are statistically significant. We include complete proofs of all theoretical contributions (§5 and §6) in App. C and E. In App. D, theoretically motivate the modeling assumptions made in §5 and validate them empirically.

heading: A COMPARISON TO LIU ET AL. (2020)
text: Liu et al. (2020) introduce Smoothed Weighted Ensembling (SWEEN), which are weighted ensembles of smoothed base models. They focus on deriving optimal ensemble weights for generalization rather than on analyzing the particularities of ensembles in the RS setting. As they do not consider multiple testing correction and the confidence of the applied Monte Carlo sampling, they obtain statements about empirical, distributional robustness of their classifiers, rather than individual certificates with high confidence. Due to this difference setting, we do not compare numerically.

heading: B ADDITIONAL RELATED WORK
text: Extensions As outlined in §3, Cohen et al. (2019) focuses on the addition of Gaussian noise to drive `2-robustness results. However, various extensions with other types of noise and guarantees have been proposed.
By using different types of noise distributions and radius calculations, Yang et al. (2020); Zhang et al. (2020) derive recipes to determine certificates for general `p-balls and specifically showcase results for p = 1, 2,∞. Numerous works have shown extensions to discrete perturbations such as `0-perturbations (Lee et al., 2019; Wang et al., 2021; Schuchardt et al., 2021), graphs (Bojchevski et al., 2020; Gao et al., 2020), patches Levine & Feizi (2020) or manipulations of points in a point cloud Liu et al. (2021).
Dvijotham et al. (2020) provide theoretical derivations for the application of both continuous and discrete smoothing measures.
Mohapatra et al. (2020) improve certificates by considering gradient information.
Beyond norm-balls certificates, Fischer et al. (2020); Li et al. (2021) show how geometric operations such as rotation or translation can be certified via randomized smoothing.
Chiang et al. (2020); Fischer et al. (2021) show how the certificates can be extended from the setting of classification to regression (and object detection) and segmentation, respectively. In the classification setting, Jia et al. (2020) extend certificates from just the top-1 class to the top-k classes. Similarly, Kumar et al. (2020a) certify the confidence of the classifier not, just the top class prediction.
Beyond the inference-time evasion attacks, Rosenfeld et al. (2020) showcase RS as a defense against data poisoning attacks.
Many of these extensions are orthogonal to the standard randomized smoothing approach and thus orthogonal to our improvements. We showcase that these improvements carry over from the `2 case to other `p norms in App. H.3.7.
Limitations Recently, some works also have investigated the limitations of randomized smoothing: Mohapatra et al. (2021) point out the limitations of training for randomized smoothing as well as the impact on class-wise accuracy. Kumar et al. (2020b); Wu et al. (2021) show that for p > 2, the achievable certification radius quickly diminishes in high-dimensional spaces.

heading: C MATHEMATICAL DERIVATION OF VARIANCE REDUCTION FOR ENSEMBLES
text: In this section, we present the algebraic derivations for §5, skipped in the main part due to space constraints.
Individual Classifier In §5 define f l(x) =: yl ∈ Rm as the sum of two random variables yl = ylp + y l c.
The behavior on a specific clean sample x is modeled by ylc with mean c ∈ Rm, the expectation of the logits for this sample over the randomization in the training process, and corresponding covariance
Σc:
E[ylc] = c
Cov[ylc] = Σc =  σ 2 c,1 · · · ρc,1mσc,1σc,m . . . ...
σ2c,m  We note that: (i) we omit the lower triangular part of covariance matrices due to symmetry, (ii) c and Σc are constant across the different f l for a fixed x and (iii) that due to (ii) and our modelling assumptions all means, (co)variances and probabilities are conditioned on x, e.g., E[ylc] = E[ylc | x] = c. However, as we define all random variables only for a fixed x, we omit this.
Similarly, we model the impact of the perturbations introduced during RS by ylp with mean zero and covariance Σp:
E[ylp] = 0
Cov[ylp] = Σp =  σ 2 p,1 · · · ρp,1mσp,1σp,m . . . ...
σ2p,m  That is, we assume cl to be the expected classification a model learns for the clean sample while the covariances Σc and Σp encode the stochasticity of the training process and perturbations, respectively. As ylp models the local behavior under small perturbations and y l c models the global training effects, we assume them to be independent:
E[yl] = c
Cov[yl] = Σ = Σc + Σp =  σ 2 p,1 + σ 2 c,1 · · · ρp,1mσp,1σp,m + ρc,1mσc,1σc,m . . . ...
σ2p,m + σ 2 c,m  . The classifier prediction arg maxq yq is determined by the differences between logits. We call the difference between the target logit and others the classification margin. During certification with RS, the first step is to determine the majority class. Without loss of generality, we assume that it has been determined to be index 1, leading to the classification margin zi = y1 − yi. If zi > 0 for all i 6= 1, the majority class logit y1 is larger than those of all other classes yi. We define z := [z2, .., zn]T , skipping the margin of y1 to itself. Under the above assumptions, the statistics of the classification margin for a single classifier are:
E[zi] = c1 − ci Var[zi] = σ 2 p,1 + σ 2 p,i + σ 2 c,1 + σ 2 c,i − 2ρp,1iσp,1σp,i − 2ρc,1iσc,1σc,i
Ensemble Now, we construct an ensemble of k of these classifiers. Using soft-voting (cf. Eq. (2)) to compute the ensemble output ȳ = 1k ∑k l=1 y l. We assume the yip and y j p to be correlated with ζpΣp for classifiers i 6= j and similarly model the correlation of ylc with ζcΣc for ζp, ζc ∈ [0, 1]. Letting y∗ := [y1 > , ...,yk > ]> denote the concatenation of the logit vectors of all classifiers, we assemble their joint covariance matrix in a block-wise manner from the classifier individual covariance matrices as:
Cov[y∗] = Σ∗ =  Σp + Σc · · · ζpΣp + ζcΣc. . . ... Σp + Σc  We then write the corresponding classification margins z̄i = ȳ1 − ȳi. or in vector notation z̄ := [z̄2, .., z̄m]
T , again skipping the margin of ȳ1 to itself. By linearity of expectation we obtain E[z̄i] = E[zi] = c1 − ci or equivalently
E[z̄] = µ̄ = ( c1 ... c1 ) − c[2:m].
We define the ensemble difference matrixD ∈ Rm−1×mk with elements di,j such that z̄ = Dy∗:
Dij =  1 k , if jmodm = 1 −1 k , else if jmodm = imodm
0, else , j ∈ [1, ...,mk], i ∈ [2, ...,m].
This allows us to write the covariance matrix of the ensemble classification margins
Cov[z̄] = Σ̄ = DΣ∗D>.
Now we can evaluate its diagonal elements or use the multinomial theorem and rule on the variance of correlated sums to obtain the variance of individual terms:
Var[z̄i] = k + 2
( k 2 ) ζp
k2 (σ2p,1 + σ 2 p,i − 2ρp,1iσp,1σp,i)︸ ︷︷ ︸
σ2p(k)
+ k + 2
( k 2 ) ζc
k2 (σ2c,1 + σ 2 c,i − 2ρc,1iσc,1σc,i)︸ ︷︷ ︸
σ2c(k)
.
Variance Reduction We can split Var[z̄i] into the components associated with the perturbation effect σ2p(k) and the clean prediction σ 2 c (k), all as functions of the ensemble element number k.
Now, we analyze these variance components independently by normalizing them with the corresponding components of an individual classifier:
σ2p(k) σ2p(1) = (1 + ζp(k − 1))(σ2p,1 + σ2p,i − 2ρp,1iσp,1σp,i) k(σ2p,1 + σ 2 p,i − 2ρp,1iσp,1σp,i) = 1 + ζp(k − 1) k k→∞−−−−→ ζp
σ2c (k) σ2c (1) = (1 + ζc(k − 1))(σ2c,1 + σ2c,i − 2ρc,1iσc,1σc,i) k(σ2c,1 + σ 2 c,i − 2ρc,1iσc,1σc,i) = 1 + ζc(k − 1) k k→∞−−−−→ ζc
We observe that both variance components go towards their corresponding correlation coefficients ζp and ζc as ensemble size grows, highlighting the importance of non-identical classifiers.
Especially for samples that are near a decision boundary, this variance reduction will lead to much more consistent predictions, in turn significantly increasing the lower confidence bound p
1 and
thereby the certified radius as per Theorem 3.1.

heading: D THEORY VALIDATION
text: In this section, we present experiments validating the modeling assumptions made in §5 to derive the reduced variance over perturbations in RS resulting from using ensembles of similar classifiers. There are three main assumptions we make:
• The (pre-softmax) output of a single model can be modeled as f l(x) = yp + yc, where yc captures the prediction on clean samples and yp the effect of perturbations with mean zero.
• The clean and perturbation components yc and yp are independent.
• The covariance of ylp and y l c between different classifiers i 6= j can be parametrized as
Cov(yip,y j p) = ζpΣp and similarly Cov(y i c,y j c) = ζcΣc.
To show that these assumptions are valid, we additionally assume multivariate Gaussian distributions and evaluate an ensemble of k = 50 ResNet20 on 100 clean samples x, each perturbed with n = 1000 different error terms , and compare the obtained distributions and covariance matrices with our modeling. We generally observe excellent agreement between observation and model at low noise level σ = 0.25 and good agreement at a high noise level of σ = 1.0.
Individual classifier predictions We model the predictions of individual classifiers f : Rd 7→ Rm on perturbed inputs x+ for a single arbitrary but fixed x and Gaussian perturbations ∼ N (0, σ2 I) with f l(x) = yp + yc.
We can interpret the clean component yc in two different ways: either as the prediction yc := f(x) on the clean sample x or as the mean prediction over perturbations of this sample yc := E ∼N(0,σ2 I)[f(x + )]. In both cases, we assume yc to be distributed with mean c ∈ Rm and covariance Σc. The first case is illustrated in Fig. 8, where we rescale both our model and the predictions with the sample-wise mean and variance over classifiers and observe excellent agreement between model and observed distribution. The second case is illustrated in the leftmost column in Fig. 9 for two different σ , where we again observe excellent agreement between the true distribution and our model. Here we denote with µf(x′) the mean prediction over perturbed samples E ∼N(0,σ2 I)[f(x+ )], corresponding to the setting of yc := E ∼N(0,σ2 I)[f(x+ )] and similarly with denote with µf(x) the mean prediction over classifiers for an unperturbed sample El[f l(x)], corresponding to the setting of yc := f(x).
The perturbation effect yp is assumed to be mean zero, either following from local linearization and zero mean perturbations when defining the clean component as yc := f(x) or directly from the definition of the clean component as yc := E ∼N(0,σ )[f(x + )]. For small perturbations (σ = 0.25), the assumption of local linearity holds, and we observe a mean perturbation effect distributed very tightly around 0 (see top right pane in Fig. 9). For larger perturbations (σ = 1.00), this assumption of local linearity begins to break down, and we observe a much wider spread of mean perturbation effect (see bottom right pane in Fig. 9). At both perturbation levels, we observe an excellent agreement of the perturbation effect model with the observed data when using the definition of yc := E ∼N(0,σ )[f(x+ )] (green in the middle column in Fig. 9). Using yc := f(x), we still observe great agreement at the lower perturbation magnitude but a notable disagreement at high perturbation levels. While both definitions for yc are fully compatible with our derivation in §5, we choose to present a version based on yc := f(x) as it allows for a better intuitive understanding.
Correlation between clean and error component We model yc and yp separately to distinguish between randomness over perturbations, which are introduced during RS, and over the models introduced during training with different random seeds. As yc models the large scale effects of training and yp the smaller local scale local effects of inference time perturbations, we assume them to be independent, corresponding to a covariance Cov(yp,i, yc,j) ≈ 0 ∀i, j. Depending on which definition of yc we use, we obtain either of the two covariance matrices in Fig. 10, where the upper left-hand block is the covariance of yc, the lower righthand block the covariance of yp and the off-diagonal blocks Cov(yp,i, yc,j). In the setting of yc := µf (x′), the two components yc and yp are perfectly independent as can be seen by the 0 off-diagonal blocks in the covariance matrix shown in Fig. 10a. In the setting of yc := f(x), very small terms can be seen in the off-diagonal blocks shown in Fig. 10b. Overall, the agreement with our model is excellent in both cases.
Inter classifier correlation structure In our analysis, we consider ensembles of similar classifiers which are based on the same architecture and trained using the same methods, just using different random seeds. Hence, we assume that the correlation between the logits of different classifiers has a similar structure but smaller magnitude than the correlation between the logits of the same classifier. Correspondingly, we parametrize the covariance between yip and y j p for classifiers i 6= j with ζpΣp and similarly between yic and y j c with ζcΣc for ζp, ζc ∈ [0, 1]. To evaluate this assumption, we
consider all (
50 2
) pairings of classifiers and compute the resulting pairwise covariance matrices for
the clean components ylc := f l(x) and perturbation effects ylp. In Fig. 11 we show the resulting covariance matrices (left column), our model (right column), and the delta between the two (center column). Again we observe excellent agreement to the point that even given the delta maps, it is hard to spot the differences in the covariance matrices.
Model fit conclusion We conclude that our model is able to capture the relevant behavior even of deep networks (here ResNet20) with great accuracy that goes beyond what is required for a theoretical analysis to allow for interesting insights into the underlying mechanics.

heading: E PROOF OF THEOREM 6.1
text: We first restate the theorem:
Theorem (Theorem 6.1 restated). For α, β ∈ [0, 1], s ∈ N+, n1 < · · · < ns, CERTIFYADP:
1. returns ĉA if at least 1− α confident that G is robust with a radius of at least r. 2. returns before stage s only if at least 1− β confident that robustness of G at radius r can
not be shown.
3. for ns≥dn(1−logα(s))e has maximum certifiable radii at least as large as CERTIFY for n.
Proof. We show statements 1, 2, and 3 individually.
Proof of 1. If a phase in CERTIFYADP returns ĉA, then via Theorem 3.1, G will be not robust with probability at most α/s. Via Bonferroni correction (Bonferroni, 1936), if CERTIFYADP returns ĉA, then G is robust with radius at least r with confidence 1− s(α/s) = 1− α. Proof of 2. If CERTIFYADP returns in phase s, we have pA < p′A = Φ(r/σ ) the minimum success probability for certification at r. By the definition of the UPPERCONFIDENCEBOUND the true success probability of G will be pA ≤ pA with confidence at least β/(s− 1). Hence, if phase j returns , robustness of G cannot be shown at radius r with confidence 1− β/(s− 1), even with access to the true success probability or infinitely many samples. Again with Bonferroni correction (Bonferroni, 1936), the overall probability that G is robust at r despite CERTIFYADP abstaining early is at most (s− 1)(β/(s− 1)) = β. Proof of 3. Finally, to prove the last part, we assume cnts[ĉA] = n, yielding the largest certifiable radii for any given n and α. Now, the largest radius provable via Theorem 3.1 with n samples at α is
α1/n (Cohen et al., 2019). Similarly, for ns samples at α/s, it is ( α s ) 1 ns . Then we have the following equivalences:
(α) 1 n = (α s ) 1 ns ⇔ αns = (α s )n ⇔ ns logα = n(logα− log s)⇔ ns = n(1− n logα(s))
Hence, if we choose ns = dn(1− logα(s))e, then we can certify at least the same maximum radius with ns samples at α/s as with n samples at α. Thus, overall, we can certify the same maximum radius at α.

heading: F ADDITIONAL RANDOMIZED SMOOTHING DETAILS
text: This section contains definitions needed for the standard certification via randomized smoothing, CERTIFY, proposed by Cohen et al. (2019).
SAMPLEUNDERNOISE(F, x, n, σ ) first samples n inputs x1, . . . , xn with x1 = x+ 1, . . . , xn = x + n for 1, . . . , n ∼ N (0, σ ). Then it counts how often F predicts which class for these x1, . . . , xn and returns the corresponding m dimensional array.
UPPERCONFBND(k, n, 1 − α) returns an upper bound on p with confidence at least 1 − α where p is an unknown probability such that k ∼ B(n, p). Here, B(n, p) is a binomial distribution with parameters n and p.
LOWERCONFBND(k, n, 1 − α) is analogous to UPPERCONFBND but returns a lower bound on p instead of an upper bound with confidence 1− α. We point the interested reader to Cohen et al. (2019) for more details about randomized smoothing.

heading: G EXPERIMENTAL DETAILS
text: In this section, we provide greater detail on the datasets, architectures, and training and evaluation methods used for our experiments.

heading: G.1 DATASET DETAILS
text: We use the MNIST, CIFAR10, and ImageNet datasets in our experiments.
MNIST (Deng, 2012) contains 60’000 training and 10’000 test set images partitioned into 10 classes for the 10 digits. We evaluate all images from the test set.
CIFAR10 (Krizhevsky et al., 2009) (MIT License) contains 50’000 training and 10’000 test set images partitioned into 10 classes. We evaluate every 20th image of the test set starting with index 0 (0, 20, 40, ..., 9980) in our experiments, following previous work (Cohen et al., 2019).
ImageNet (Russakovsky et al., 2015) contains 1’287’167 training and 50’000 validation images, partitioned into 1000 classes. We evaluate every 100th image of the validation set starting with index 0 (0, 100, 200, ..., 49900) in our experiments, following previous work (Cohen et al., 2019).
G.2 ARCHITECTURE DETAILS
We use different versions of ResNet (He et al., 2016) for our experiments. Concretely, we evaluate ensembles of ResNet20 and ResNet110 on CIFAR10 and ensembles of ResNet50 on ImageNet.
ResNet110 has about 6.35 times as many parameters as ResNet20 (see Table 6). ResNet50 instantiated for ImageNet has substantially more parameters than ResNet110 instantiated for CIFAR10 because of the significantly larger input dimension of ImageNet samples.

heading: G.3 TRAINING METHODS
text: To obtain high certified radii via RS, the base model F must be trained to cope with the added Gaussian noise . To this end, Cohen et al. (2019) propose data augmentation with Gaussian noise during training, referred to as GAUSSIAN in the following. Building on this, Salman et al. (2019) suggest SMOOTHADV, a combination of adversarial training (Madry et al., 2018; Kurakin et al., 2017; Rony et al., 2019) with the data augmentation from GAUSSIAN. (Here, we always consider the PGD version.) While improving accuracy, this training procedure is computationally very expensive. MACER (Zhai et al., 2020) achieves similar performance with a cheaper training procedure, adding a loss term directly optimizing a surrogate certification radius. Jeong & Shin (2020), called CONSISTENCY in the following, replace this term with a more easily optimizable loss, further decreasing training time and improving performance. Depending on the setting, the current state-of-the-art results are either achieved by SMOOTHADV, MACER, CONSISTENCY, or a combination thereof.
Further, Salman et al. (2020) present denoised smoothing, where the base classifier f = h ◦ f ′ is the composition of a denoiser h that removes the Gaussian noise from the input and an underlying classifier f ′ that is not specially adapted to Gaussian noise, e.g., accessed via an API.

heading: G.4 TRAINING DETAILS
text: All models are implemented in PyTorch (Paszke et al., 2019) (customized BSD license).
For GAUSSIAN training (No license specified), we use the published code by Cohen et al. (2019). For each network, we chose a different random seed, otherwise using identical parameters. Overall, we trained 50 ResNet20 and 10 ResNet110 for each σ ∈ {0.25, 0.5, 1.0} on CIFAR10. For ImageNet, we trained 3 ResNet50 for σ = 1.0, and took the results for σ ∈ {0.25, 0.50} from their GitHub. For CONSISTENCY training (MIT License), we use the code published by Jeong & Shin (2020) and CONSISTENCY instantiation built on GAUSSIAN training. Similarly to GAUSSIAN, we use a different random seed and otherwise identical parameters for all networks. We chose the parameters reported to yield the largest ACR by Jeong & Shin (2020) for any given σ , except for σ = 0.25 on ImageNet, for which no parameters were reported. In detail, for CIFAR10 we generally use η = 0.5, using λ = 20 for σ = 0.25 and λ = 10 for σ = 0.5 and σ = 1.0 (all for both ResNet20 and ResNet110). In this way, we train 50 ResNet20 and 10 ResNet110 for each σ ∈ {0.25, 0.5, 1.0}. For ImageNet, we use η = 0.1 and λ = 5 for σ = 1.0. For σ = 0.25 and σ = 0.5, we use η = 0.5 with λ = 10 respectively λ = 5.
For SMOOTHADV training (MIT License), we use the PGD based instantiation of the code published by Salman et al. (2019). Note that due to the long training times, we only train SMOOTHADV models ourselves on CIFAR10. The individual models we train only differ by random seed, all using PGD attacks. For σ = 0.25, we train with T = 10 steps, an = 255/255, 10 epochs of warm-up, and mtrain = 8 noise terms per sample during training. For σ = 0.5, we train with T = 2 steps, an = 512/255, 10 epochs of warm-up, and mtrain = 2 noise terms per sample during training. Finally, σ = 1.0, we train with T = 10 steps, an = 512/255, 10 epochs of warm-up, and mtrain = 2 noise terms per sample during training.
We use the same training schedule and optimizer for all models, i.e., stochastic gradient descent with Nesterov momentum (weight = 0.9, no dampening), with an `2 weight decay of 0.0001. For CIFAR10, we use a batch size of 256 and an initial learning rate of 0.1, reducing it by a factor of 10 every 50 epochs and training for a total of 150 epochs. For ImageNet, we use the same settings, only reducing the total epoch number to 90 and decreasing the learning rate every 30 epochs.
All single MACER (No license specified) trained models for CIFAR10 are taken directly from Zhai et al. (2020).
For our experiments on denoised smoothing, we use the 4 white box denoisers with DNCNN-WIDE architecture trained with learning schedules 1, 3, 4, and 5 for STAB ResNet110, and the ResNet110 trained on unperturbed samples for 90 epochs from Salman et al. (2020) (MIT License).
To rank the single models for CIFAR10, we have evaluated them on a disjunct hold-out portion of the CIFAR10 test set. Concretely, we use the test images with indices 1, 21, 41, ..., 9981 to rank the single models for GAUSSIAN, CONSISTENCY, and SMOOTHADV trained models. The performances of individual models values we report are those of the models with the best score on this validation set for CIFAR10 and the best score on the test set for ImageNet (favoring individual models in this setting). For ensembles of size k < 10 and k < 50 for ResNet110 and ResNet20, respectively, we ensemble the k models according to their performance on the hold-out set. We note that other combinations might yield stronger ensembles, but an exhaustive search of all combinatorially many possibilities is computationally infeasible.
CIFAR10 models were trained on single GeForce RTX 2080 Ti and ImageNet models on quadruple 2080 Tis. We report the epoch-wise and total training times for individual models in Table 7 (when trained sequentially one at a time).

heading: G.5 EXPERIMENT TIMINGS
text: To evaluate the ensembles of ResNet20 and ResNet110 on CIFAR10, we use single GeForce RTX 2080 Ti, and to evaluate ensembles of ResNet50 on ImageNet, we use double 2080 Tis.
For both datasets, 500 samples are evaluated per experiment as discussed in App. G.1. Below we list the time required for certification using CERTIFY and no K-Consensus aggregation:
• ResNet110 on CIFAR10 (1 GPU): 4.12h per single model; 41.2h per ensemble of 10
• ResNet20 on CIFAR10 (1 GPU): 0.825h per single model; 41.3h per ensemble of 50
• ResNet50 on ImageNet (2 GPUs): 8.75h per single model; 26.2h per ensemble of 3
The timing for different ensemble sizes scales linearly between full size and individual model timings. The time required for certification using CERTIFYADP or K-Consensus aggregation or both can be obtained by dividing the timings reported above with the speed-up factor TimeRF reported for the corresponding experiments.

heading: H ADDITIONAL EXPERIMENTS
text: In the following, we present numerous additional experiments and more detailed results for some of the experiments presented in §7.

heading: H.1 ADDITIONAL RESULTS ON CIFAR10
text: In this section, we present more experimental results on CIFAR10. Concretely, Table 8 presents the performance of ensembles of ResNet20 (k ∈ {1, 5, 50}) and ResNet110 (k ∈ {1, 10}) trained using a wide range of of methods for σ = 0.25, σ = 0.5, and σ = 1.0. We again consistently observe that ensembles significantly outperform their constituting models, implying that ensembles are effective for various training methods and model architectures, as well as different σ . In particular, they lead to a new state-of-the-art both in ACR and at most radii at σ = 0.25, σ = 0.5, and σ = 1.0.
Fig. 12 visualizes the evolution of ACR and certified accuracy with ensemble size for σ = 0.25. In Fig. 13, we visualize the evolution of certified accuracy over radii for GAUSSIAN trained models.

heading: H.2 ADDITIONAL RESULTS ON IMAGENET
text: Extending the results from §7, Table 9 also provides ensemble results on ImageNet for σ = 0.25, and σ = 0.5 with CONSISTENCY trained ResNet50. Similarly to CIFAR10, we consistently observe that ensembles significantly outperform their constituting models, implying that ensembles are effective on various datasets. In particular, we achieve a new state-of-the-art both in ACR and certified accuracy at most radii for σ = 0.25, σ = 0.50, and σ = 1.00.
Additionally, we report the performance of all individual classifiers which constitute the ensembles, where k = 3 combines all three and k = 2 the first two for GAUSSIAN and CONSISTENCY, and the last two for SMOOTHADV. Note that the SMOOTHADV models are taken directly from Salman et al. (2019) and were trained with a different parameter (256, 512, and 1024). We observe that if there is a large discrepancy between model performance at some radii (e.g., the SMOOTHADV models at
small radii), the ensemble might perform worse than the strongest constituting model. Overall and when performance is homogenous, however, ensembles once again obtain significant performance increases over individual models.

heading: H.3 ADDITIONAL ABLATION STUDY OF ENSEMBLES
text: This section presents a range of experiments on different different aspects of variance reduction through ensembles for randomized smoothing. In particular, we cover aggregation approaches, the effect of ensemble size, the variability of our results, certifications with an equal number of inferences, ensembles of differently trained models, and ensembles of denoisers.

heading: H.3.1 AGGREGATION APPROACHES
text: We experiment with various instantiations of the general aggregation approach described in §4. In particular, we consider soft-voting, hard-voting, soft-voting after softmax, and weighted soft-voting.
Soft-voting We simply average the logits to obtain
f̄(x) = 1
k k∑ l=1 f l(x).
Hard-voting We process the outputs of single models with the post-processing function γHV (yl) = 1j=arg maxi(yli) that provides a one-hot encoding of the arg max of f l(x) before averaging to obtain
f̄(x) = 1
k k∑ l=1 1j=arg maxi(f l(x)i).
Soft-voting after softmax We process the outputs of the single models by applying the softmax function as post-processing function γsoftmax to obtain
f̄(x) = 1
k k∑ l=1 exp(f l(x))∑ i exp(f l(x)i) .
Weighted soft-voting We consider soft-voting with a classifier-wise weights wl learned on a separate holdout set and obtain
f̄(x) = 1
k k∑ l=1 wlf l(x).
We compare these approaches in Table 10 and observe that the two soft-voting schemes perform very similarly and outperform the other approaches. We decide to use soft-voting without soft-max for its conceptual simplicity for all other experiments.

heading: H.3.2 ENSEMBLE SIZE
text: In Table 11, we present extended results on the effect of ensemble size for different training methods and σ . Generally, we observe that ensembles of sizes 3 and 5 already significantly improve over the performance of a single model, while further increasing the ensemble size mostly leads to marginal gains. This aligns well with our theory from §5.
H.3.3 VARIABILITY OF RESULTS
In Table 12, we report mean and standard deviation of ACR and certified accuracy for 50 GAUSSIAN and CONSISTENCY trained ResNet20 with different random seeds, either combined to 10 ensembles with k = 5 or evaluated individually. Generally, we observe a notably smaller standard deviation of the ensembles compared to individual models. We visualize this in Fig. 14, where we show the ±3σ region.

heading: H.3.4 EQUAL NUMBER OF INFERENCES
text: In Table 13, we compare the certified accuracy of an ensemble of k = 10 ResNet110 obtained with n = 100′000 samples
CONSISTENCY 1 0.5202 ± 0.0049 63.57 ± 0.86 54.22 ± 1.03 44.30 ± 0.895 0.5445 ± 0.0017 64.90 ± 0.66 56.86 ± 0.49 47.94 ± 0.68
on CIFAR10 with that of an individual model obtained with n = 1′000′000, i.e. at the same computational cost.
For both GAUSSIAN and CONSISTENCY training, we observe that at all radii which can mathematically be certified using n = 100′000 samples the ensemble it significantly outperforms the individual model at equal computational cost. Note that the individual models with n = 1′000′000 samples have the same or just a marginally larger certified accuracy at all radii ≤ 1.75 than the same individual model with n = 100′000 samples. Only at radius 2.00 where certification is mathematically impossible with n = 100′000 samples does using additional samples increase the certified accuracy. Note however, that, in contrast to our ensembles, this does not influence the actual accuracy of the underlying model in any way, but simply allows to derive tighter confidence bounds. Table 14 shows similar results for ImageNet and ensembles of k = 3 ResNet110, implying that these observations hold for various datasets.

heading: H.3.5 ENSEMBLES OF DIFFERENTLY TRAINED MODELS
text: In Table 15 we show how ensembles built from differently trained base models (SMOOTHADV, MACER, and CONSISTENCY) perform. We observe that the ensemble of all three models is strictly better than any individual model, suggesting that ensembles also work for models trained with different training methods. However, the improvements are significantly smaller than the improvements of ensembling similar classifiers in the same setting (see also Table 8 for comparison). This is again due to a very heterogeneous model performance at some radii, where the strongest models are weighed down by weaker ones.
H.3.6 ENSEMBLES OF DENOISERS
It was shown (Salman et al., 2020) that standard neural networks can be robustified with respect to Gaussian noise by training a suitable denoiser and evaluating the original model on the denoised sample. Here, we aim to overcome the drawback of any ensemble, the need to train multiple diverse models, by instead training different denoisers and building an ensemble by combining these denoisers with a single underlying model. We show the general applicability of our ensembling approach by
applying it to the setting of denoised smoothing (Salman et al., 2020). We consider 4 pre-trained denoisers and combine them with a single standard ResNet110 (see App. G.3 for details). Comparing the strongest resulting model with an ensemble of all 4 at σ = 0.25, we observe a 17% ACR improvement. Training strong denoisers with current methods is challenging. However, based on these results we are confident that advances in denoiser training combined with ensembles can present an efficient way to obtain strong provable models.

heading: H.3.7 ADDITIONAL NORMS
text: Below, we demonstrate that our approach generalizes beyond the `2-norm setting and generalizes well to improve the robust accuracy against `0-, `1- and `∞-norm bounded attacks.
`∞-norm robustness is typically derived directly from `2-norm bounds (Yang et al., 2020; Salman et al., 2019), and hence our results for the `2-norm are directly applicable.
In Table 17 and Table 18, we show the effectiveness of our approach for `0-norm bounded perturbations using the training and certification method from Lee et al. (2019). For both, we use the default hyper-parameter α (which is not the same as the α in standard Randomized Smoothing from Cohen et al. (2019)), i.e. α = 0.8 for MNIST and α = 0.2 for ImageNet. We observe that even an ensemble of just 3 models outperforms the individual models in every setting.
In Table 19, we similarly show the effectiveness of our approach for `1-norm bounded perturbations using training with uniform noise and the uniform noise based certification method from Yang et al. (2020) with σ = λ = 1.0. We observe again that even an ensemble of just 3 models outperforms the single model in every setting.

heading: H.4 ADDITIONAL COMPUTATIONAL OVERHEAD REDUCTION EXPERIMENTS
text: This section contains additional experiments for the computational overhead methods introduced in §6. Whenever we use CERTIFYADP, we set β = 0.001.
Adaptive sampling and K-consensus Table 20 is an extensive version of Table 3. In Table 21, we show the effect of adaptive sampling andK-Consensus aggregation for GAUSSIAN trained ResNet110 on CIFAR10, analogously to Table 20. Certifying an ensemble with CERTIFYADP and K-Consensus aggregation achieves approximately the same certified accuracy as CERTIFY while the certification time is up to 67 times faster. In particular, for both training methods and for each radius, the certifying the ensemble of k = 10 networks with CERTIFYADP and K-Consensus aggregation is faster than
certifying an individual model with CERTIFY while the certified accuracy is significantly improved in most cases. In addition, we confirm our previous observation that the sample reduction is most prominent for small radii (for a given σ ) while KCR increases with radius.

heading: H.4.1 ADDITIONAL ADAPTIVE SAMPLING EXPERIMENTS
text: This section contains other experiments with CERTIFYADP, highlighting various aspects of adaptive sampling for ensembles and also individual models.
Adaptive sampling for ensembles on CIFAR10 For comparison with Table 3 (respectively Table 20 and Table 21), where we show the combined effect of K-Consensus aggregation and adaptive sampling, we show the results for just applying adaptive sampling to the same setting in Table 22, using {nj} = {100, 1′000, 10′000, 120′000}, σ = 0.25, and ensembles of 10 CONSISTENCY trained ResNet110. We observe that the speed-ups for small radii are comparable to applying both K-Consensus aggregation and adaptive sampling, as only a few samples have to be evaluated. The additional benefit due to K-Consensus aggregation grows as the later certification stages are entered more often for larger radii, highlighting the complementary nature of the two methods. Note that for radius r = 0.75, 100 samples, and even 1′000 are never sufficient for certification but can still yield early abstentions.
Adaptive sampling for ImageNet ensembles In Table 23, we demonstrate the applicability of our adaptive sampling approach to ensembles of ResNet50 on ImageNet. In particular, we achieve sampling reduction factors of up to 42 without incurring a significant accuracy penalty. For a fixed σ , we empirically observe larger speed-ups at smaller radii. This effect, however, depends on the underlying distribution of success probabilities and would be expected to invert for a model that has true success probabilities close to pA ≈ 0.5 for most samples.
Effect of sampling stage sizes Additionally, in Table 24, we investigate how the choice of the number of phases and the number of samples in each phase influence the sampling gain for various radii and models. We observe that as few as 100 samples can be sufficient to certify or abstain from up to 80% of samples, making schedules with many phases attractive. However, when the true success probability is close to the one required for certification, many samples are required to obtain sufficiently tight confidence bounds on the success probability to decide either way. This leads to many samples being discarded for schedules with many phases. Additionally, the higher per phase confidence required due to Bonferroni correction requires more samples to be drawn to be able to obtain the same maximum lower confidence bound to the success probability. This latter effect becomes especially pronounced for larger radii, where higher success probabilities are required.
Ablation of α and β In Fig. 15, we illustrate the effect of different α and β on the first phase of certification with CERTIFYADP where we consider just 100 perturbations. Increasing α reduces the confidence 1− α required for certification and hence makes early certification more likely for the same underlying true success probabilities (see Fig. 15a). Increasing β reduces the confidence 1 − β that certification will not be possible required for early termination and hence makes early termination more likely for the same underlying true success probabilities (see Fig. 15b). Overall, we observe that evaluating even only 100 perturbations allows us to, with high probability, abstain from or certify samples with a wide range of true success probabilities. We only have to continue with the costly certification process for samples in a narrow band of true success probabilities (see Fig. 15c), greatly reducing the expected certification cost.

heading: H.4.2 ADDITIONAL K-CONSENSUS EXPERIMENTS
text: In Table 25, we present a more detailed version of Table 5. We emphasize that forK = 5 andK = 10 for ResNet110 and ResNet20, respectively, we achieve the same ACR as with the full ensemble while reducing the certification time by a factor of 1.59 and 2.01, respectively.

Review summaries:
The paper aims to boost randomized smoothing (RS). Specifically, the paper demonstrates that an ensemble of diverse base models can enhance RS both theoretically and empirically. The key insight is that reducing variance of ensembles over the introduced perturbations can lead to more consistent classifications for inputs. The paper also introduce two simple  yet effective techniques to speed up the certification. Extensive experiments are conducted to thoroughly evaluate the proposed boosted RS. 

This paper integrates model ensembles with randomized smoothing to improve the certified accuracy. The methodology is motivated theoretically by showing the effect of model ensemble on reducing the variance of smooth classifiers. Moreover, it proposes an adaptive sampling algorithm to reduce the computation required for certifying with randomized smoothing. Extensive experiments were conducted on CIFAR10 and ImageNet datasets.

In this paper, the authors propose using the aggregation of an ensemble of similar models as the base classifier in the randomized smoothing (RS). They show that the use of ensembles helps reduce the variance of the base classifier under noisy inputs, thus, improving the performance of RS. Both theoretical arguments and numerical experiments are included to support their idea. Further, the authors provide practical algorithms that significantly reduce the computational costs of their method.