Year: 2021

Abstract:
Randomized smoothing is currently considered the state-of-the-art method to obtain certifiably robust classifiers. Despite its remarkable performance, the method is associated with various serious problems such as “certified accuracy waterfalls”, certification vs. accuracy trade-off, or even fairness issues. Input-dependent smoothing approaches have been proposed with intention of overcoming these flaws. However, we demonstrate that these methods lack formal guarantees and so the resulting certificates are not justified. We show that the input-dependent smoothing, in general, suffers from the curse of dimensionality, forcing the variance function to have low semi-elasticity. On the other hand, we provide a theoretical and practical framework that enables the usage of input-dependent smoothing even in the presence of the curse of dimensionality, under strict restrictions. We present one concrete design of the smoothing variance and test it on CIFAR10 and MNIST. Our design solves some of the problems of classical smoothing and is formally underlined, yet further improvement of the design is still necessary.

text: Randomized smoothing is currently considered the state-of-the-art method to obtain certifiably robust classifiers. Despite its remarkable performance, the method is associated with various serious problems such as “certified accuracy waterfalls”, certification vs. accuracy trade-off, or even fairness issues. Input-dependent smoothing approaches have been proposed with intention of overcoming these flaws. However, we demonstrate that these methods lack formal guarantees and so the resulting certificates are not justified. We show that the input-dependent smoothing, in general, suffers from the curse of dimensionality, forcing the variance function to have low semi-elasticity. On the other hand, we provide a theoretical and practical framework that enables the usage of input-dependent smoothing even in the presence of the curse of dimensionality, under strict restrictions. We present one concrete design of the smoothing variance and test it on CIFAR10 and MNIST. Our design solves some of the problems of classical smoothing and is formally underlined, yet further improvement of the design is still necessary.

heading: 1 INTRODUCTION
text: Deep neural networks are one of the dominating recently used machine learning methods. They achieve state-of-the-art performance in a variety of applications like computer vision, natural language processing, and many others. The key property that makes neural networks so powerful is their expressivity (Gühring et al., 2020). However, as a prize, they possess a weakness - a vulnerability against adversarial attacks (Szegedy et al., 2013; Biggio et al., 2013). The adversarial attack on a sample x is a point x′ such that the distance d(x′, x) is small, yet the model f ’s predictions on x and x′ differ. Such examples are often easy to construct, for example by optimizing for a change in prediction f(x) (Biggio et al., 2013). Even worse, these attacks are present even if the model’s prediction on x is unequivocal.
This property is highly undesirable because in several sensitive applications, misclassifying a sample just because it does not follow the natural distribution might lead to serious and harmful consequences. A well-known example would be a sticker placed on a traffic sign, which could possibly confuse the self-driving car and cause an accident (Eykholt et al., 2018). To prevent this behaviour, the robustness of classifiers against adversarial examples has begun to be a strongly discussed topic. Though many methods claim to provide robust classifiers, just some of them are certifiably robust, i.e. the robustness is mathematically guaranteed. The certifiability turns out to be essential since more sophisticated attacks can break empirical defenses (Carlini & Wagner, 2017).
Currently, the dominating method to achieve the certifiable robustness is randomized smoothing (RS). This clever idea to get rid of adversarial examples using randomization of input was introduced by Lecuyer et al. (2019) and Li et al. (2019) and fully formalized and improved in Cohen et al. (2019). Let f be a classifier classifying inputs x ∈ RN as one of the classes C ∈ C. Assume now a random deviation ∼ N (0, σ2I). The smoothed classifier g, made of f , is defined as: g(x) = arg maxC P(f(x+ ) = C), for C ∈ C. In other words, the smoothed classifier classifies a class that has the highest probability under the sampling of f(x + ). Consequently, an adversarial attack x′ on f is less dangerous for g, because g does not look directly at x′, but rather at its whole neighborhood, in a weighted manner. This way, we can get rid of local artifacts that f possesses – thus the name “smoothing”. It turns out, that g enjoys strong robustness properties against attacks
bounded by a specifically computed l2-norm threshold, especially if f is trained under a Gaussian noise augmentation (Cohen et al., 2019).
Unfortunately, since the introduction of the RS, several serious problems were reported to be connected to the technique. Cohen et al. (2019) mention two of them. First is the usage of lower confidence bounds to estimate the leading class’s probability. With a high probability, this leads to smaller reported certified radiuses in comparison with the true ones. Moreover, it yields a theoretical threshold, which upper-bounds the maximal possible certified radius and causes the “certified accuracy waterfalls”, significantly decreasing the certified accuracy. This problem is particularly pronounced for small levels of the used smoothing variance σ2, which motivates to use larger variance. Second, RS possesses a robustness vs. accuracy trade-off problem. The bigger σ we use as the smoothing variance, the smaller clean accuracy will the smoothed classifier have. This motivates to use rather smaller levels of σ. Third, as pointed out by Mohapatra et al. (2020a), RS smoothens the decision boundary of f in such a way, that bounded or convex regions begin to shrink as σ increases, while the unbounded and anti-convex regions expand. This, as the authors empirically demonstrate, creates a disbalance in class-wise accuracies of g and causes serious fairness issues. Therefore, again, the smaller values of σ are more preferable. See Appendix A for a detailed discussion.
Clearly, the usage of a global, constant σ is suboptimal. For the samples close to the decision boundary, we want to use small σ, so that f and g have similar decision boundaries and the expressivity of f is not lost (where not necessary). On the other hand, far from the decision boundary of f , where the probability of the dominant class is close to 1, we need bigger σ to avoid the severe under-certification (see Appendix A). Together, using a non-constant σ(x) rather than constant σ, a suitable smoothing variance could be used to achieve optimal robustness. Yet there are some works introducing this concept (see Appendix C), most of them lack mathematical reasoning about the correctness of their method, which, as we show, turns out to be critical.
To support this argumentation, we present a toy example. We train a network on a 2D dataset of a circular shape with the classes being two complementary sectors, where one is of a very small angle. In Figure 1 we show the difference between constant and input-dependent σ. Using non-constant σ(x) defined in Equation 1, we obtain an improvement both in terms of the certified radiuses as well as clean accuracy. For more details see Appendix A.
The main contributions of this work are fourfold. First, we generalize the methodology of Cohen et al. (2019)’s work for the case of input-dependent RS (IDRS), obtaining useful and important insights about how to use the Neyman-Pearson lemma in this, general, case. Second and most importantly, we show that the IDRS suffers from the curse of dimensionality in the sense that the semi-elasticity coefficient r of σ(x) (that is | log(σ(x0))−log(σ(x1))| ≤ r‖x0−x1‖ ∀x0, x1 ∈ RN ) in a high-dimensional setting is restricted to be very small. This means, that even if we wanted to vary the σ(x) significantly with varying x, we can’t. The maximal reasonable speed of change of σ(x) turns out to be almost too small to handle, especially in high dimensions. Third, in contrast,
we also study the conditions on σ(x) under which it is applicable in high-dimensional regime and prepare a theoretical framework necessary to build an efficient certification algorithm. We are the first to do so for σ(x) functions, which are not locally constant (as in Wang et al. (2021)). Finally, we provide a concrete design of the σ(x) function and test it extensively and compare it to the classical RS on the CIFAR10 and MNIST datasets. We discuss to what extent the method treats the issues mentioned above.

heading: 2 INPUT-DEPENDENT RS AND THE CURSE OF DIMENSIONALITY
text: Let C be the set of classes, f : RN −→ C a classifier (reffered to as the base classifier), σ : RN −→ R a non-negative function and P(C) a set of distributions over C. Then we call Gf : RN −→ P(C) the smoothed class probability predictor, if Gf (x)C = P(f(x + ) = C), where ∼ N (0, σ(x) 2I) and gf : RN −→ C is called smoothed classifier if gf (x) = arg maxC Gf (x)C , for C ∈ C. We will omit the subscript f in gf often, since it is usually clear from the context, to which base classifier the g corresponds. Furthermore, let A := g(x) refer to the most likely class under the random variable f(N (x, σ2I)), while B denote the second most likely class. Define pA = Gf (x)A and pB = Gf (x)B as the respective probabilities. It is important to note that in practice, it is impossible to estimate pA and pB precisely. Instead, pA is estimated as a lower confidence bound (LCB) of the relative occurence of class A in f ’s predictions given certain number of Monte-Carlo samples n and a confidence level α and the estimate is denoted as pA. We use the exact Clopper-Pearson interval for estimation of the LCB. Similarly for pB . We work with l2-norms denoted as ‖x‖. First of all, we summarize the main steps in the derivation of certified radius using any method that relies on a use of Neyman-Pearson lemma.
1. For a potential adversary x′ specify the worst-case classifier f∗, such that P(f∗(N (x, σ2I)) = A) = pA, while P(f∗(N (x′, σ2I)) = B) is maximized. 2. Express the probability Gf∗(x′)B as a function depending on x ′. 3. Determine the conditions on x′ (possibly related to ‖x− x′‖) for which this probability is ≤ 0.5. From these condtions, derive the certified radius.
Cohen et al. (2019) proceeded in this way to obtain a tight certified radius R = σ2 (Φ −1(pA) − Φ−1(pB)). Unfortunately, their result is not directly applicable to the input-dependent case. The constant σ simplifies the derivation of f∗ that turns out to be a linear classifier. This is not the case for non-constant σ(x) anymore. Therefore, we need to generalize the methodology of Cohen et al. (2019). We put pB = 1−pA for simplicity (yet it is not necessary to assume this, see Appendix B.5). Let x0 be the point to certify, x1 the potential adversary point, δ = x1−x0 the noise and σ0 = σ(x0), σ1 = σ(x1) the standard deviations used in x0 and x1 respectively. Furthermore, let fi be a density and Pi a probability measure corresponding to N (xi, σ2i I), i ∈ {0, 1}. Lemma 1. Out of all possible classifiers f such that Gf (x)B ≤ pB = 1 − pA, the one, for which Gf (x+ δ)B is maximized is the one, which predicts classB in a region determined by the likelihood ratio:
B = { x ∈ RN : f1(x)
f0(x) ≥ 1 r
} ,
where r is fixed, such that P0(B) = pB . Note, that we use B to denote both the class and the region of that class.
We use this lemma to compute the decision boundary of the worst-case classifier f∗. Theorem 2. If σ0 > σ1, then B is a N -dimensional ball with the center at S> and radius R>:
S> = x+ σ20
σ20 − σ21 δ, R> =
√ σ20σ 2 1
(σ20 − σ21)2 ‖δ‖2 + 2N σ
2 0σ 2 1
σ20 − σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 σ20 − σ21 log(r).
If σ0 < σ1, then B is the complement of a N -dimensional ball with the center at S< and radius R<:
S< = x− σ20
σ21 − σ20 δ, R< = √ 2σ40 − σ20σ21 (σ21 − σ20)2 ‖δ‖2 + 2N σ 2 0σ 2 1 σ21 − σ20 log ( σ1 σ0 ) − 2σ 2 0σ 2 1 σ21 − σ20 log(r).
As we depict on Figure 2, both resulting balls are centered on the line connecting x0, x1. Moreover, the centers of the balls are always further from x0, than x1 is from x0 (even in the case σ0 < σ1). In both cases, it depends on pA (since r is fixed such that P0(B) = pB) and the ball can, but might not cover x0 and/or x1. Note that if σ0 = σ1, what can happen even in input-dependent regime, the worst-case classifier is the half-space described by Cohen et al. (2019).
To compute the probability of a ball under a probability measure with an isotropic Gaussian density is far more challenging than to compute the probability of a half-space. In fact, there is no closedform formula for such probability. However, this probability is connected to the non-central chisquare distribution (NCCHSQ). More precisely, the probability of an N -dimensional ball centered at z with radius r under N (0, I) can be expressed as a cumulative distribution fucntion (cdf) of NCCHSQ with N degrees of freedom, non-centrality parameter ‖z‖2 and argument r2. With this knowledge, we can express P0(B) and P1(B) in terms of the cdf of NCCHSQ as follows. Theorem 3.
P0(B) = χ2N
( σ20
(σ20 − σ21)2 ‖δ‖2, R2<,> σ20
) ,P1(B) = χ2N ( σ21
(σ20 − σ21)2 ‖δ‖2, R2<,> σ21
) ,
where the sign < or > is chosen according to the inequality between σ0 and σ1.
Note, that both Theorem 2 and Theorem 3 work well also for δ = 0. In this case, we encounter a ball centered at x0 = x1 and all the cdf functions become cdf functions of central chi-squared.
We expressed probabilities of the worst-case class B’s decision region using the cdf of NCCHCSQ. Now, how do we do the certification? We start with the certification just for two points, x0 and x1. We question, under which circumstances can x1 be certified from the point of view of x0. Having x0, pA and σ0 > σ1, we can obtain such R, that P0(B) = χ2N ( ‖δ‖2σ20/(σ20 − σ21)2, R2 ) = 1 − pA =
pB simply by putting it into the quantile function: R2 = χ2N,qf ( ‖δ‖2σ20/(σ20 − σ21)2, 1− pA ) . Then,
we can substitute into P1(B) = χ2N ( ‖δ‖2σ21/(σ20 − σ21)2, R2σ20/σ21 ) . This way, we obtain P1(B) and can judge, whether P1(B) < 1/2 or not. Similar computation can be done if σ0 < σ1. Denote a := ‖δ‖. We can express the P1(B) more simply for σ0 > σ1 as
ξ>(a) := P1(B) = χ2N (
σ21 (σ20 − σ21)2 a2, σ20 σ21 χ2N,qf
( σ20
(σ20 − σ21)2 a2, 1− pA )) and for σ0 < σ1 as
ξ<(a) := P1(B) = 1− χ2N (
σ21 (σ21 − σ20)2 a2, σ20 σ21 χ2N,qf
( σ20
(σ21 − σ20)2 a2, pA
)) .
With this in mind, if we have x0, x1, pA, σ0, σ1, then we can certify x1 w.r.t x0 simply by choosing the correct sign (<,>) and computing ξ<(‖x0 − x1‖) or ξ>(‖x0 − x1‖), comparing them with 0.5. The sample plots of these ξ functions can be found in Appendix B.
Now, we are ready to discuss the curse of dimensionality. The problem that arises is that having a high dimension N and σ0, σ1 differing a lot from each other, ξ functions are already big at 0, even
for considerably small pB . For fixed ratio σ0/σ1 and probability pB , with increasing dimension, the ξ(0) increases and soon becomes bigger than 0.5. This, together with monotonicity of ξ function yields that any x1 cannot be certified w.r.t. x0, if σ0, σ1 are used. The more dissimilar the σ0 and σ1 are, the smaller the dimension N needs to be for this situation to occur. If we want to certify x1 in a reasonable distance from x0, we need to use similar σ0, σ1. This restricts the variability of σ(x) function. We will formalize the curse of dimensionality in the following theorems. More on why the curse of dimensionality is present is in Appendix B.2. Theorem 4 (curse of dimensionality). Let x0, x1, pA, σ0, σ1, N be as usual. Then, the following two implications hold:
1. If σ0 > σ1 and
log ( σ21 σ20 ) + 1− σ 2 1 σ20 < 2 log(1− pA) N ,
then x1 is not certified w.r.t. x0.
2. If σ0 < σ1 and
log ( σ21 σ20 N − 1 N ) + 1− σ 2 1 σ20 N − 1 N < 2 log(1− pA) N ,
then x1 is not certified w.r.t. x0. Corollary 5 (one-sided simpler bound). Let x0, x1, pA, σ0, σ1, N be as usual and assume now σ0 > σ1. Then, if
σ1 σ0 <
√ 1− 2 √ − log(1− pA)
N ,
then x1 is not certified w.r.t x0.
Note, that both Theorem 4 and Corrolary 5 can be adjusted to the case where we have a separate estimate pB and do not put pB = 1− pA (see Appendix B.5). We must emphasize, that the bounds obtained in Theorem 4 are very tight. In other words, if the ratio σ1σ0 is just slightly bigger than the minimal possible threshold determined in Theorem 4, ξ>(0) becomes smaller than 0.5 and similarly for ξ>(0). The reason for this is, that the only two estimates used in the proof of Theorem 4 are the estimates on the median, which are very tight and constant with respect to N and the Chernoff bound, which is generally considered to be tight too and improves for larger N . The tightness is depicted on Figure 3, where we plot the minimal possible threshold σ1/σ0 given by Theorem 4 and minimal threshold for which ξ>(0) < 0.5 as a function of N .
To get a better feeling about the concrete numbers, we provide a simple table, which shows the theoretical threshold values provided by Theorem 4. If σ1/σ0 is smaller than the threshold, we are not able to certify any pair of x0, x1 using σ0, σ1.
Results from Table 1 are very restrictive. Assume we have a CIFAR10 sample with pA = 0.999. For such a probability, constant σ = 0.5 is more than sufficient to guarantee the certified radius of more than 1. However, in the non-constant regime, to certify R ≥ 1, we first need to guarantee that no sample within the distance of 1 from x0 uses σ1 < 0.953σ0. To even strengthen this statement, note that one needs to guarantee σ1 to be even much closer to σ0 in practice. Why? The results of Theorem 4 lower-bound the ξ functions at 0. However, since ξ functions are strictly increasing (as shown in Appendix F), one usually needs σ0 and σ1 to be much closer to each other to guarantee ξ being smaller than 0.5 at a 0. This not only forces the σ(x) function to have really small semielasticity but also makes it problematic to define a stochastic σ(x). For more, see Appendix B.2.
To fully understand, how the curse of dimensionality affects the usage of IDRS, we mention two more significant effects. First, with increasing dimension, the average distance between samples tends to grow as √ N . This enables bigger distance to change σ(x). On the other hand, the average level of σ(x) used (like ∼ 0.12, 0.25, . . . ) needs to be adjusted also as √ N with increasing dimension. The bigger average level of σ(x) we use, the more is the semi-elasticity of σ(x) restricted by Theorem 4 and Theorem 7. All together, these two effects combine in a final trend that for σ0 and σ1 being variances used in two random test samples, |σ0/σ1 − 1| is restricted to go to 0 as 1/ √ N. For detailed explanation, see Appendix B.4.

heading: 3 HOW TO USE INPUT-DEPENDENT SMOOTHING PROPERLY
text: As we discuss above, usage of IDRS is challenging. How can we use σ(x) and obtain valid, mathematically justified certified radiuses? Fix some design σ(x). If σ(x) is not trivial, to obtain a certified radius at x0, we need to go over all the possible adversaries x1 in the neighborhood of x0, compute σ1 and ξ<,>(a). Then, the certified radius is the infimum over ‖x0−x1‖ for all uncertified x1 points. Of course, this is a priori infeasible. Fortunately, the ξ functions possess a property that helps to simplify this procedure. For convenience, we extend the notation of ξ functions such that ξ(a, σ1) additionally denotes the dependance on the σ1 value. Theorem 6. Let x0, x1, pA, σ0 be as usual and let ‖x0 − x1‖ = R. Then, the following two statements hold:
1. Let σ1 ≤ σ0. Then, for all σ2 : σ1 ≤ σ2 ≤ σ0, if ξ>(R, σ2) > 0.5, then ξ>(R, σ1) > 0.5.
2. Let σ1 ≥ σ0. Then, for all σ2 : σ1 ≥ σ2 ≥ σ0, if ξ<(R, σ2) > 0.5, then ξ<(R, σ1) > 0.5.
Theorem 6 serves as a kind of monotonicity property. The main gain is that now, for each distance R from x0, we need to pick just two adversaries – the one with biggest σ1 (if bigger than σ0) and the one with the smallest σ1 (if smaller than σ0). If we cannot certify some point x1 at the distance R from x0, then we will for sure not be able to certify at least one of the two adversaries with the most extreme σ1 values.
This does, however, not suffice for most of the reasonable σ(x) designs, since it might be still too hard to determine the two most extreme σ1’s at some distance from x0. Therefore, we need to assume that our σ(x) has a bounded semi-elasticity coefficient r. Then we have a guarantee that σ(x0) exp(−ra) ≤ σ(x1) ≤ σ(x0) exp(ra). Thus, we can assume the worst-case extreme σ1’s for every distance a from x0. Using this, we guarantee the following certified radius. Theorem 7. Let σ(x) be r-semi-elastic function and x0, pA, N, σ0 as usual. Then, the certified radius at x0 guaranteed by our method is
CR(x0) = max {0, sup {R ≥ 0; ξ>(R, σ0 exp(−rR)) < 0.5 and ξ<(R, σ0 exp(rR)) < 0.5}} .
Note that Theorem 7 can be adjusted to the case where we have a separate estimate pB and do not put pB = 1 − pA (see Appendix B.5). Since the bigger the semi-elasticity constant of σ(x) is, the worse certifications we obtain, it is important to estimate the constant tightly. Even with a good estimate of r, we still get smaller certified radiuses in comparison with using the σ(x) exactly, but that is a prize that is inevitable for the feasibility of the method.
The practical algorithm is then very easy - we just need to pick sufficiently dense space of possible radiuses and determine the smallest, for which either ξ>(R, σ0 exp(−rR)) or ξ<(R, σ0 exp(rR)) becomes bigger than a half. The only non-trivial part is, how to evaluate the ξ functions. For small values of R, the exp(−rR) is very close to 1 and from the definition of ξ functions it is obvious that this results in extremely big inputs to the cdf and quantile function of NCCHSQ. To avoid numerical problems, we employ a simple hack where we assume thresholds for σ1 such that for R small enough, these thresholds are used instead of σ0 exp(±rR)). Unfortunately, the numerical stability still disables the usage of this method on really high-dimensional datasets like ImageNet. For more details on implementation, see Appendix D.
4 THE DESIGN OF σ(x) AND EXPERIMENTS
The only missing ingredient to finally being able to use IDRS is the σ(x) function. As we have seen, this function has to be r-semi-elastic for rather small r and ideally deterministic. Yet it should at least roughly fulfill the requirements imposed by the motivation – it should possess big values for points far from the decision boundary of f and rather small for points close to it. Adhering to these restrictions, we use the following function:
σ(x) = σb exp r 1 k  ∑ xi∈Nk(x) ‖x− xi‖ −m  , (1)
for σb being a base standard deviation, r being the required semi-elasticity, {xi}di=1 the training set, Nk(x) the k nearest neighbors of x and m the normalization constant. Intuitively, if a sample is far from all other samples, it will be far from the decision boundary, unless the network overfits to this sample. On the other hand, the dense clusters of samples are more likely to be positioned near the decision boundary, since such clusters have a high leverage on the network’s weights, forcing the decision boundary to adapt well to the geometry of the cluster. To use such a function, however, we first prove that it is indeed r-semi-elastic. Theorem 8. The σ(x) defined in equation 1 is r-semi-elastic.
We test our IDRS and our σ(x) functions extensively. For both CIFAR10 (Krizhevsky, 2009) and MNIST (LeCun et al., 1999) datasets, we analyze series of different experimental setups, including experiments with an input-dependent train-time Gaussian data augmentation. We present a direct comparison of our method with the constant σ method using evaluation strategy from Cohen et al. (2019) (all other experiments, including ablation studies, and the discussion on the hyperparameter selection are presented in Appendix E). Here, we compare Cohen et al. (2019)’s evaluations for σ = 0.12, 0.25, 0.50 with our evaluations, setting σb = σ, r = 0.01, k = 20, m = 5, 1.5 (for CIFAR10 and MNIST, respectively), applied on models trained with Gaussian data augmentation, using constant standard deviation roughly equal to the average test-time σ(x) or test-time σ. For CIFAR10, these levels of train-time standard deviation are σtr = 0.126, 0.263, 0.53 and for MNIST σtr = 0.124, 0.258, 0.517. In this way, the levels of σ(x) we use in the direct comparison are spread from the values roughly equal to Cohen et al. (2019)’s constant σ to higher values. The results are depicted in Figure 4.
From Figure 4 we see that we outperform the constant σ for small levels of σ, such as 0.12 or 0.25. On higher levels of σ, we are, in general, worse (see explanation in Appendix B.3). The most visible improvement is in mitigation of the truncation of certified accuracy (certified accuracy waterfall). To comment on the other two issues, we provide Tables 2 and 3 with the clean accuracies and classwise accuracy standard deviations. These results are averages of 8 independent runs and in Table 2, the displayed error values are equal to empirical standard deviations.
From Tables 2 and 3, we draw two conclusions - first, it is not easy to judge about the robustness vs. accuracy trade-off, because the differences in clean accuracies are not statistically significant in
any of the experiments (not even for CIFAR10 and σ = 0.25, where the difference is at least pronounced). However, the general trend in Table 2 indicates that the clean accuracies tend to slightly decrease with increasing rate. Nevertheless, the differences are not large enough to compensate negatively the fact that we outperform constant σ in terms of the certified accuracies. Second, the standard deviations of the class-wise accuracies, which serve as a good measure of the impact of the shrinking phenomenon and subsequent fairness, don’t significantly change after applying the non-constant RS.

heading: 5 RELATED WORK
text: Since the vulnerability of deep neural networks against adversarial attacks has been noticed by Szegedy et al. (2013); Biggio et al. (2013) a lot of effort has been put into making neural nets more robust. There are two types of solutions – empirical and certified defenses. While empirical defenses suggest heuristics to make models robust, certified approaches additionally provide a way to compute a mathematically valid robust radius.
One of the most effective empirical defenses, adversarial training (Goodfellow et al., 2014; Kurakin et al., 2016; Madry et al., 2017), is based on a very intuitive idea to use adversarial examples for training. Unfortunately, together with adversarial training, other promising empirical defenses were subsequently broken by more sophisticated adversarial methods (for instance Carlini & Wagner (2017); Athalye & Carlini (2018); Athalye et al. (2018), among many others).
Among many certified defenses (Tsuzuku et al., 2018; Anil et al., 2019; Hein & Andriushchenko, 2017; Wong & Kolter, 2018; Raghunathan et al., 2018; Mirman et al., 2018; Weng et al., 2018), one of the most successful yet is RS. While Lecuyer et al. (2019) introduced the method within the context of differential privacy, Li et al. (2019) proceeded via the knowledge of Rényi divergences. Possibly the most prominent work on RS is that of Cohen et al. (2019), where authors fully established RS and proved tight certification guarantees.
Later, a lot of authors further worked with RS. The work of Yang et al. (2020) generalizes the certification provided by Cohen et al. (2019), to certifications with respect to the general lp norms and provide the optimal smoothing distributions for each of the norms. Other works point out different problems or weaknesses of RS like the curse of dimensionality (Kumar et al., 2020; Hayes, 2020; Wu et al., 2021), robustness vs. accuracy trade-off (Gao et al., 2020) or a shrinking phenomenon(Mohapatra et al., 2020a), which yields serious fairness issues (Mohapatra et al., 2020a).
The work of Mohapatra et al. (2020b) improves RS further by introducing the first-order information about g. In this work, authors not only estimate g(x), but also ∇g(x), making more restrictions on the possible base models f that might have created g. Zhai et al. (2020) and Salman et al. (2019) improve the training procedure of f to yield better robustness guarantees of g. Salman et al. (2019) directly use adversarial training of the base classifier f . Finally, Zhai et al. (2020) introduce socalled soft smoothing, which enables to compute gradients directly for g and construct a training method, which optimizes directly for the robustness of g via the gradient descent.
To address several issues connected to randomized smoothing, there have been already four works that introduce the usage of IDRS. Wang et al. (2021) divide RN into several regions Ri, i ∈ {1, . . . ,K} and optimize for σi, i ∈ {1, . . . ,K} locally, such that σi is a most suitable choice for the region Ri. Yet this work partially solves some problems of randomized smoothing, it also possesses some practical and philosophical issues (see Appendix C). Alfarra et al. (2020); Eiras et al. (2021); Chen et al. (2021), suggest to optimize for locally optimal σi, for each sample xi from the test set. A similar strategy is proposed by these works in the training phase, with the intention of obtaining the base model f that is most suitable for the construction of the smoothed classifier g. They demonstrate, that by using this input-dependent approach, one can overcome some of the main problems of randomized smoothing. However, as we demonstrate in Appendix C, their methodology is not valid and therefore their results are not trustworthy.

heading: 6 CONCLUSIONS
text: We show in this work that input-dependent randomized smoothing suffers from the curse of dimensionality. In the high-dimensional regime, the usage of input-dependent σ(x) is being put under strict constraints. The σ(x) function is forced to have very small semi-elasticity. This is in conflict with some recent works, which have used the input-dependent randomized smoothing without mathematical justification and therefore claim invalid results. It seems that input-dependent randomized smoothing has limited potential of improvement over the classical, constant-σ RS. Moreover, due to numerical instability, the computation of certified radiuses on high-dimensional datasets like ImageNet remains to be an open challenge.
On the other hand, we prepare a ready-to-use mathematically underlined framework for the usage of the input-dependent RS and show, that it works well for small to medium-sized problems. We also show via extensive experiments, that our concrete design of the σ(x) function reasonably treats the truncation issue connected to constant-σ RS and is partially capable of mitigating the robustness vs. accuracy one. The most intriguing and promising direction for the future work lies in the development of new σ(x) functions, which are capable of even better treatment of the mentioned issues.

heading: 7 REPRODUCIBILITY STATEMENT
text: Both our theoretical and practical results and experiments are reproducible. In the theoretical part, we provide all the relevant proofs, insights and all the important reasoning. We use public datasets for our experiments and our code contains just public, well-known libraries. The code will be publicly available after the review process. We do not use seeds for stochastic algorithms, but the level of variance is not large enough to obtain qualitatively different results. Upon request, we are willing to provide models trained by us to obtain exactly the same results.

heading: A.1 TOY EXAMPLE
text: To better demonstrate our ideas, we prepared a two-dimensional simple toy dataset. This dataset can be seen in Figure 5. The dataset is generated in polar coordinates, having uniform angle and the distance distributed as a square root of suitable chi-square distribution. The classes are positioned in a circle sectors, one in a sector with a very sharp angle. The number of training samples is 500 for each class, number of test samples is 100 for each class (except demonstrative figures, where we increased it to 300). The model that was trained on this dataset was a simple fully connected three-layer neural network with ReLU activations and a maximal width of 20.

heading: A.2 UNDERCERTIFICATION CAUSED BY THE USE OF LOWER CONFIDENCE BOUNDS
text: As we mention in Section 1, one can not usually obtain exact values of pA and pB . However, it is obvious, that for vast majority of evaluated samples, pA < pA and pB > pB . Given the nature of our certified radius, it follows that R < R, where R denotes the certified radius coming from the certification procedure with pA and pB , while R here stands for the certified radius corresponding to true values pA, pB .
It is, therefore, clear, that we face a certain level of under-certification. But how serious undercertification it is? Assume the case with a linear base classifier. Imagine, that we move the point x further and further away from the decision boundary. Therefore, pA −→ 1. At some point, the
probability will be so large, that with high probability, all n samplings in our evaluation of pA will be classified as A, obtaining p̂A = 1 - the empirical probability. The lower confidence bound pA is therefore bounded by having p̂A = 1. Thus, from some point, the certification will yield the same pA regardless of the true value of pA. So in practice, we have an upper bound on the certified radius R in the case of the linear boundary. In Figure 7 (left), we see the truncation effect. Using σ = 1, from a distance of roughly 4, we can no longer achieve a better certified radius, despite its theoretical value equals the distance. Similarly, if we fix a distance of x from decision boundary and vary σ, for very small values of σ, the value of Φ−1(pA) will no longer increase, but the values of σ will pull R towards zero. This behaviour is depicted in Figure 7 (right).
We can also look at it differently - what is the ratio between Φ−1(pA) and Φ−1(pA) for different values of pA? Since R = σΦ−1(pA) and R = σΦ−1(pA), the ratio represents the “undercertification rate”. In Figure 6 we plot Φ−1(pA)
Φ−1(pA) as a function of pA for two different ranges of values.
The situation is worst for very small and very big values of pA. In the case of very big values, this can be explained due to extreme nature of Φ−1. For small values of pA, it can be explained as a consequence of a fact, that even small difference between pA and pA will yield big ratio between Φ−1(pA) and Φ−1(pA) due to the fact, that these values are close to 0.
If we look at the left plot on Figure 8 we see, that the certified accuracy plots also possess the truncations. Above some radius, no sample is certified anymore. The problem is obviously more serious for small values of σ. On the right plot of Figure 8, we see, that samples far from the decision boundary are obviously under-certified. We can also see, that certified radiuses remain constant, even though in reality they would increase with increasing distance from the decision boundary.
All the observations so far motivate us to use rather large values of σ in order to avoid the truncation problem. However, as we will see in the next sections, using a large σ carries a different, yet equally serious burden.

heading: A.3 ROBUSTNESS VS. ACCURACY TRADE-OFF
text: As we demonstrate in the previous subsection, it is be useful to use large values of σ to prevent the under-certification. But does it come without a prize? If we have a closer look at Figure 8 (right), we might notice, that the accuracy on the threshold 0, i.e. “clean accuracy”, decreases as σ increases. This effect has been noticed in the literature (Cohen et al., 2019; Gao et al., 2020; Mohapatra et al., 2020a) and is called robustness vs. accuracy tradeoff.
There are several reasons, why this problem occurs. Generally, changing σ changes the decision boundary of g and we might assume, that due to the high complexity of the boundary of f , the decision boundary of g becomes smoother. If σ is too large, however, the decision boundary will be so smooth, that it might lose some amount of the base classifier’s expressivity. Another reason for the accuracy drop is also the increase in the number of samples, for which the evaluation is abstained. This is because using big values of σ makes more classes “within the reach of our distribution”, making the pA and pA small. If pA < 0.5 and we do not estimate pB but set pB = 1 − pA, then
we are not able to classify the sample as class A, yet we cannot classify it as a different class either, which forces us to abstain. To demonstrate these results, we computed not only the clean accuracies of Cohen et al. (2019) evaluations but also the abstention rates. Results are depicted in the Table 4.
From the table, it is obvious, that the abstention rate is possibly even bigger cause of accuracy drop than the “clean misclassification”. This problem can be partially solved if one estimated pB together
with pA too. In this way, using big σ yields generally small estimated class probabilities, but since pA ≥ pB , the problematic pB ≥ pA occur just very rarely. Another option is to increase the number of Monte-Carlo samplings for the classification decision, what is almost for free.
Yet another reason for the decrease in the accuracy is the so-called shrinking phenomenon, which we will discuss in the next subsection.
In contrast with the truncation effect, the robustness vs. accuracy trade-off motivates the usage of smaller values of σ in order to prevent the accuracy loss, which is definitely a very serious issue.

heading: A.4 SHRINKING PHENOMENON
text: How exactly does the decision boundary of g change, as we change the σ? For instance, if f is a linear classifier, then the boundary does not change at all. To answer this question, we employ the following experiment: For our toy base classifier f on our toy dataset, we increase σ and plot the heatmap of f, g, together with its decision boundary. This experiment is depicted on Figure 9. As we see from the plots, increasing σ causes several effects. First of all, the heatmap becomes more and more blurred, what proves, that stronger smoothing implies stronger smoothness.
Second, crucial, effect is that the bigger the σ, the smaller the decision boundary of a submissive class is. The shrinkage becomes pronounced from σ = 1. Already for σ = 4, there is hardly any decision boundary anymore. Generally, as σ −→ ∞, g will predict the class with the biggest volume in the input space (in the case of bounded input space, like in image domain, this is very well defined). For extreme values of sigma, the pA will practically just be the ratio between the volume of A and the actual volume of the input space (for bounded input spaces).
Following from these results, but also from basic intuition, it seems, that an undesired effect becomes present as σ increases - the bounded/convex regions become to shrink, like in Figure 9, while the unbounded/big/anti-convex regions expand. This is called shrinking phenomenon. Mohapatra et al. (2020a) investigate this effect rather closely. They define the shrinkage and vanishing of regions formally and prove rigorously, that if σ −→ ∞, bounded regions, or semi-bounded regions (see Mohapatra et al. (2020a)) will eventually vanish. We formulate the main result in this direction.
Theorem 9. Let us haveK the number of classes and the dimensionN . Assume, that we have some bounded decision region D for a specific class roughly centered around 0. Further assume, that this is the only region where the class is classified. Let R be a smallest radius such that D ⊂ BR(0). Then, this decision region will vanish at most for σ ≥ R
√ K√ N .
Proof. The idea of the proof is not very hard. First, the authors prove, that the smoothed region will be a subset of the smoothed BR(0). Then, they upper-bound the vanishing threshold of such a ball in two steps. First, they show, that if 0 is not classified as the class, then no other point will be (this is quite an intuitive geometrical statement. The BR(0) has the biggest probability under N (x, σ2I) if x ≡ 0). Second, they upper-bound the threshold for σ, under which BR(0) will have probability below 1K (since they use slightly different setting as Cohen et al. (2019)) for the N (x, σ
2I). Using some insights about incomplete gamma function, which is known to be also the cdf of central chisquare distribution, and some other integration tricks, they obtain the resulting bound.
Besides Theorem 9, authors also claim many other statements abound shrinking, including shrinking of semi-bounded regions. Moreover, they also conduct experiments on CIFAR10 and ImageNet to support their theoretical findings. They also point out serious fairness issue that comes out as a consequence of the shrinkage phenomenon. For increasing levels of σ, they measure the class-wise clean accuracy of the smoothed classifier. If f is trained with gaussian data augmentation (what is known to be a good practice in randomized smoothing), using σ = 0.12, the worst class cat has test accuracy of 67%, while the best class automobile attains the accuracy of 92%. The figures, however, change drastically, if we use σ = 1 instead. In this case, the worst predicted class cat has accuracy of poor 22%, while ship has reasonable accuracy 68%. As authors claim, this is a consequence of the fact that samples of cat are situated more in bounded, convex regions, that suffer from shrinking, while samples of ship are mostly placed in expanded regions of anti-convex shape that will expand as the σ grows. In addition, the authors also show, that the gaussian data augmentation or adversarial training will reduce the shrinking phenomenon just partially and for moderate and high values of σ, this effect will be present anyway.
We must emphasize, that this is a serious fairness issue, that has to be treated before randomized smoothing can be fully used in practice. For instance, if we trained a neural network to classify humans into several categories, fairness of classification is inevitable and the neural network cannot be used until this issue is solved.
Similarly as the robustness vs. accuracy tradeoff, this issue also motivates to use rather smaller values of σ. We see, that it is not possible to address all three problems consistently because they disagree on whether to use smaller, or bigger values of σ .

heading: A.5 EXPERIMENTS ON HIGH-DIMENSIONAL TOY DATASET
text: In this subsection, we present the results of our motivational experiment on a synthetic dataset. Before reading this section, please read our main text, because we will use the necessary notation of the paper.
The dataset we evaluated our method on is a generalization of the dataset visualized on Figure 1. The data points from one class lie in a cone of small angle and the points are generated such that the density is higher near the vertex of the cone (which is put in origin). Points from other class are generated from a spherically symmetrical distribution (where points sampled into the cone are excluded) with density again highest in the center (note, that the density peak is more pronounced than in the case of normal distribution, where the density around the center resembles uniform distribution). This dataset is chosen so that the σ(x) function designed in Equation 1 well corresponds to the geometry of the decision boundary. Moreover it is chosen so that the conic decision region
will shrink rather fast with increasing σ. The motivation of this example is to show that if the σ(x) function is well-designed, our IDRS can outperform the constant RS considerably.
The setup of our experiment is as follows: We evaluate dimensions N = 2, 6, 18, 60, 180, 400. The σ used for constant smoothing is σ = 0.5, 0.5, 1.0, 1.0, 2.0, 2.0 respectively. The σb used is 0.4, 0.5 for N = 2, 0.4 for N = 6, 0.8 for N = 18, 0.8, 1.0 for N = 60, 1.9 for N = 180 and 1.95 for N = 400. The rates are r = 0.2, 0.1, 0.05, 0.03, 0.01, 0.005 respectively. The training was executed without data augmentation (because samples from different classes are very close to each other). Moreover, we have set maximal σ(x) threshold for numerical purposes, because some samples were outliers and were way too far from other samples (and if the σ(x) is way too big, the method encounteres numerical problems). In this case we set σ(x) ≤ 5σb, but we are aware that also much bigger thresholds would have been possible. We present our comparisons in Figure 10 and Table 5.
From both the Figure 10 an Table 5 it is clear that the IDRS can outperform the constant σ RS considerably, if we use really suitable σ(x) function. We manage to improve significantly the certified radiuses without losing a single correct classification. On the other hand, in cases where σb < σ, we outperform constant σ both in clean accuracy and in certified radiuses. This example is synthetic and designed in our favour. The main message is not how perfect our design of σ(x) is, but the fact, that if σ(x) is designed well, the IDRS can bring real advantages, even in moderate dimensions.

heading: B MORE ON THEORY

heading: B.1 GENERALIZATION OF RESULTS FROM LI ET AL. (2019)
text: In our main text, we mostly focus on the generalization of the methods from Cohen et al. (2019). This is because these methods yield tight radiuses and because the application of Neyman-Pearson lemma is beautiful. However, the methodology from Li et al. (2019) can also be generalized for the input-dependent RS. To be able to do it, we need some auxiliary statements about the Rényi divergence.
Lemma 10. The Rényi divergence between two one-dimensional normal distributions is as follows:
Dα(N (µ1, σ21)||N (µ0, σ20)) = α(µ1 − µ2)2
2σ2α +
1
1− α log
( σα
σ1−α1 σ α 0
) ,
provided, that σ2α := (1− α)σ21 + ασ20 ≥ 0.
Proof. See Van Erven & Harremos (2014).
Note, that this proposition induces some assumptions on how σ0, σ1, α should be related. If σ0 > σ1, then the required inequality holds for any 1 6= α > 0. If σ0 < σ1, then α is restricted and we need to keep that in mind.
Lemma 11. Assume, we have some one-dimensional distributions P1,P1, . . . ,PN and Q1,Q2, . . . ,QN defined on common space for pairs with the same index. Then, assuming product space with product σ-algebra, we have the following identity:
Dα(P1 × P2 × · · · × PN ||Q1 ×Q2 × · · · × QN ) = N∑ i=1 Dα(Pi||Qi).
Proof. See Van Erven & Harremos (2014).
Using these two propositions, we are now able to derive a formula for Rényi divergence between two multivariate isotropic normal distributions:
Lemma 12.
Dα(N (x1, σ21I)||N (x0, σ20I)) = α‖x0 − x1‖2
2σ21 + 2α(σ 2 0 − σ11)
+N log ( σα σ1 ) 1− α −N α 1− α log ( σ0 σ1 ) .
Proof. Imporant property that is needed here is, that isotropic gaussian distributions factorize to one-dimensinal independent marignals. In other words:
N (x1, σ21I) = N (x11, σ21)×N (x12, σ21)× · · · × N (x1N , σ21),
and analogically for x0. Therefore, using Lemma 11 we see:
Dα(N (x1, σ21I)||N (x0, σ20I)) = N∑ i=1 Dα(N (x1i, σ21)||N (x0i, σ20)).
Now, it suffices to plug in the formula from Proposition 10 to obtain the required result:
Dα(N (x1i, σ21I)||N (x0i, σ20I)) = α(x1i − x2i)2
2σ2α +
1
1− α log
( σα
σ1−α1 σ α 0
)
= α(x1i − x2i)2
2σ21 + 2α(σ 2 0 − σ11)
+ log ( σα σ1 ) 1− α − α 1− α log ( σ0 σ1 ) Now it suffices to sum up over i and the result follows.
To obtain the certified radius, we also need a result from Li et al. (2019), which gives a guarantee that two measures on the set of classes will share the modus if the Rnyi divergence between them is small enough.
Lemma 13. Let P = (p1, p2, . . . , pK) and Q = (q1, q2, . . . , qK) two discrete measures on C. Let pA, pB correspond to two biggest probabilities in distribution P. Let M1(a, b) = a+b2 and M1−α(a, b) = ( a1−α+b1−α 2 ) 1 1−α If
Dα(Q||P) ≤ − log(1− 2M1(pA, pB) + 2M1−α(pA, pB)),
then the distributions P and Q agree on the class with maximal assigned probability.
Proof. This lemma can be proved by directly computing the minimal required Dα to be able to disagree on the maximal class probabilities via a constrained optimization problem (with variables pi, qi, i ∈ {1, . . . ,K}), solving KKT conditions. For details, consult Li et al. (2019).
Having explicit formula for the Rnyi divergence, we can mimic the methodology of Li et al. (2019) to obtain the certified radius:
Theorem 14. Given x0, pA, pB , σ0, N , the certified radius squared for all x1 such that fixed σ1 is used is:
R2 = sup α∈Sσ0,σ1
2σ21 + 2α(σ 2 0 − σ11)
α
( N α
1− α log ( σ0 σ1 ) −N
log ( σα σ1 ) 1− α
− log(1− 2M1(pA, pB) + 2M1−α(pA, pB)) ) ,
where Sσ0,σ1 = R+, if σ0 > σ1 and Sσ0,σ1 = ( 0, σ21
σ21−σ20
] if σ0 < σ1.
Proof. Let us fix x1 and assume, that α ∈ Sσ0,σ1 . Then, due to post-processing inequality for Renyi divergence, it follows that
Dα(f(x1 +N (0, σ21I))||f(x0 +N (0, σ20I))) ≤ Dα(x1 +N (0, σ21I)||x0 +N (0, σ20I))
= α‖x0 − x1‖2
2σ21 + 2α(σ 2 0 − σ11)
+N log ( σα σ1 ) 1− α −N α 1− α log ( σ0 σ1 ) .
Due to Lemma 13, it suffices that the following inequality holds for some α ∈ Sσ0,σ1 :
α‖x0 − x1‖2
2σ21 + 2α(σ 2 0 − σ11)
+N log ( σα σ1 ) 1− α −N α 1− α log ( σ0 σ1 ) ≤
− log(1− 2M1(pA, pB) + 2M1−α(pA, pB)).
This can be rewritten w.r.t. ‖x0 − x1‖2:
‖x0 − x1‖2 ≥ 2σ21 + 2α(σ 2 0 − σ11)
α
( N α
1− α log ( σ0 σ1 ) −N
log ( σα σ1 ) 1− α
− log(1− 2M1(pA, pB) + 2M1−α(pA, pB)) ) .
The resulting certified radius squared is now simply obtained by taking maximum over ‖x0 − x1‖2 s.t. ∃α ∈ Sσ0,σ1 such that the preceding inequality holds.
Note, that this theorem is formulated assuming, that except in x0, we use σ1 everywhere. It would require some further work to generalize this for general σ(x) functions, but to demonstrate the next point, it is not even necessary. Looking at the expression, we can observe that
N α
1− α log ( σ0 σ1 ) −N
log ( σα σ1 ) 1− α
depends highly on N and even for a ratio of σ0σ1 close to 1, we already obtain very strong negative values for high dimensions. The expression log(1 − 2M1(pA, pB) + 2M1−α(pA, pB)) is far less sensitive w.r.t pA and for large dimensions of N it is easily “beaten” by the first expression. Therefore, the higher the dimension N is, the bigger pA or the closer to 1 the σ0σ1 has to be in order to obtain even valid certified radius (not to speak about big). This points out that also the method of Li et al. (2019) suffers from the curse of dimensionality, as we know it must have done. This method is not useful for big N , because the conditions on pA, σ0, σ1 are so extreme, that barely any inputs would yield a positive certified radius. This fact is depicted in the Figure 11.
The key reason why this happens if done via Rényi divergences is that while the divergence Dα(N (x1, σ21I)||N (x0, σ20I)) grows independently of dimension as ‖x0−x1‖ grows, it drastically increases for big N even if x1 = x0! This reflects the effect, that if σ0 6= σ1, then the more dimensions we have, the more dissimilar are N (x1, σ21I) and N (x0, σ20I)). We can think of it as a consequence of standard fact from statistics that the more data we have, the more confident statistics against the null hypothesis σ0 = σ1 will we get if the null hypothesis is false. Since isotropic normal distributions can be actually treated as a sample of one-dimensional normal distributions, this is in accordance with our multivariate distributions setting.

heading: B.2 THE EXPLANATION OF THE CURSE OF DIMENSIONALITY
text: In the Section 2 we show that input-dependent RS suffers from the curse of dimenisonality. Now we will elaborate a bit more on this phenomenon and try to explain why it occurs. First, it is obvious from the Subsection B.1, that also the generalized method of Li et al. (2019) suffers from the curse of dimensionality, because the Rényi divergence between two isotropic Gaussians with different variances grows considerably with respect to dimension. This suggests that the input-dependent RS might suffer from the curse of dimensionality in general. To motivate this idea even further, we present this easy observation:
Theorem 15. Denote RC to be a certified radius given for pA and σ0 at x0 assuming the constant σ0 and following the certification of Cohen et al. (2019) 1. Assume, that we do the certification for each x1 by assuming the worst case-classifier as in Theorem 2. Then, for any x0, any function σ(x) and any pA, the following inequality holds:
R ≤ RC
Proof. Fix x1 and σ1. From Theorem 2 we know that the worst-case classifier f∗ defines a ball B such that P0(B) = 1− pA. From this it obviously follows, that the linear classifier fl and the linear spaceBl that assume constant σ0 also for x1 and is the worst-case for σ0 such that P0(Bl) = 1−pA is not worst-case for the case of using σ1 instead. Therefore, P1(Bl) ≤ P1(B). Moreover, let PC1 be a probability measure corresponding to N (x1, σ0I), i.e. the probability measure assuming constant σ0. It is easy to see that PC1 (Bl) > 0.5 ⇐⇒ P1(Bl) > 0.5 because the probability of a linear half-space under isotropic normal distribution is bigger than half if and only if the mean is contained in the half-space.
Assume, for contradiction that R > RC . From that, it exists a particular x1 such that PC1 (Bl) > 0.5 > P1(B), because otherwise there would be no such point, which would cause R > RC . However, PC1 (Bl) > 0.5 =⇒ P1(Bl) > 0.5, thus P1(Bl) > P1(B) and that is contradiction.
This theorem shows, that we can never achieve a better certified radius at x0 using σ0 and having probability pA than that, which we would get by Cohen et al. (2019)’s certification. Of course, this does not mean, that using non-constant σ is useless, since σ0 can vary. The question is, how much do we lose using non-constant σ. To get a better intuition, we plot the functions ξ< and ξ> under different setups in Figure 12, together with P1(Bl) from the proof of Theorem 15. From the top row we can deduce that dimensionN has a very significant impact on the probabilities and therefore also on the certified radius. We particularly point out the fact, that even ξ>(0), ξ<(0) can have significant margin w.r.t. to the probability coming out of linear classifier.2 Already for N = 90, we are not able to certify pA = 0.99 for rather conservative value of σ0σ1 . From middle row we see, that decreasing σ0 σ1
can mitigate this effect strongly. For instance, for σ0 = 1, σ1 = 0.95 the difference between P1(B) and P1(Bl) is almost negotiated. Bottom row compares ξ>(a), ξ<(a) and the respective linear classifier probabilities. We can see, that the case σ0 < σ1 might cause stronger restrictions on our certification (yet we deduce it just form the picture).
1The “C” in the subscript of certified radius might come both from “constant” and “Cohen et. al.” 2Notice the similarity with Rényi divergence, which also has positive value even for x0 = x1 if σ0 6= σ1
and then grows rather reasonably with distance.
What is the reason for ξ>(a), ξ<(a) being so big even at 0? The problem is following: Assume σ0 > σ1. If x0 = x1, the worst-case classifier coming from Lemma 2 will be a ball B centered right at x0, such that P0(B) = 1− pA. If we look at P1(B), we see, that we have the same ball centered directly at the mean, but the variance of the distribution is smaller. Using spherical symmetry of the isotropic gaussian distribution, this is equivalent to evaluating the probability of a bigger ball. If we fix σ0σ1 and look at the ratio of probabilities P1(B) P0(B) with increasing N , the curse of dimensionality comes into the game. For N = 2, the ratio is not too big. However, if N = 3072, like in CIFAR10, this ratio is far bigger. This can be intuitively seen from a property of chi-square distribution (which is present in the case x0 = x1), that while expectation is N , the standard deviation is “just” √ 2N ,
i.e. √ V ar(χ2N )
E(χ2N ) −→ 0 as N −→∞.
B.3 WHY DOES THE INPUT-DEPENDENT SMOOTHING WORK BETTER FOR SMALL σ VALUES?
As can be observed in Section 4 and Appendix E, the bigger the σb = σ we use, the harder it is to keep up to standards of constant smoothing. An interesting question is, why is the usage of small σb = σ helpful for the input-dependent smoothing?
Assume fixed σ, say σb = σ = 0.12. The theoretical bound on the certified radius given 100000 Monte-Carlo samplings and 0.001 confidence level using constant smoothing is about 0.48. Having
σ(x) ∼ 0.12, we cannot expect much bigger certified radius. Therefore, if we follow Theorem 7, the values of exp(−rR) and exp(rR) in the critical distance ∼ 0.5 will be much closer to 1, than the values of exp(−rR) and exp(rR) if we used σb = σ = 0.50 instead, where the critical values of R could be much bigger than 0.5. Therefore, the “gain” in P1(B) imposed by the curse of dimensionality, compared to P1(B) assuming constant σ will not be that severe yet. This means, that the loss in certified radius caused by the curse of dimensionality will be much less pronounced on the “active” range of certified radiuses (those for which the constant smoothing still works), compared to using big σb = σ. To support this idea, we demonstrate it on Figure 13, where we depict the certified radius as a function of distance from decision boundary, assuming f to be a linear classifier, using σb = σ = 0.12 and σb = σ = 0.50 for comparison.

heading: B.4 HOW DOES THE CURSE OF DIMENSIONALITY AFFECT THE TOTAL POSSIBLE
text: VARIABILITY OF σ(x)?
Fix certain type of task, say RGB image classification with images of similar object, but consider many possible resolutions (dimensions N ). Given two random images from the test set (x0, x1), what is the biggest reasonable value of |σ(x0)/σ(x1)−1|? Theoretically, the expression is bounded by | exp(±r‖x0 − x1‖) − 1|, given that r is the semi-elasticity constant of σ(x). However, the average distance between two samples from a test set of constant size, but increasing dimension scales as √ N . Therefore, with constant r, this upper-bound increases.
The increasing distance between samples is, therefore, a countereffect to the curse of dimensionality. In simple words, we have “more distance to change σ(x0) to σ(x1)”. Even if the r decreased just as 1/ √ N , the increasing distances would cancel the effect of the curse of dimensionaity and as a result, the maximal reasonable value of |σ(x0)/σ(x1)− 1| would remain roughly constant w.r.t. N . However, we need to take into account another effect. As the dimension increases, also the average distance of samples from the decision boundary increases. This is because the distances in general grow with dimension and if we assume that the number of intersections of a line segment between x0 and x1 with the decision boundary of the network remains roughly constant then the average distance from the decision boundary grows as √ N too. In order to compensate for this, we need to adjust the basic level of σ(x) (which we later call σb and can be understood as the general offset of our σ(x)) as √ N too. This is because the maximal attainable certified radius given fixed confidence level α and the number of Monte-Carlo samples is a constant multiple of σ(x).
However, with increased σ, we need to decrease the semi-elasticity rate r in order to obtain full certifications (see also Appendix B.3 for intuition behind this).
As a sketch of proof, we provide a simple computation, which tells us the approximate asymptotic behavior of |σ(x0)/σ(x1)− 1|. By Theorem 5 it holds:
exp(−rc √ N) ≥ √ 1− 2 √ − log(pB)
N ,
if we want to be able to predict a certified radius of c √ N (though this is just a necessary condition. For sufficiency, the LHS must be much closer to 1). After simple manipulation, we obtain:
r ≤ − 1 2c √ N log
( 1− 2 √ − log(pB)
N
) ∼ − 1
2c √ N
(−1)2 √ − log(pB)
N = √ − log(pB) cN .
So the rate scales as 1/N . Now we have:
|σ(x0)/σ(x1)− 1| ≤ | exp(±r‖x0 − x1‖)− 1| ≤ exp(r‖x0 − x1‖)− 1 ≤
exp( √ − log(pB) cN C √ N)− 1 ∼ √ −C log(pB) c √ N .

heading: B.5 DOES THE CURSE OF DIMENSIONALITY APPLY IN MULTI-CLASS REGIME?
text: In the main text, we presented a setup, where pB is set to be 1−pA. This is equivalent to pretending that we have just 2 classes. By not estimating the proper value of pB we lose some amount of power and the resulting certified radius is smaller than it could have been, did we have the pB as well. This is most pronounced for datasets with many classes. The natural question, therefore, is, whether we could avoid the curse of dimensionality by properly estimating the pB together with pA. The answer is no. The problem is that the the theory in Section 2 already implicitly works with the estimate of pB in a form of 1 − pA. The theory would work also with any other estimate of pB . Assuming constant pB , instead of constant pA, as we did in Section 2, will, therefore, yield the same conclusions. Moreover, there is neither theoretical, nor practical reason, why should pB decrease with increasing dimension.
This insight even applies to the question of the usage of input-dependent RS in practice. The assumption pB = 1−pA is no more important in Section 3 than in Section 2. Therefore, we can apply our method also for the pB obtained directly by Monte-Carlo sampling for the class B (or by any other estimation method).

heading: C CONCURRENT WORK
text: As we mention in Section 1, the idea to use input-dependent RS is not new. It has popped out in years 2020 and 2021 in at least four works from three completely distinct groups of authors, even though none of these works has been successfully published yet. We find it necessary to comment on all of these works because of two orthogonal reasons. First, it is a good practice to compare our work with the concurrent work to see what are pros and cons of these similar approaches and to what extend the approaches differ. Second, we are convinced, that three of these four works claim results, which are not mathematically valid. We find this to be a particularly critical problem in a domain such as certifiable robustness, which is by definition based on rigorous, mathematical certifications.

heading: C.1 THE WORK OF WANG ET AL. (2021)
text: This work, submitted for the ICLR conference 2021 is the only work that seems to be mathematically functional. In this work, authors have two main contributions – first, they propose a two-phase training, where in the second phase, for each sample xi, roughly the optimal σi is being found and then this sample xi is being augmented with this σi as an augmentation standard deviation. Authors call this method pretrain to finetune. Second, they provide a specific version of input-dependent RS. Essentially, they try to overcome the mathematical problems connected to the usage of nonconstant σ(x) by splitting the input space in so called robust regions Ri, where the constant σi is guaranteed to be used. All the certified balls are guaranteed to lie within just one of these robust regions, making sure that within one certified region, constant level of σ is used. Authors test this method on CIFAR10 and MNIST and show, that the method can outperform existing state-of-the-art approaches, mainly on the more complex CIFAR10 dataset.
However, we make several points, which make the results of this work, as well as the proposed method less impressive:
• The computational complexity of both their train-time and test-time algorithms seems to be quite high.
• The final smoothed classifier depends on the order of the incoming samples. As a consequence, it is not clear, whether the method works well for any permutation of the wouldbe tested samples. This creates another adversarial attack possibility - to attack the final smoothed classifier by manipulating the test set so that the order of samples is inappropriate for the good functionality of the final smoothed classifier.
• Even more, the fact, that the smoothed classifier depends on the order of the would-be tested samples makes it necessary, that the same smoothed classifier is used all the time for some test session in a real-world applications. For instance, a camera recognizing faces to approve an entry to a high-security building would need to keep the same model for its whole functional life, because restarting the model would enable attackers to create attacks on the predictions from the previous session. This might lead to significant restrictions on the practical usability of this method.

heading: C.2 THE WORKS OF ALFARRA ET AL. (2020) AND EIRAS ET AL. (2021)
text: In these works, similarly as in the work of Chen et al. (2021), authors suggest to optimize in each test point x for such a σ(x), that maximizes the certified radius given by Zhai et al. (2020), which is an extension of Cohen et al. (2019)’s certified radius for soft smoothing. The optimization for σ(x) differs but is similar in some respect (as will be discussed).
Besides, all three works further propose input-dependent training procedure, for which σ(x) - the standard deviation of gaussian data augmentation is also optimized. Altogether, both authors claim strong improvements over all the previous impactful works like Cohen et al. (2019); Zhai et al. (2020); Salman et al. (2019). The only significant difference between the works of Alfarra et al. (2020) and Eiras et al. (2021) (which have strong author intersections) is that in Eiras et al. (2021), authors build upon Alfarra et al. (2020)’s work and move from the isotropic smoothing to the smoothing with some specific anisotropic distributions.
As mentioned, authors first deviate from the setup of Cohen et al. (2019) and turn to the setup introduced by Zhai et al. (2020), i.e. they use soft smoothed classifier G defined as
GF (x)C = Eδ∼N (0,σ2I)F (x+ δ)C .
The key property of soft smoothed classifiers is that the Cohen et al. (2019)’s result on certified radius holds for them too. Theorem 16 (certified radius for soft smoothed classifiers). Let G be the soft smoothed probability predictor. Let x be s.t. G(x)A ≥ EA ≥ EB ≥ G(x)B . Then, the smoothed classifier g is robust at x with radius
R = σ
2 (Φ−1(EA)− Φ−1(EB)) = σ
Φ−1(EA) + Φ −1(1− EB))
2 ,
where Φ−1 denotes the quantile function of standard normal distribution.
Proof. Is provided in Zhai et al. (2020).
Note, that it is, similarly as in the hard randomized smoothing version of this theorem, essential to provide lower and upper confidence bounds for G(x)A and G(x)B , otherwise we cannot use this theorem with the required probability that the certified radius is valid. Denote G(x, σ) to be the soft smoothed classifier using σ in x. Authors propose to use the following theoretical σ(x) function:
σ(x) = arg max σ>0
σ 2 (Φ−1(G(x, σ(x))A)− Φ−1(G(x, σ(x))B)). (2)
It is of course not possible to optimize for this particular function since it is not known. It is also not feasible to run the Monte-Carlo sampling for each σ, because that is too costly and moreover due
Algorithm 1 Data dependent certification (Alfarra et al., 2020) function OPTIMIZESIGMA(F, x, β, σ0,M,K):
for k = 0, . . . ,K do sample δ1, . . . , δM ∼ N (0, I)
φ(σk) = 1 M M∑ i=1 F (x+ σδi)
ÊA(σk) = maxC φ(σk)C ÊB(σk) = maxC 6=A φ(σk)C R(σk) = σk 2 (Φ
−1(ÊA(σk))− Φ−1(ÊB(σk))) σk+1 ← σk + β∇σkR(σk)
σ∗ = σK return σ∗
to stochasticity, it would lead to discontinuous function. Treatment of this problem is probably the most pronounced difference between the works of Alfarra et al. (2020) and Chen et al. (2021).
Alfarra et al. (2020) use the following easy observation: N (0, σ2I) ≡ σN (0, I). Assume we have δi, i ∈ {1, . . . ,M} be i.i.d. sample from N (0, I). Obviously, G(x, σ(x))A ∼ 1M M∑ i=1 F (x+ σδi)A, since this is just the empirical mean of the theoretical expectation. Then, Expression 2 can be approximated as:
σ(x) = arg max σ>0
σ
2
( Φ−1 ( 1
M M∑ i=1 F (x+ σδi)A
) − Φ−1 ( 1
M M∑ i=1 F (x+ σδi)B
)) . (3)
Here, M is the number of Monte-Carlo samplings used to approximate this function. Note, that this function is a random realization of stochastic process in σ which is driven by the stochasticity in the sample δi, i ∈ {1, . . . ,M}. To find the maximum of this function, authors furhter propose to use simple gradient ascent, which is possible due to the simple differentiable form of Expression 3. This differentiability is one of the main motivations to switch from hard to soft randomized smoothing. Now, we are able to state the exact optimization algorithm of Alfarra et al. (2020):
Note, that being done in this way, this algorithm can be viewed as a stochastic gradient ascent. After obtaining σ∗ ≡ σ(x), authors further run the Monte-Carlo sampling to estimate the certified radius exactly as in Cohen et al. (2019), but with σ(x) instead of some global σ. Using this algorithm, authors achieve significant improvement over the Cohen et al. (2019)’s results, particularly getting rid of the first problem mentioned in Appendix A, the truncation issue. For the results, we refer to Alfarra et al. (2020). We will now give several comments on this algorithm and this method.
To begin with, in this optimization, authors do not adjust the estimated expectations and therefore don’t use lower confidence bounds, but rather raw estimates. This is not incorrect, since these estimates are not used directly for the estimation of certified radius, but it is inconsistent with the resulting estimation. In other words, authors optimize for a slightly different function than they then use. The difference is, however, not very big apart from extreme values of EA, where the difference might be really significant.
To overcome slightly this inconsistence, authors further (without comment) use clamping of the ÊA(σk) and ÊB(σk) on the interval [0.02, 0.98]. I.e. if ÊA(σk) > 0.98, it will be set to 0.98 and this is also taken into account in the computation of gradients. This way, authors get rid of the inconvenient issue, that if G(x)A ∼ 1, then ÊA(σk) ∼ 1 for σk ∼ 0, what might cause very big value of Φ−1(ÊA(σk)), yielding strong inconsistency with what would be obtained, if lower confidence bound was used instead.
However, the clamping causes even stronger inconsistence in the end. Note, that if G(x)A ∼ 1, then the true value of EA(σk) would be really close to 1, yielding high values of Φ−1(EA(σk)). This value would be far better approximated by the lower confidence bound than with the clamping, since
the lower confidence bound of 1 for M = 100000 and α = 0.001 is more than 0.9999, while the clamped value is just 0.98. This makes small values of σ highly disadvantageous, since σ2 −→ 0 as σ −→ 0, yet Φ−1(ÊA(σk)) is being stuck on Φ−1(0.98). In other words, this way authors artificially force the resulting σ(x) to be big enough, s.t. E(σ(x))A ≤ 0.98. This assumption is not commented in the article and might result in intransparent behaviour.
Second of all, authors use M = 1 for their experiments. This can be interpreted as using batch size 1 in classical SGD. We suppose that this small batch size is suboptimal since it yields an insanely high variance of the gradient.
Third of all, during the search for σ(x), it is not taken into account, whether the prediction is correct or not. This is, of course, a scientifically correct approach, since we cannot look at the label of the test sample before the very final evaluation. However, it is also problematic, since the function in Expression 2 might attain its optimum in such a σ(x), which leads to misclassification. This could have been avoided if constant σ was used instead.
To further illustrate this issue, assume F (x) = 1(B1(0)), i.e. F (x) predicts class 1 if and only if ‖x‖ ≤ 1, otherwise predicts class 0. Assume we are certifying x ≡ 0 and assume that σ0 in Algorithm 1 is initialized such that class 0 is already dominating. Then, we will have positive gradient ∇σkR(σk) in all steps, because F (σδi) is obviously non-increasing, so the number of points classified as class 1 for fixed sample δi, i ∈ 1, . . . ,M is decreasing, yielding ÊA(σk) nondecreasing in σk, while σk2 strictly increasing in σk. This way, the σk will diverge to∞ for k −→∞. However, point x ≡ 0 is classified as class 1, yielding misclassification which is, moreover, assigned very high certified radius.
This issue is actually even more general - the function in Expression 2 does in most cases (assuming infinite region RN ) not possess global maximum, because usually
lim σ−→∞
σ 2 (Φ−1(G(x, σ(x))A)− Φ−1(G(x, σ(x))B)) =∞.
This can be seen, for instance, easily for the F (x) = 1(B1(0)), but it is the case for any hard classifier, for which one region becomes to have dominating area as the radius around some x0 goes to infinity. This is because, if some region becomes to be dominating (for instance if all other regions are bounded), then σ2 grows, while Φ
−1(G(x, σ(x))A) − Φ−1(G(x, σ(x))B) either grows too, or stagnates, making the whole function strictly increasing with sufficiently high slope.
This issue also throws the hyperparameter K under closer inspection. What is the effect of this hyperparameter on the performance of the algorithm? From the previous paragraph, it seems, that this parameter serves not only as the “ scaled number of epochs”, but also as some stability parameter, which, however, does not have theoretical, but rather practical justification.
Another issue is, that the function in Expression 2 might be non-convex and might possess many different local minima, from which not all (or rather just a few) are actually reasonable. Therefore, the Algorithm 1 is very sensitive to initialization σ0.
However, probably the biggest issue of all is connected to the impossibility result showed in Section 2, which shows, that the Algorithm 1 actually yields invalid certified radiuses. Why it is so?
First of all, we must justify, that our impossibility result is applicable also for the soft randomized smoothing. This is because classifiers of type F (x)C = 1(x ∈ RC) forRC being decision region for class C are among applicable classifiers s.t. G(x, σ)A = EA. With such classifiers, however, there is no difference between soft and hard smoothing and moreover EA ≡ pA from our setup. This way we can construct the worst-case classifiers F ∗ exactly as in our setup and therefore the same worstcase classifiers and subsequent adversarial examples are applicable here as well. In other words, for fixed value of soft smoothed G(x, σ)A = EA we can denote pA = EA and find the worst-case hard classifier F defined as indicator of the worst-case ball, which will yield EB ≡ P1(B) from Theorem 3 in some queried point x1.
As we have seen in previous paragraphs, the resulting σ(x) yielded in Algorithm 1 is very instable and stochastic - it depends heavily on F, σ0,K, β,M and of course δi, i ∈ {1, . . . ,M} for each iteration of the for cycle. Now, for instance for CIFAR10 and pA = 0.99, we have the minimal possible ratio σ0σ1 equal to more than 0.96. It is hard to believe, that such instable, highly stochastic and non-regularized (except for K,β) method will yield σ(x) sufficiently slowly varying such that
within the certified radius around x0, there will be no x1 for which σ1 deviates more than by this strict threshold from σ0. This is even more pronounced on ImageNet, where the minimal possible ratio σ1σ0 is above 0.99 for any pA or EA.
Even without the help of curse of dimensionality, we can construct a counterexample for which the algorithm will not yield valid certified radius. Assume again F (x) = 1(B1(0)) and assume modest dimension N = 2. Assume we try to certify point x0 ≡ [50, 0]. Then, the theoretical σ-dependent function from Equation 2 is depicted on Figure 14.
We can see, that the resulting σ(x0) will be as big as our regularizers K and β in Algorithm 1 will allow. Therefore, if we run the algorithm for K high-enough, surely the resulting certified radius will be far bigger than 50. However, if we certify the point x0 ≡ [0, 0] and we start with σ0 = 0.25, for instance, then the σ-dependent certified radius in Expression 2 will be decreasing in this σ0, yielding σ(x0) < 0.25, which will result in classification of class 1. This point [0, 0] lies within the “certified” range of [50, 0], yet it is not classified the same, because, obviously, [50, 0] is classified as class 0. This is therefore a counterexample to the validity of Alfarra et al. (2020)’s certification method and their results.
Note that even though our counterexample is a bit “extreme” and one could argue that in practice such a situation would not occur, we must emphasize, that this counterexample is constructed even without the help of the curse of dimensionality. In practice, it fully suffices, that for x0 and some certified radius R in x0, there exists x1 within the range of this certified radius, s.t. σ1 is quite dissimilar to σ0. If such situation occurs, then R surely is not a valid certified radius.

heading: C.3 THE WORK OF CHEN ET AL. (2021)
text: The methodology of Chen et al. (2021) is rather similar to that of Alfarra et al. (2020). The biggest difference consists in the optimization of Expression 2.
Instead of stochastic gradient descent, they use more sophisticated version of grid search - so called multiple-start fast gradient sign search. Simply speaking, this method first generates a set of pairs (σ0, s)i, i ∈ {1, . . . ,K} and then for each of the i runs, it runs a j-while cycle, where in each step j, it increases σ2j = σ 2 0 + js to σ 2 j+1 = σ 2 0 + (j + 1)s and checks, whether the σ-dependent empirical certified radius in Expression 3 increases or not. If yes, they continue until j is above some threshold T , if not, they break and report σi as the σj from the inner step where while cycle was broken. After obtaining σi for i ∈ {1, . . . ,K}, they choose σ(x) to be the one, that maximizes Expression 3. More concretely, their multiple-start fast gradient sign search algorithm looks as follows:
It is not entirely clear from the text of Chen et al. (2021), how exactly are lk, sk sampled, but it is written there, that the interval for l is [1, 16] and for s it is (−1, 1). Moreover, the authors don’t
Algorithm 2 Instance-wise multiple-start FGSS (Chen et al., 2021) function OPTIMIZESIGMA(F, x, σ0,M,K, T ):
generate (lk, sk), k ∈ {1, . . . ,K} for k = 1, . . . ,K do
sample δ1, . . . , δM ∼ N (0, lkσ20I) φ( √ lkσ0) = 1 M M∑ i=1 F (x+ δi) ÊA( √ lkσ0) = maxC φ( √ lkσ0)C R( √ lkσ0) = √ lkσ0Φ −1(ÊA( √ lkσ0))
mk = R( √ lkσ0) while lk ∈ [1, T ] do sample δ1, . . . , δM ∼ N (0, (lk + sk)σ20I)
φ( √ lk + skσ0) = 1 M M∑ i=1 F (x+ δi) ÊA( √ lk + skσ0) = maxC φ( √ lk + skσ0)C R( √ lk + skσ0) = √ lk + skσ0Φ −1(ÊA( √ lk + skσ0))
if R( √ lk + skσ0) ≥ R( √ lkσ0) then
lk ← lk + sk mk = R( √ lk + skσ0)
else break
σ(x) = max k∈{1,...,K} mk return σ(x)
provide the code and from the text, it seems, that they don’t use lower confidence bounds during the evaluation of certified radiuses, what we consider to be a serious mistake (if really the case). However, we add some comments to this method regardless of the lower confidence bounds.
Generally, this method possesses most of the disadvantages mentioned in Section C.2. They use the same function for optimization, the Expression 2 and its empirical version 3. This means, that the method suffers from having several local optima, having no global optimum in general (and in most cases with limit infinity). Similarly like before, here is also no control over the correctness of the prediction, i.e. many or all local optima might lead to misclassification.
On the other hand, in this paper authors use M = 500 (the effective batch size), which is definitely more reasonable than M = 1 as in Alfarra et al. (2020). Furthermore, they use multiple initializations, making the optimization more robust and improving the chances to obtain global, or at least very good local minimum.
However, the main problem, the curse of dimensionality yielding invalid results is even more pronounced here. Unlike the “continuous approach” in Alfarra et al. (2020), here authors for each x0 sample just some discrete grid (more complex, since there are more initializations) of possible values of σ(x). For instance, if s = 1, then the smallest possible ratio between two consecutive l’s in the Algorithm 2 is √ 15/4 ∼ 0.97, making it impossible to certify some x1 w.r.t. x0 if for both s = 1 and l0 6= l1 on ImageNet and also for a lot of samples on CIFAR10. Of course, the fact that s is randomly sampled from (−1, 1) makes this counter-argumentation more difficult, but it is, again, highly unlikely that this highly stochastic method without control over σ(x) would yield function with sufficiently small semi-elasticity. Therefore also the impressive results of Chen et al. (2021) are, unfortunately, scientifically invalid.
D IMPLEMENTATION DETAILS
Even though our algorithm is rather easy, there are some perks that should be discussed before one can safely use it in practice. First, we show the actual Algorithm 3
Algorithm 3 Pseudocode for certification and prediction of my method based on Cohen et al. (2019) # evaluate g at x0 function PREDICT(f, σ0, x0, n, α):
counts←− SampleUnderNoise(f, x0, n, σ0) ĉA, ĉB ←− two top indices in counts nA, nB ←− counts(cA),counts(ĉB) if BinomPValue(nA, nA + nB , 0.5) ≤ α then return ĉA else return ABSTAIN
# certify the robustness of g around x0 function CERTIFY(f, σ0, x0, n0, n, α):
counts0←− SampleUnderNoise(f, x0, n0, σ0) ĉA ←− top index in counts0 counts←− SampleUnderNoise(f, x0, n, σ0) pA ←− LowerConfBound(counts[ĉA], n, 1− α) if pA > 1/2 then return prediction ĉA and radius
ComputeCertifiedRadius(σ0, r,N, pA, num steps) else return ABSTAIN
function COMPUTECERTIFIEDRADIUS(σ0, r,N, pA, num steps) radiuses←− linspace(num space) for R in radiuses do
σ11 ←− σ0 exp(−rR) σ12 ←− σ0 exp(rR) xi bigger← ξ>(R, σ11) xi lower← ξ<(R, σ12) if max{xi bigger,xi lower} > 0.5 then BREAK
return R
Note that the function ComputeCertifiedRadius is a bit more complicated than depicted in Algorithm 3. We don’t use a simple for-loop, but rather quite an efficient search method.
Theoretically speaking, this algorithm works perfectly. However, in practice, it is a bit problematic. The issue is, that since we use σ11 and σ12, which are extremely close to σ0 for small tested radiuses R, the NCCHSQ functions will get extremely high inputs, making the results numerically instable. To prevent this, we use a simple trick. Since the more extreme σ1 will we assume in evaluation at particular distance R, the worse for us, we can prevent numerical issues simply by putting σt < σ0 and σT > σ0 to be maximal and minimal used σ’s in our evaluation, i.e. the true σ used will be min{σt, σ0 exp(−rR)} and max{σT , σ0 exp(rR)}. This way, we avoid numerical issues, because we can put σt, σT to be s.t. 1σ20−σ2t is not too big and in the same time maintainting the correct certification thanks to the Lemma 6. The problems of this workarounds are first that it decreases the certification power, since it assumes σ1’s that are even worse than the theoretically guaranteed worst-case possibilities and second, more importantly, that it requires some engineering to design the σt, σT designs. It is submoptimal to put one constant value for these thresholds, because the numerical problems occur at different ratio thresholds of σ1/σ0 for different class probabilities pA and the dimension N . This requires to design a specific σt(pA) and σT (pA) functions for each dimension N which we want to apply. For instance, we use σt(pA) = 0.9993 + 0.001 log10(pB), and σT (pA) = 1/σt(pA) for CIFAR10, while for MNIST we use σt(pA) = 0.9988 + 0.001 log10(pB), and σT (pA) = 1/σt(pA). To design such functions, one needs to plot plot real probability of a ball with fixed variances as fcn of dist function, which computes the ξ functions, for particular N and several values of pA and look, whether it computes correctly. As an example, we provide such a plots for well and ill working setups on Figure 15.
Another performance trick is to not evaluate ξ for each Ri, where Ri is i-th grid point of evaluation, but rather evaluate sequentially Ri2 , i.e. just every i2-th point, until we reach value > 0.5 and then to search just the interval [(i− 1)2, i2], where i is the first iteration for which ξ>(Ri2 , σ1) ≥ 0.5.

heading: E MORE TO EXPERIEMNTS AND ABLATIONS
text: Before we present our further results, we must emphasize that our certification procedure is barely any slower than that of Cohen et al. (2019). More specifically, given 100000 iterations of montecarlo sampling, certification of one sample using Cohen et al. (2019)’s algorithm on CIFAR10 takes ∼ 15 seconds on our machine, while certification of a sample using our Algorithm 3 takes 15− 20 seconds depending on the σb, r setup. If at least one of σb and r is not small, then our method runs practically instantly. If both parameters are small, then one evaluation can take up to 5 seconds depending on the exact value of parameters and on the pA. Note, that this part of the certification is dimension-independent and therefore can run in the same time also on much higher-dimensional problems.
Besides the actual certification, we have to compute σ(x) for each of the test examples. This part of the algorithm is being executed before the actual certification and usually takes around 1 minute on our machine and on CIFAR10.
All in all, even in the really worst-case scenario, our method runs at most 1/3-times longer than the old method on CIFAR10. On MNIST, the ratio between our run and the original run is higher, since MNIST is smaller-dimensional problem. However, since our part of evaluation is practically independent of the setup (except the values of σb and r, which, however, can yield just some upperbounded amount of slow-down), our algorithm does not bring any added asymptotic time complexity.

heading: E.1 HOW TO CHOOSE THE HYPERPARAMETERS?
text: Our design of σ(x) function defined in Equation 1 uses several hyperparameters. These are: r for the rate, m for the scaling, σb for the base sigma and k for the k-nearest neighbors. How do we choose these hyperparameters?
The m parameter depends on our goals. We can set it so that σ(x) achieves lowest values at σb by setting it so that it is roughly equal to the minimal distance from k nearest neighbors across, for instance, training samples. Other possibility is to set it so that it is roughly equal to the average distance from k nearest neighbors, to ensure that the average σ(x) will roughly correspond to the σb.
The k parameter needs to be set with two objectives in mind. Firstly, it would be unwise to set it too small, because then the distance from k nearest neighbors would be too noisy. On the other hand, we don’t want it too high, because then it will not be changing fast enough with changing the position of x. The suitable value can be obtained by looking at histograms of average distances from k nearest neighbors and choosing the k for which the histogram is enough scattered, but it is not too small.
The r parameter needs to be chosen so that we can have some significant advantage over constant smoothing, but it cannot be too big, because otherwise the curse of dimensionality would apply. The value can be decided either by trial and error, or by plotting the certified radius given linear classifier, or from Theorem 4, setting the rate low-enough so that within the expected certified radius range, the ratio σ1σ0 can’t move anywhere near the theoretical thresholds implied by Theorem 4.
The σb is the base σ and should be used according to the level of smoothing variance we want to use. More discussion on this can be found in Cohen et al. (2019).

heading: E.2 COMPARISON WITH COHEN ET AL. (2019) METHODOLOGY ON CIFAR10 DATASETS
text: Here, we compare Cohen et al. (2019)’s evaluations for σ = 0.12, 0.25, 0.50 with our evaluations directly on models trained by Cohen et al. (2019), setting σb = σ, r = 0.005, 0.01 and 0.015 for σb = σ = 0.12, k = 20 and m = 5. In this way, the levels of σ(x) used in direct comparison will rise from the values roughly equal to Cohen et al. (2019)’s constant σ to higher values. The results are depicted in Figure 16.
Note, that this evaluation is being done on the models trained directly by Cohen et al. (2019) and therefore the variance of Gaussian data augmentation is not entirely consistent with the optimal variance that should be used for non-constant σ, which should be either the same, σ(x) or constant, but in average equal to σ(x). The results are similar as in the Section 4. Note, that for σb = σ = 0.50, the curse of dimensionality becomes most pronounced, as explained in Appendix B. Further, we provide the Tables 6, 7, where the clean accuracies and class-wise standard deviations are displayed.
The results are, again, similar as in the section 4.

heading: E.3 COMPARISON WITH COHEN ET AL. (2019) METHODOLOGY ON MNIST DATASETS
text: Here, we present similar comparison as in Subsection E.2, but on MNIST and with models trained by us. Again, the setup is similar as in the Section 4. We compare σ = σb = 0.12, 0.25, 0.50 with test-time rates r = 0.005, 0.01, 0.02, 0.05 and train-time level of σ again equal to σ = σb. It is important to note, that we use different normalization constant m in the MNIST case. In CIFAR10, we set m = 5, in MNIST, the suitable m is 1.5. This way we assure, that the smallest σ(x) values in the test set will roughly equal the σb = σ. The certified accuracy plots are depicted on Figure 17. We also add the clean accuracy table and class-wise clean accuracies standard deviation table (8, 9).
All the results are, again, very similar to those presented in Section 4, even though the gain in certified accuracies is marginally worse, since our evaluations run on models trained with in average smaller train-time data-augmentation standard deviation σb = σ.
E.4 INVESTIGATION OF THE EFFECT OF TRAINING WITH INPUT-DEPENDENT GAUSSIAN AUGMENATATION
It has been shown by many works, that apart from a good test-time certification method, also the appropriate training plays a very important role in the final robustness of our smoothed classifier g. Already Cohen et al. (2019) realize this and propose to train with gaussian data augmentation with constant σ. They experiment with different levels of σ during training and conclude that training with the same level of σ that will be later used in the test time is usually the most suitable option.
The question of best-possible training to boost the certified robustness didn’t stay without the interest of different researchers. Both Zhai et al. (2020) and Salman et al. (2019) try to improve the way of training and propose two different, yet interesting and effective training methods. While Zhai et al. (2020) manage to incorporate the adversarial robustness into the training loss function, therefore training directly for the robustness, Salman et al. (2019) propose to use adversarial training to achieve more robust classifiers.
Both Alfarra et al. (2020) and Chen et al. (2021) already propose to use training with input-dependent σ as the variance of gaussian data augmentation. Both of them proceed similarly as during test time - to obtain training σ(x), they optimize for such, that would maximize the certified accuracy of training samples.
In this section, we propose and test our own training method. We propose to use again gaussian data augmentation with input-dependent σ(x), but we suggest to use the simple σ(x) defined in Equation 1. In other words, we suggest using the same σ(x) during training as during testing (up to parametrization, which might differ).
Note, that, unlike the certification, the training procedure does not require any mathematical analysis nor certification. It is totally up to us how we train the base classifier f and the way of training does not influence the validity of subsequent certification guarantees during test time. However, it is good to have a reasonable training procedure, because otherwise, we would achieve a satisfactory model neither in terms of clean accuracy nor in terms of adversarial robustness.
In the subsequent analysis, we evaluate and compare our certification procedures on models trained with different training parametrizations. For this particular section, we run the comparison only on the CIFAR10 dataset. For each test-time σb, r, we evaluate our method with these parameters on base models f trained with the same σb, but different level of training rate trr. The training rate trr plays exactly the same role as the evaluation rate r but is used exclusively during training. Note, that this makes our σ(x) different during training and testing since it is parametrized with different rates.
On the Figure 18 we plot evaluations on CIFAR10 of our method for rate 0.01, all levels of σb = 0.12, 0.25, 0.50 and each of these test-time setups is evaluated on 4 different levels of train-time rate trr = 0.0, 0.01, 0.04, 0.1.
From the results, we judge, that our training procedure works satisfactorily well. It can generally outperform the constant σ training, yet the standard accuracy vs. robustness trade-off is present in some cases. If we train with small train-time rate, the improvement of the certified accuracies is not pronounced (the case for σb = σ = 0.50 is slightly misleading, since such a configuration is just a result of the variance of clean accuracy w.r.t different traning runs) enough, but we also don’t lose almost any clean accuracy. Increasing the rate to trr = 0.04 results in much more pronounced improvements in high certified accuracies, yet also comes at a prize of clean accuracy drop, especially for large σ levels. Even bigger training rate, such as trr = 0.1 seems to be too big and does not bring almost any improvement over the rate trr = 0.04, yet loses a large amount of clean accuracy.
These results suggest, that the input-dependent training with a carefully chosen training rate for σ(x) can lead to significant improvements in certifiable robustness. However, it is important to note, that
the optimal trr seems to be dependent on the σb, therefore for each value of σb, some effort has to be invested to find the optimal hyperparameters.
Besides, we were also interested, whether using an input-dependent σ(x) during training influences the class-wise accuracy balance. In Table 10 we report the standard deviations of class-wise accuracies.
We can observe, that unlike the pure input-dependent evaluation, the input-dependent training is partially capable of mitigating the effects of the shrinking. For instance, the trr = 0.04 for σb = 0.12 provides obvious improvement in establishing class-wise balance. Similarly successful are trainings with trr = 0.01 for σb = 0.12 and both trr = 0.01, 0.04 for σb = 0.25. Also for σ = 0.50 the mitigation is present for small-enough training rates. However, we must emphasize, that if we use too big training rate, the disbalance between class accuracies will be re-established and in some cases even magnified. Therefore, we must be careful to choose the appropriate training rate for the σb, r.

heading: E.5 WHY DO WE NOT COMPARE WITH THE CURRENT STATE-OF-THE-ART?
text: Briefly speaking – we could, but we don’t consider it necessary. Since we claim just one type of improvement over the Cohen et al. (2019)’s model (experiment-wise) and don’t claim new state-ofthe-art training method, we didn’t find it necessary to measure our strengths with methods of Salman et al. (2019) and Zhai et al. (2020). It is obvious that we would outperform these methods in the question of certified accuracy waterfalls anyway, since these methods focus on the training phase. Since we do not outperform Cohen et al. (2019) neither in terms of the clean accuracies nor in terms of class-wise accuracies, it is not our belief that we would outperform the two modern methods in these metrics. Moreover, we find the comparison with Cohen et al. (2019) most structured, since we extend the theory built by them.

heading: E.6 ABLATIONS
text: Even though our results so far might look impressive, we can’t claim that it is fully due to our particular method until we exclude the possibility, that some different effects play an essential role in the improvement over Cohen et al. (2019)’s work.
To investigate, whether our particular method dominates the contribution to the performance boost, we conduct several ablation studies - first, we study the variance of our evaluations and trainings, second, we study the effect of input-dependent test-time randomized smoothing, and third, we study the effect of input-dependent train-time data augmentation.
E.6.1 VARIANCE OF THE EVALUATION
To find out, whether there is a significant variance in the evaluation of certified radiuses, we conduct a simple experiment - we train a single model on CIFAR10 and evaluate our method on this model for the very same setup of parameters multiple times. This way, the only present stochasticity is in the Monte-Carlo sampling, which influences the evaluation of certified radiuses. We pick the parameters as follows: σb = 0.50, r = 0.01, trr = 0.0, since the σb = 0.50 turns out to have biggest variance in the training. The results are depicted in Figure 19.
From the results, it is obvious that the variance in the evaluation phase is absolutely negotiable. Therefore, there is no need to run the same evaluation setup more times.
E.6.2 VARIANCE OF THE TRAINING
To estimate the variance of the training, we train several models for one specific training setup and evaluate them with the same evaluation setup (knowing, that there is no variance in the evaluation phase, this is equivalent to measuring directly the training variance). We pick our classical non-constant σ(x) for the evaluation, but we train with constant variance data augmentation. The concrete parameters we work with are: σb ∈ {0.12, 0.25, 0.50}, r = 0.01, trr = 0.0 and we run 9 trainings for each of these parameter configurations. Then we run full certification to not only see the variance in clean accuracy, but also the variance in the certified radiuses. The results are depicted in Figure 20.
From the figures we see, that the variance of the training is strongly σb-dependent. Most volatile clean accuracy is present for the case σb = 0.50. However, fortunately, the biggest variability is present for the clean accuracy and the curves seem to be less scattered in the areas of high certified radiuses. The concrete standard deviations of clean accuracies are in Table 11. The standard deviations of clean accuracies for MNIST dataset and the same parameters are in Table 12.
Since the differences in accuracies of different methods are very subtle, it is hard to obtain statistically trustworthy results. For instance, given, that the standard deviation 0.4% is the true standard deviation of the σb = 0.25 runs, we would need 16 runs to decrease it to a standard deviation of 0.1%, which might be considered to be precise-enough. To do the same in the case of σb = 0.50 on CIFAR10, we would roughly need 400 runs to decrease the standard deviation below 0.1%. Therefore, the results we provide in the subsequent subsections, being the average of “just” 8 runs, have to be taken just modulo variance in the results, which might still be considerable.

heading: E.6.3 EFFECT OF INPUT-DEPENDENT EVALUATION
text: In this ablation study, we compare the certification method for particular σb, r = 0.01, trr = 0.0 with the constant-σ certification method with Cσb, r = 0.0, trr = 0.0, where C is an appropriate constant. The motivation behind such an experiment is, that our σ(x) is generally bigger than σb, but originally, we compare this method to constant σ = σb evaluation. Therefore, in average, samples in our method enjoy bigger values of σ(x). Natural question is, whether we cannot obtain the same performance boost using just the constant σ method with Cσb > σb set to such value, which roughly corresponds to the average of σ(xi) for xi, i ∈ {1, . . . , T} being the test set. The problem of using bigger Cσb is, that we encounter performance drop and more severe case of shrinking, but we need to check, to what extent is the performance drop present in the input-dependent σ(x) method. Comparing the performance drops of larger constant Cσb and input-dependent σ(x), which is in average larger (but in average the same as the Cσb), we will be able to answer, to what degree is the usage of input-dependent σ(x) really justified. If we remind ourselves, that
σ(x) = σb exp r 1 k  ∑ xi∈Nk(x) ‖x− xi‖ −m  ,
then we see, that the constant C we are searching for is the average (or rather median) value of
exp r 1 k  ∑ xi∈Nk(x) ‖x− xi‖ −m  .
Fortunately, empirically, the mean and median of the above expression are roughly equal for both CIFAR10 and MNIST, so we are not forced to choose between them. For r = 0.01,m = 5, we choose the rounded value of C = exp(0.05) on CIFAR10. For r = 0.01,m = 1.5 as in MNIST, the constant is set to C = 1.035. In the end, the values of Cσ used in this experiment are Cσ = 0.126, 0.263, 0.53 for CIFAR10 and Cσ = 0.124, 0.258, 0.517 for MNIST. To obtain a fair comparison, though, we evaluate the input-dependent σ(x) evaluation strategy on models trained with constant Cσb standard deviation of gaussian augmentation. This is because this level of σ is equal to the mean value of the σ(x) and we believe, that such a training data augmentation standard deviation is more consistent with our σ(x) function. We provide the plots of single-run evaluations of certified accuracies for CIFAR10 in Figure 21 and for MNIST in Figure 22. The models on which we evaluate differ because for the increased constant σ evaluations we needed to also use an increased level of the training data augmentation variance.
From the figures, it is obvious, that our method is not able to outperform the constant σ method using the same mean σ in terms of certified accuracy, not even for our strongest σb = 0.12. This fact might not be in general bad news, if we demonstrated, that our method suffers from less pronounced accuracy drop or less pronounced disbalance in class-wise accuracies. To find out, we measure
average accuracies of the evaluation strategies from 8 runs for each, as well as average class-wise accuracy standard deviations from 8 runs. The results are provided in Tables 13 and 14 for CIFAR10 and Tables 15 and 16 for MNIST.
As for CIFAR10, except for σ = σb = 0.25, the differences in accuracies between different evaluation strategies are so small, that we cannot consider them to be statistically significant. Even though the difference for σ = σb = 0.25 is high, it is still not possible to draw some definite conclusions, especially for the difference between the input-dependent σ(x) and the increased constant Cσb evaluations. In general, it is not easy to judge, whether our method possesses some advantage (or disadvantage) over the increased Cσb method in terms of clean accuracy. Similar conclusions can be drawn in the context of the shrinking phenomenon. Here, the differences are also very small, but
unlike in the comparison with Cohen et al. (2019) models, where we evaluate our input-dependent σ(x) method on classifiers trained with inconsistent data-augmentation variance, here we observe the general trend, that our method is able to outperform the increased constant Cσb evaluation. This is good news and it confirms our suspicion, that the bad results from Subsection E.2 could come from the train-test σ inconsistency.
The results on MNIST suggest similar conclusions for the accuracy vs. robustness tradeoff. Similarly, the σ = σb = 0.12, 0.50 are not telling much, and for σ = σb = 0.25, the differences are still rather small (yet the standard deviation of the results should be ∼ 0.0001, so it is rather on the edge). The conclusions for the shrinking phenomenon are a bit more pesimistic than in the case of CIFAR10. Here we don’t see any improvement over the constant σ, not even the one with increased σ level.

heading: E.6.4 EFFECT OF INPUT-DEPENDENT TRAINING
text: In this last ablation study, we compare our input-dependent data augmentation for particular σb, r and particular training rate trr with the constant Cσb data augmentation, where the training rate trr is set to 0. The strategy for choosing the constant C is exactly the same as in the first experiment. Particularly, we evaluate our method with r = 0.01 and σb = 0.12, 0.25, 0.50, trained with the same level of σb and training rate trr = 0.01 with the evaluations using r = 0.01 and σb = 0.12, 0.25, 0.50 during test time, while during train time using training rate trr = 0.0, but using the constant σ = 0.126, 0.263, 0.53 for CIFAR10 and σ = 0.124, 0.258, 0.517 for MNIST. This way, we compensate for the “increased levels of σ(x)” with respect to σb. We present our comparisons in the Figure 23 for CIFAR10 and 24 for MNIST, providing the evaluations with r = 0.01, σb = 0.12, 0.25, 0.50 and the same σb and trr = 0.0 during train time as a reference.
The certified accuracy results for the CIFAR10 and the MNIST differ a bit. For CIFAR10 training with rate r = 0.01 is does not overperform the constant Cσb training. For σb = 0.12, the constant Cσb training clearly outperforms the input-dependent training. For σb = 0.25, these two training strategies seem to have almost identical performances. For σb = 0.50, the input-dependent σ(x) strategy outperforms the constant σ ones, but now we know, that it is purely due to the variance in the training. On the other hand on MNIST, we either have very similar performance or even slightly outperform the constant σ training.
Looking at the accuracy and standard deviation Tables 17, 18, 19 and 20, we can deduce the following. In terms of clean accuracy, the input-dependent training strategy performs worst in most of the cases, even though the differences in performance might not be statistically significant. We see, that we would need far more evaluations to see some clear pattern. However, these results are definitely not good news for the use of input-dependent σ(x) during training.
In terms of the class-wise accuracy standard deviation, we again see countering results for CIFAR10 and MNIST datasets. For CIFAR10 the input-dependent σ(x) clearly outperforms the smaller constant σb training method, particularly for σb = 0.50. However, the constant Cσb method seem to outperform even the input-dependent σ(x). For MNIST, the smaller constant σb outperforms both other methods, while they are rather similar.
Together with findings from previous sections, these results suggest, that usage of this particular design of input-dependent σ(x) might not be worthy until a more precise evaluation is conducted. However, the combination of input-dependent test-time evaluation with constant, yet increased traintime augmentation is possibly the strongest combination that can be achieved using input-dependent sigma at all (especially for CIFAR10).

heading: F PROOFS
text: Lemma 17 (Neyman-Pearson). Let X,Y be random vectors in RN with densities x, y. Let h : RN −→ {0, 1} be a random or deterministic function. Then, the following two implications hold:
1. If S = { z ∈ RN : y(z)x(z) ≤ t } for some t > 0 and P(h(X) = 1) ≥ P(X ∈ S), then
P(h(Y ) = 1) ≥ P(Y ∈ S). 2. If S = { z ∈ RN : y(z)x(z) ≥ t } for some t > 0 and P(h(X) = 1) ≤ P(X ∈ S), then
P(h(Y ) = 1) ≤ P(Y ∈ S).
Proof. See Cohen et al. (2019).
Lemma 1: Out of all possible classifiers f such that Gf (x)B ≤ pB = 1 − pA, the one, for which Gf (x+ δ)B is maximized is the one, which predicts classB in a region determined by the likelihood ratio:
B = { x ∈ RN : f1(x)
f0(x) ≥ 1 r
} ,
where r is fixed, such that P0(B) = pB . Note, that we use B to denote both the class and the region of that class.
Proof. Let f be arbitrary classifier. To invoke the Neyman-Pearson Lemma 17, define h ≡ f (with the only difference, that h goes to {0, 1} instead of {A,B}). Moreover, let S ≡ B and X ∼ N (x, σ20I), Y ∼ N (x + δ, σ21I). Let also f∗ classify S as B. Then obviously, P(X ∈ S) = P0(B) = pB . Since Gf (x)B ≤ pB , we have P(h(X) = 1) ≤ pB . Using directly the second part of Neyman-Pearson Lemma 17, this will yield P(Y ∈ S) ≥ P(h(Y ) = 1). Rewritten in the words of our setup, Gf∗(x+ δ) ≥ Gf (x+ δ).
Theorem 2: If σ0 > σ1, then B is a N -dimensional ball with the center at S> and radius R>:
S> = x+ σ20
σ20 − σ21 δ, R> =
√ σ20σ 2 1
(σ20 − σ21)2 ‖δ‖2 + 2N σ
2 0σ 2 1
σ20 − σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 σ20 − σ21 log(r).
If σ0 < σ1, then B is the complement of a N -dimensional ball with the center at S< and radius R<:
S< = x− σ20
σ21 − σ20 δ, R< = √ 2σ40 − σ20σ21 (σ21 − σ20)2 ‖δ‖2 + 2N σ 2 0σ 2 1 σ21 − σ20 log ( σ1 σ0 ) − 2σ 2 0σ 2 1 σ21 − σ20 log(r).
Proof. From spherical symmetry of isotropic multivariate normal distribution, it follows, that without loss of generality we can take δ ≡ (a, 0, . . . , 0). With little abuse of notation, let a refer to (a, 0, . . . , 0) as well as ‖(a, 0, . . . , 0)‖. With this, the B is a set of all x, for which:
f1(x) f0(x) ≥ 1 r ⇐⇒
1
(2π)N/2σN0 exp
( − 1
2σ20 N∑ i=1 x2i
) ≤ r
(2π)N/2σN1 exp
( − 1
2σ21
[ (x1 − a)2 +
N∑ i=2 x2i
]) ⇐⇒
1
2σ21
( (x1 − a)2 +
N∑ i=2 x2i
) − 1
2σ20 N∑ i=1 x2i ≤ N log ( σ0 σ1 ) + log(r) ⇐⇒
(σ20 − σ21) N∑ i=2 x2i + (σ 2 0 − σ21)x21 − 2σ20x1a+ σ20a2 ≤ 2Nσ20σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 log(r)
Now assume σ0 > σ1 and continue:
(σ20 − σ21) N∑ i=2 x2i + (σ 2 0 − σ21)x21 − 2σ20x1a+ σ20a2 ≤ 2Nσ20σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 log(r) ⇐⇒
N∑ i=2 x2i + x 2 1 − 2σ20 σ20 − σ21 ax1 + a2σ20 σ20 − σ21 ≤ 2N σ 2 0σ 2 1 σ20 − σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 σ20 − σ21 log(r) ⇐⇒
( x1 −
σ20 σ20 − σ21 a
)2 + N∑ i=2 x2i ≤ σ20σ 2 1 (σ20 − σ21)2 a2 + 2N σ20σ 2 1 σ20 − σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 σ20 − σ21 log(r)
Such inequality defines exactly the ball from the statement of the theorem. On the other hand, if σ0 < σ1:
(σ21 − σ20) N∑ i=2 x2i + (σ 2 1 − σ20)x21 + 2σ20x1a− σ20a2 ≥ 2Nσ20σ21 log ( σ1 σ0 ) − 2σ20σ21 log(r) ⇐⇒
N∑ i=2 x2i + x 2 1 + 2σ20 σ21 − σ20 ax1 − a2σ20 σ21 − σ20 ≥ 2N σ 2 0σ 2 1 σ21 − σ20 log ( σ1 σ0 ) − 2σ 2 0σ 2 1 σ21 − σ20 log(r) ⇐⇒
( x1 +
σ20 σ21 − σ20 a
)2 + N∑ i=2 x2i ≥ 2σ40 − σ20σ21 (σ21 − σ20)2 a2 + 2N σ20σ 2 1 σ21 − σ20 log ( σ1 σ0 ) − 2σ 2 0σ 2 1 σ21 − σ20 log(r)
This is exactly the complement of a ball from the second part of the statement of the theorem.
Theorem 3:
P0(B) = χ2N
( σ20
(σ20 − σ21)2 ‖δ‖2, R2<,> σ20
) ,P1(B) = χ2N ( σ21
(σ20 − σ21)2 ‖δ‖2, R2<,> σ21
) ,
where the sign < or > is choosed according to the inequality between σ0 and σ1.
Proof. Assume first σ0 > σ1. Let us shift the coordinates, such that x + σ20
σ20−σ21 δ ←− 0. Now, the x
will have coordinates − σ 2 0 σ20−σ21 δ. Assume X ∼ N (− σ 2 0 σ20−σ21 δ, σ20I). To obtain
P0(B) = P(X ∈ B) = P ( ‖X‖2 < σ 2 0σ 2 1
(σ20 − σ21)2 ‖δ‖2 + 2N σ
2 0σ 2 1
σ20 − σ21 log ( σ0 σ1 ) + 2σ20σ 2 1 σ20 − σ21 log(r) ) ,
we could almost use NCCHSQ, but we don’t have correct scaling of variance of X . However, for any regular square matrix Q, it follows that P(X ∈ B) = P(QX ∈ QB), where QB is interpreted as set projection. Therefore, if we choose Q ≡ 1σ0 I , we will get
P0(B) = P ( ‖X/σ0‖2 <
σ21 (σ20 − σ21)2 ‖δ‖2 + 2N σ 2 1 σ20 − σ21 log ( σ0 σ1 ) + 2σ21 σ20 − σ21 log(r) ) .
Now, since X/σ0 ∼ N (− σ0σ20−σ21 δ, I), we can use the definition of NCCHSQ to obtain the final:
P0(B) = χ2N (
σ20 (σ20 − σ21)2 ‖δ‖2, σ 2 1 (σ20 − σ21)2 ‖δ‖2 + 2N σ 2 1 σ20 − σ21 log ( σ0 σ1 ) + 2σ21 σ20 − σ21 log(r) ) .
To obtain P1(B), we will do similar calculation, yet we need to compute the offset:
x+ σ20 σ20 − σ21 δ − x− δ = σ 2 1 σ20 − σ22 a.
Thus, after shifting coordinates in the same way, our alternative X will be distributed like X ∼ N (− σ 2 1
σ20−σ22 a, σ1I). Now, the same idea as before will yield the required formula.
In the case of σ0 < σ1, we do practically the same thing, yet now, we have to keep in mind, that B will not be a ball, but its complement, therefore we will obtain “1−” in the formulas.
Lemma 18. Functions ξ>(a), ξ<(a) are continuous on the whole R+. Particularly, they are continuous at 0.
Proof. Assume for simplicity σ0 > σ1 and fix x0, whose position is irrelevant and fix x1 such that ‖x0− x1‖ = am, where am is the point, where we prove the continuity. Note, that χ2N (λ, x) can be interpreted (as we have seen) as a probability of an offset ball with radius √ x and offset √ λ. Assume we have a sequence {ai}∞i=1 , ai −→ am. Define xi to be a point lying on the line defined by x0, x1 s.t. ‖x0 − xi‖ = ai. define Bai to be the worst-case ball corresponding to ai, xi. Now, without loss of generality, we can assume, that all B’s are open. We have already seen from Lemma 2, that centers of Bi converge to the center of Bm. Define Xi = 1(Bi), Xm = 1(Bm).
First we need to prove, that ri, radiuses of the balls converge. Assume for contradiction, that ri do not converge. Without loss of gererality, let rs = lim sup
i−→∞ ri > rm. If rs = ∞, it is trivial to
see, that P0(Bi) 6= pB for some i for which ri is too big. If rs < ∞, consider {aik}∞k=0 to be the subsequence for which rs is monotonically attained. Define Bs the ball with center x1 and radius rs and Xs = 1(Bs). Then, Xik
a.s.−→ Xs for k −→ ∞ and from dominated convergence theorem, P0(Bik) −→ P0(Bs). However, P0(Bs) > P0(Bm) = pB , what is contradiction, since obviously P0(Bik) 6= pB for some k. Since σ1 is fixed, the P1i , probability measures corresponding to N (xi, σ1I) are actually the same probability measure up to a shift. Therefore, P1i(Bi) can be treated as P2(B̄i), where P2 is simply measure corresponding toN (0, σ1I) and B̄i is simplyBi shifted accordingly s.t. P1i(Bi) = P2(B̄i) (and assume, without loss of generality, that for each i, they are shifted such that their centers lie on a fixed line). Now, since we know, that “position” of both the centers of the balls and the x1 is continuous w.r.t a, as can be seen from Lemma 2 (and the radiuses are still ri and converge), we see, that even 1(B̄i) −→ 1(B̄m) almost surely. Now, we can simply use dominated convergence theorem using P2 to obtain P2(B̄i) −→ P2(B̄m) and thus P1i(Bi) −→ P1m(Bm), what we wanted to prove. Note that the proof for the case σ0 < σ1 is fully analogous, yet instead of class B probabilities, we work with class A probabilities to still work with balls and not with less convenient complements.
Lemma 19. If λ1 > λ2, then χ2N (λ21, x2) ≤ χ2N (λ22, x2).
Proof. Let us fix N (0, I) and respective measure P and respective density f . From symmetry, the NCCHCSQ defined as distribution of ‖X‖2 for an offset normal distribution can be as well defined as ‖X − s‖2 under centralized normal distribution. Define B1 a ball with center at (λ1, 0, . . . , 0) and radius x and B2 a ball with center at (λ2, 0, . . . , 0) and radius x. Denote C(B) as the center
of a ball B. From definition of NCCHSQ it now follows, that P(Bi) = χ2N (λ2i , x2), i ∈ {1, 2}. Therefore, it suffices to show P(B1) ≤ P(B2). Define D1 = B1\B2 and D2 = B2\B1. Then we know:
P(B1) = ∫ B1 f(z)dz = ∫ B1∩B2 f(z)dz + ∫ B1\B2 f(z)dz = ∫ B1∩B2 f(z)dz + ∫ D1 f(z)dz
P(B2) = ∫ B2 f(z)dz = ∫ B2∩B1 f(z)dz + ∫ B2\B1 f(z)dz = ∫ B2∩B1 f(z)dz + ∫ D2 f(z)dz.
Thus,
P(B1) ≤ P(B2) ⇐⇒ ∫ D1 f(z)dz ≤ ∫ D2 f(z)dz.
Let S = C(B1)+C(B2)2 . Define a central symmetry M with center S. Let z1 ∈ D1. Then z1 can be decomposed as z1 = C(B1) + d, ‖d‖ ≤ x. Then, z2 := M(z1) = C(B2)− d from symmetry. This way, we see, that D1 = M(D2) and D2 = M(D1) under a bijection M which does not distort the geometry and distances of the euclidean space. Therefore, it suffices to show:
∀ z ∈ D2 : f(z) ≥ f(M(z)),M(z) ∈ D1.
From the monotonicity of f(y) w.r.t ‖y‖ it actually suffices to show ‖z‖ ≤ ‖M(z)‖ ∀ z ∈ D2. Fix some z ∈ D2. By the fact that M is central symmetry and z −→M(z), it is obvious, that z = S + p, M(z) = S − p, where p is some vector. Now, using law of cosine, we can write:
‖z‖2 = ‖S‖2 + ‖p‖2 − 2‖S‖‖p‖ cos(α),
where α is angle between S and −p. On the other hand:
‖M(z)‖2 = ‖S‖2 + ‖−p‖2 − 2‖S‖‖−p‖ cos(π − α).
It is obvious from these equations, that
‖z‖ ≤ ‖M(z)‖ ⇐⇒ α ≤ π/2 ⇐⇒ pTS ≤ 0.
Here, the crucial observation is, that D1 and D2 are separated by a hyperplane perpendicular to S (vector), such that S (point) is in this hyperplane. From this it follows:
y ∈ D2 =⇒ yTS ≤ ‖S‖2, y ∈ D1 =⇒ yTS ≥ ‖S‖2.
Now, since z = S+p and z ∈ D2, this implies ‖S‖2 ≥ zTS = ‖S‖2 +pTS and thus pTS ≤ 0.
Lemma 20. Functions ξ>(a), ξ<(a) are non-decreasing in a.
Proof. First assume σ0 > σ1 and analyse ξ>(a). From Lemma 18, we can without loss of generality assume, that a > 0, since ξ>(0) is simply the limit for a −→ 0 and cannot change the monotonicity status.
Now, fix a > 0 and define x0, x1 s.t. ‖x0 − x1‖ = a. Denote Ba to be the worst-case ball corresponding to x0, x1 and P1 as usual. Choose < σ 2 1
σ20−σ21 a, which is the distance between x1 and
C(Ba), as can be seen from Lemma 2.
Now, assume x2 lies on line defined by x0, x1 with ‖x0−x2‖ = a+ . LetBa+ be the corresponding worst-case ball and P2 as usual. First observe, that P2(Ba) ≥ P1(Ba), since ‖x1 − C(Ba)‖ > ‖x2 − C(Ba)‖. Here we use χ2N (λ1, x) ≤ χ2N (λ2, x) if λ1 > λ2, what is proved in Lemma 19. Second, note, that since Ba+ is the worst-case ball for x2, it follows P2(Ba+ ) ≥ P2(Ba). Thus, P1(Ba) ≤ P2(Ba+ ), but that is exactly ξ>(a) ≤ ξ>(a+ ).
To prove ξ>(a) ≤ ξ>(a + ) also for > σ 2 1
σ20−σ21 a, it suffices to consider finite sequence of points
ai starting at a and ending at a+ that are “close enough to each other” such that the respective i that codes the shift ai −→ ai+1 satisfy i < σ 2 1
σ20−σ21 ai.
Now, assume σ0 < σ1 and analysie ξ<(a). The proof is similar, but we have to be a bit careful about some details. Again, fix a > 0, define x0, x1 accordingly and all other objects as before, except now, let us denote Aa = BCa and Aa+ = B C a to be the class A balls which are complements to the anti-balls B. Again, P2(Aa) ≤ P1(Aa), since ‖x1 −C(Aa)‖ < ‖x2 −C(Aa)‖. We again used the monotonicity from Lemma 19. Moreover, P2(Aa+ ) ≤ P2(Aa), since Ba+ is the worst-case set for x2. Therefore, P1(Aa) ≥ P2(Aa+ ), but after reverting to B’s, it follows ξ>(a) ≤ ξ>(a+ ). Here, we don’t even need to care about ‖ ‖, since the centers of A’s are on the opposite half-lines from x0 than x1 and x2.
To prove the main theorem, we need a simple bound on a median of central chi-squared distribution, shown in Robert (1990) in a more general way. Lemma 21. For all c ≥ 0,
N − 1 + c ≤ χ2N,qf (c, 0.5) ≤ χ2N,qf (0.5) + c.
Proof. See Robert (1990).
Theorem 4 (the curse of dimensionality): Let x0, x1, pA, σ0, σ1, N be as usual. Then, the following two implications hold:
1. If σ0 > σ1 and
log ( σ21 σ20 ) + 1− σ 2 1 σ20 < 2 log(1− pA) N ,
then x1 is not certified w.r.t. x0. 2. If σ0 < σ1 and
log ( σ21 σ20 N − 1 N ) + 1− σ 2 1 σ20 N − 1 N < 2 log(1− pA) N ,
then x1 is not certified w.r.t. x0.
Proof. We will first prove first statement, thus let us assume σ0 > σ1. Then P1(B) = ξ>(‖x0 − x1‖). From monotonicity of ξ showed in Lemma 20, we know ξ>(‖x0 − x1‖) ≥ ξ>(0). We will show ξ>(0) > 0.5. We have, using definition of ξ> plugging in a = 0:
ξ>(0) = χ 2 N ( σ20 σ21 χ2N,qf (1− pA) ) .
Note, that here, we work with central chi-square cdf and quantile function. In order to show ξ>(0) > 0.5, it suffices to show
σ20 σ21 χ2N,qf (1− pA) ≥ N,
because it is well-known, that median of central chi-square distribution is smaller than mean, which is N , i.e. from strict monotonicity of cdf, we will get χ2N (N) > 0.5. To show the above inequality, we will use Chernoff bound on chi-squared, which states the following: If 0 < z < 1, then χ2N (zN) ≤ (z exp(1− z))N/2. Putting z ≡ σ21 σ20 , using chernoff bound we get:
χ2N ( σ21 σ20 N ) ≤ [ σ21 σ20 exp ( 1− σ 2 1 σ20 )]N 2 ! < 1− pA.
The last inequality is required to hold. If it holds, then necessarily χ2N,qf (1− pA) > σ21 σ20 N and thus σ20 σ21 χ2N,qf (1− pA) > N. Manipulating the required inequality, we will get exactly
log ( σ21 σ20 ) + 1− σ 2 1 σ20 < 2 log(1− pA) N ,
what is the assumption of 1.
Similarly we will also prove the statement 2. Assume σ0 > σ1. Like in part 1, we will just prove ξ<(0) > 0.5 by using chernoff bound. This time, however, we have:
ξ<(0) = 1− χ2N ( σ20 σ21 χ2N,qf (pA) ) ,
i.e. we need to prove
χ2N ( σ20 σ21 χ2N,qf (pA) ) < 1 2 .
The second part of Chernoff bound states: If 1 < z, then χ2N (zN) ≥ 1− (z exp(1− z))N/2. Let us choose z ≡ σ 2 1
σ20
N−1 N . Then Chernoff bound yields:
χ2N ( σ21 σ20 N − 1 N N ) ≥ 1− [ σ21 σ20 N − 1 N exp ( 1− σ 2 1 σ20 N − 1 N )] ! > pA.
If this holds, then
σ21 σ20 N − 1 N N > χ2N,qf (pA) ⇐⇒ σ20 σ21 χ2N,qf (pA) < N − 1.
Now, using Lemma 21 for the easy case of central chi-squared, we see: χ2N (N − 1) < 0.5 and thus
χ2N ( σ20 σ21 χ2N,qf (pA) ) < 1 2 ,
what we wanted to prove.
Corollary 5 (one-sided simpler bound): Let x0, x1, pA, σ0, σ1, N be as usual and assume now σ0 > σ1. Then, if
σ1 σ0 <
√ 1− 2 √ − log(1− pA)
N ,
then x1 is not certified w.r.t x0.
Proof. We will simply prove
σ1 σ0 <
√ 1− 2 √ − log(1− pA)
N =⇒ log ( σ21 σ20 ) + 1− σ 2 1 σ20 < 2 log(1− pA) N .
Assume expression log(1−y)+y; 1 > y > 0. From Taylor series, it is apparent, that log(1−y)+y < −y 2
2 . Therefore, if − y2 2 < 2 log(1−pA) N , then also log(1− y) + y < 2 log(1−pA) N . Solving for y in the first inequality, we get sufficient condition y > 2 √ − log(1−pA)
N . Plugging 1− σ21 σ20 into y we get:
1− σ 2 1
σ20 > 2
√ − log(1− pA)
N ,
which is very easily manipulated to the inequality from theorem statement.
Theorem 6: Let x0, x1, pA, σ0 be as usual and let ‖x0 − x1‖ = R. Then, the following two statements hold:
1. Let σ1 ≤ σ0. Then, for all σ2 : σ1 ≤ σ2 ≤ σ0, if ξ>(R, σ2) > 0.5, then ξ>(R, σ1) > 0.5. 2. Let σ1 ≥ σ0. Then, for all σ2 : σ1 ≥ σ2 ≥ σ0, if ξ<(R, σ2) > 0.5, then ξ>(R, σ1) > 0.5.
Proof. We will first prove the first statement. Denote, as usual in the proofs Bi the worst-case ball for σi, Pi the probability associated toN (x1, σ2i I). Since ξ>(R, σ2) > 0.5 and since it is essentially P2(B2), we see, that the probability of a ball under normal distribution is bigger than half. This is obviously possible just if x1 ∈ B2. From the fact, that B1 is the worst-case ball for σ1 we see ξ>(R, σ1) = P1(B1) ≥ P1(B2). It suffices to show P1(B2) ≥ P2(B2). This follows, since σ1 ≤ σ2 and x1 ∈ B2. We know, that we can rescale the space such that
P1(B2) = P2 ( σ2 σ1 (B2 − x1) + x1 ) ,
using the fact that σ just scales the normal distribution. The set σ2σ1 (B2 − x1) + x1 is just an image of B2 via homothety with center x1 and rate σ2σ1 . So it suffices to prove
P2 ( σ2 σ1 (B2 − x1) + x1 ) ≥ P2(B2).
However, obviously σ2σ1 (B2 − x1) + x1 ⊃ B2 from convexity of a ball. If, namely, x1 + z ∈ B2, then from convexity also x1 + σ1σ2 z ∈ B2 and this maps back to x1 + z, thus x1 + z is in an image. Applying monotonicity of P, we obtain the result.
Now we will prove the second statement and as usual, let A1, A2,P1,P2 be as usual (A is now the ball connected to class A). As always, P1(A1) ≤ P1(A2), so it suffices to show P2(A2) < 0.5 =⇒ P1(A2) < 0.5. Now, we need to distinguish two cases. If x1 ∈ A2, the proof is completely analogical to the first part, but now reasoning on A’s rather than B’s. In this case, we will even get stronger P1(A2) ≤ P2(A2) just like in the first part. If x1 6∈ A2, then it is easy to see, that indeed P2(A2) < 0.5, yet it is also obvious to see that P1(A2) < 0.5. This finishes the proof of the lemma.
Theorem 7: Let σ(x) be r-semi-elastic function and x0, pA, N, σ0 as usual. Then, the certified radius at x0 guaranteed by our method is
CR(x0) = max {0, sup {R ≥ 0; ξ>(R, σ0 exp(−rR)) < 0.5 and ξ<(R, σ0 exp(rR)) < 0.5}} .
Proof. This follows easily from Theorem 6 Lemma 22. Let us have f(x) = M∑ i=1 fi(x)1(x ∈ Ri), where {Ri}Mi=1 is finite set of regions that divide the RN and {fi}Mi=1, fi : Ri −→ R is finite set of 1-Lipschitz continuous (1-LC) functions. Moreover assume, that f(x) is continuous. Then, f(x) is 1-LC.
Proof. We do not assume any nice behaviour from our decision regions, what can make the situation quite ugly. For instance, regions might not be measurable. However, it will not be a problem for us.
Fix x1, x2. Consider line segment S = x1 + α(x2 − x1), α ∈ [0, 1]. Let us instead of points in S work with numbers in [0, 1] via the α encoding. Consider the following coloring C of [0, 1]: Each point a ∈ [0, 1] will be assigned one of M colors according to which region the x1 + a(x2 − x1) belongs. Define d1 ≡ 0 and d2 = sup{z ∈ [0, 1], C(z) = C(d1)}. Thus, d2 is the supremum of all numbers colored the same color as 0. Then,
|f(x1 + d2(x2 − x1))− f(x1 + d1(x2 − x1))| ≤ (d2 − d1)‖x2 − x1‖.
Why? Let {zj}∞j=1 be a non-decreasing sequence s.t. C(zj) = C(0) and zj −→ d2. Since f is continuous, obviously f(x1 + d2(x2 − x1)) = lim
j−→∞f(x1 + zj(x2 − x1)). Now, since norm and absolute value are both continuous functions and since from 1-LC of fC(0) on RC(0) we have
∀j ∈ N : |f(x1 + zj(x2 − x1))− f(x1 + d1(x2 − x1))| (zj − d1)‖x2 − x1‖ ≤ 1,
we also necessarily have
|f(x1 + d2(x2 − x1))− f(x1 + d1(x2 − x1))| (d2 − d1)‖x2 − x1‖ ≤ 1.
If d2 = 1, we finish the construction. If not, distinguish two cases. First assume C(d2) = C(d1). In this case, take some color C s.t. ∃{zj}∞j=1 : zj+1 ≤ zj ∀j ∈ N and zj −→ d2 and zj > d2 ∀j ∈ N, and fix one such {zj}∞j=1. Obviously, C 6= C(0), since d2 is upper-bound on points of color C(0). Then, define d3 = sup{z ∈ [0, 1], C(z) = C} and also define {zj}∞j=1 : zj+1 ≥ zj ∀j ∈ N and zj −→ d3. From continuity of f , we again have f(x1 +d2(x2−x1)) = lim
j−→∞f(x1 +zj(x2−
x1)) and similarly f(x1 + d3(x2 − x1)) = lim j−→∞f(x1 + zj(x2 − x1)). Again from continuity of absolute value and norm and 1-LC of all the partial functions we have:
∀j ∈ N : |f(x1 + zj(x2 − x1))− f(x1 + zj(x2 − x1))| (zj − zj)‖x2 − x1‖ ≤ 1
and |f(x1 + d3(x2 − x1))− f(x1 + d2(x2 − x1))|
(d3 − d2)‖x2 − x1‖ ≤ 1.
Now assume C(d2) 6= C(d1). Then, we can take as C directly C(d2) and do the same as in the last paragraph (note, that this case could have implicitly come up in the previous construction too, but we would need to not take C = C(0) and we find this case distinction to be more elegant).
If d3 = 1, we finish the construction. If not, we continue in exactly the same manner as before. Since the number of colors M is finite, we will run out of colors in finite number of steps and thus, eventually there will be l ≤M s.t. dl = 1. The final 1-LC is now trivially obtained as follows:
|f(x1 + 1(x2 − x1))− f(x1 + 0(x2 − x1))| = ∣∣∣∣ l−1∑ i=1 f(x1 + di+1(x2 − x1))− f(x1 + di(x2 − x1)) ∣∣∣∣
≤ l−1∑ i=1 ∣∣∣f(x1 + di+1(x2 − x1))− f(x1 + di(x2 − x1))∣∣∣ ≤ l−1∑ i=1 ∣∣(di+1 − di)‖x2 − x1‖∣∣ = ‖x2 − x1‖
l−1∑ i=1 (di+1 − di) = ‖x2 − x1‖
Theorem 8: The σ(x) defined in Equation 1 is r-semi-elastic.
Proof. Our aim is to prove, that
log(σ(x)) = log(σb) + r 1 k  ∑ xi∈Nk(x) ‖x− xi‖ −m 
is r lipschitz continuous. Obviously, this does not depend neither on log(σb), nor on −rm, so we will focus just on rk ∑ xi∈Nk(x) ‖x− xi‖. Obviously, this function is r lipschitz continuous if and only
if 1k ∑
xi∈Nk(x) ‖x− xi‖ is 1-LC.
Let us fix y ∈ RN . We will first prove ‖x− y‖ is 1-LC. Let us fix x1, x2. From triangle inequality we have ∣∣‖x1 − y‖ − ‖x2 − y‖∣∣ ≤ ‖x1 − x2‖, what is exactly what we wanted to prove.
Now fix y1, y2, . . . , yk and x1, x2. Then∣∣∣∣1k k∑ i=1 ‖x1 − yi‖ − 1 k k∑ i=1 ‖x2 − yi‖ ∣∣∣∣ = 1k ∣∣∣∣ k∑ i=1 ‖x1 − yi‖ − ‖x2 − yi‖ ∣∣∣∣
≤ 1 k k∑ i=1 ∣∣∣‖x1 − yi‖ − ‖x2 − yi‖∣∣∣ ≤ 1 k k∑ i=1 ‖x1 − x2‖ = 1
Finally note, that using the k nearest neighbors out of finite training dataset will divide RN in a finite number of regions, where each region is defined by the set of k nearest neighbors for x in that region. Note, that the average distance from k nearest neighbors is obviously continuous. Then, using Lemma 22, the claim follows.

Review summaries:
The authors study input-dependent adversarial smoothing i.e. the setting where the smoothing distribution has noise level $\sigma$ dependent on $x$. They begin by deriving the robustness certificate in this setting, and show that it suffers from the curse of dimensionality -- in order to maintain non-vacuous certificates in high dimensions, $\sigma(x)$ must be prohibitively smooth. They then provide empirical experiments using an example function $\sigma(x)$ on CIFAR-10 and MNIST and show that performance is more or less the same compared to the fixed noise level baseline (Cohen et al. 2019).

This paper studies the input-dependent randomized smoothing technique - it adds isotropic Gaussian noise with different variance on different input x to give robustness certificates.  The hope is that this can perform better than applying a uniform noise.  The main theoretical claim in this paper is that when adding input dependent isotropic Gaussian noise, the variance of the noise should not change too much among different inputs. Based on this, the authors propose an input dependent smoothing mechanism based on k-nearest neighbor and evaluate its performance on MNIST and CIFAR datasets.

This work extends the theoretical results from Cohen et.al to the input dependent setup. In particular, the derive the theoretical conditions for the validity of data-dependent randomized smoothing through a similar procedure to the one in Cohen et.al. At last, a function that optimizes the smoothing parameters per input was proposed and validated experimentally on MNIST and CIFAR10.

The authors propose an input-dependent randomized smoothing method with non-constant input-dependent $\sigma(x)$.  ### contributions  - The authors generalize randomized smoothing with non-constant input-dependent $\sigma(x)$. - The proposed method successfully addresses some intrinsic limitations/flaws in the previous (input-independent) RS.