Double Sampling Randomized Smoothing

Linyi Liuiuc
Jiawei Zhanguiuc
Tao Xiepku
Bo Liuiuc

uiucUniversity of Illinois Urbana-Champaign, Illinois, USA
pkuPeking University, Beijing, China

Linyi<EMAIL_ADDRESS>Tao<EMAIL_ADDRESS>Bo<EMAIL_ADDRESS>
Certified Deep Learning, Robustness, Verification

Neural networks (NNs) are known to be vulnerable against adversarial perturbations, and thus there is a line of work aiming to provide robustness certification for NNs, such as randomized smoothing, which  samples smoothing noises from a certain distribution to certify the robustness for a smoothed classifier.
However, as shown by previous work, the certified robust radius in randomized smoothing suffers from scaling to large datasets (“curse of dimensionality”).
To overcome this hurdle, we propose a  () framework, which exploits the sampled probability from an additional smoothing distribution to tighten the robustness certification of the previous smoothed classifier.
under mild assumptions, we prove that can certify $\Theta(\sqrt{d})$ robust radius under $\ell_2$ norm where $d$ is the input dimension, implying  that
may be able to break the
curse of dimensionality of randomized smoothing.
We instantiate for a generalized family of Gaussian smoothing and propose an efficient and sound computing method based on customized dual optimization considering sampling error.
Extensive experiments on MNIST, CIFAR-10, and ImageNet verify our theory and show that certifies larger robust radii than existing baselines consistently under different settings.
Code is available at <https://github.com/llylly/DSRS>.

§ INTRODUCTION

Neural networks (NNs) have achieved great advances on a wide range of tasks, but
have been shown
vulnerable against adversarial examples [Szegedy et al., 2014, Goodfellow et al., 2015, Eykholt et al., 2018, Wang et al., 2021, Qiu et al., 2020, Li et al., 2020, Zhang et al., 2021].
A plethora of empirical defenses are proposed to improve the robustness; however, most of these are broken by  strong adversaries again [Carlini & Wagner, 2017, Athalye et al., 2018, Tramèr et al., 2020].
Recently, there are great efforts in developing certified defenses for NNs under certain adversarial constraints [Wong & Kolter, 2018, Raghunathan et al., 2018, Li et al., 2020, Xu et al., 2020, Li et al., 2019].

Upper: Standard certification for randomized smoothing leverages information from only one distribution ($P_A$) to compute robustness certification.
leverages information from two distributions ($P_A$ and $Q_A$) to compute certification for the smoothed classifier, yielding significantly larger certified radius.

Randomized smoothing [Cohen et al., 2019, Li et al., 2019] has emerged as a popular technique to provide certified robustness for large-scale datasets.
Concretely, it samples noise from a certain smoothing distribution to construct a smoothed classifier, and thus certifies the robust radius for the smoothed classifier.
Compared to other techniques [Wong & Kolter, 2018, Mirman et al., 2018, Gowal et al., 2019, Zhang et al., 2020],
randomized smoothing is efficient and agnostic to the model, and is applicable to a wide range of ML models, including large ResNet [He et al., 2016] on the ImageNet dataset.

To improve the certified robust radius, existing studies [Cohen et al., 2019, Lee et al., 2019, Li et al., 2019, Yang et al., 2020] have explored different smoothing distributions.
However, the improvement is limited.
For example, $\ell_2$ certified robust radius does not increase on large datasets despite that the input dimension $d$ increases [Cohen et al., 2019], resulting in a low $\ell_\infty$ certified radius on large datasets, theoretically shown as an intrinsic barrier of randomized smoothing (“curse of dimensionality” or “$\ell_\infty$ barrier”) [Yang et al., 2020, Blum et al., 2020, Kumar et al., 2020, Hayes, 2020, Wu et al., 2021].

Given these challenges toward tight robustness certification, a natural question arises:
Q1: Is it possible to circumvent the barrier of randomized smoothing by certifying with additional “information"?
Q2: What type of information is needed to provide tight robustness certification?
To answer these questions, we propose a   () framework to leverage the sampled noises from an additional smoothing distribution as additional information to tighten the robust certification.
In theory, we show that
(1) ideally, if the decision region of the base classifier is known, can provide tight robustness certification;
(2) more practically, if the inputs, which can be correctly classified by the base classifier, satisfy the concentration property within an input-centered ball with constant mass under standard Gaussian measure, the standard Neyman-Pearson-based certification [Li et al., 2019, Cohen et al., 2019, Salman et al., 2019, Yang et al., 2020] can  certify only a dimension-independent $\ell_2$ radius, whereas with generalized Gaussian smoothing can certify radius $\Omega(\sqrt d)$ (under $\ell_2$ norm),  which would increase with the dimension $d$, leading to tighter certification.
Under more general conditions, we provide numerical simulations to verify our theory.
Our results provide a positive answer to Q1 and sufficient conditions for Q2, i.e., may be able to circumvent the
barrier of randomized smoothing.

Motivated by the theory, we leverage a type of generalized Gaussian [Zhang et al., 2020] as the smoothing distribution and truncated generalized Gaussian as an additional distribution.
For this type of concretization, we propose an efficient and sound computation method to compute the certifiably robust radius for practical classifiers considering sampling error.
Our method formulates the certification problem given additional information as a constrained optimization problem and  leverages specific properties of the dual problem to decompose the effects of different dual variables to solve it.
is fully scalable since the computational time is nearly independent of the size of the dataset, model, or  sampling.
Our extensive experimental evaluation on MNIST, CIFAR-10, and ImageNet shows that
(1) under large sampling size ($2\times 10^5 - 8\times 10^5$),
the certified radius of consistently increases as suggested by our theory;
(2) under practical sampling size ($10^5$), DSRS can certify consistently higher robust radii than existing baselines, including standard Neyman-Pearson-based certification.

As further discussed in <Ref>, we believe that as a framework opens a wide range of future directions for selecting or optimizing different forms of additional information to tighten the certification of randomized smoothing.

We summarize the main technical contributions as follows:

* We propose a general robustness certification framework ,  which leverages additional information by sampling from another smoothing distribution.

* We prove that under practical concentration assumptions, certifies $\Omega(\sqrt d)$ radius under $\ell_2$ norm with $d$ the input dimension, suggesting a possible way to circumvent the intrinsic barrier of randomized smoothing.

* We concretize by generalized Gaussian smoothing mechanisms and propose a method to efficiently compute the certified radius for given classifiers.

* We conduct extensive experiments, showing that provides consistently tighter robustness certification than existing baselines, including standard Neyman-Pearson-based certification across different models on MNIST, CIFAR-10, and ImageNet.

Related Work.
For the certification method of randomized smoothing, most existing methods leverage only the true-class prediction probability to certify.
In this case, the tightest possible robustness certification is based on the Neyman-Pearson lemma [Neyman & Pearson, 1933] as first proposed by Cohen et al., 2019 for certifying $\ell_2$ radius under Gaussian smoothing.
Several methods extend this certification to accommodate different smoothing distributions and different $\ell_p$ norms [Dvijotham et al., 2020, Yang et al., 2020, Zhang et al., 2020, Levine & Feizi, 2021].
In randomized smoothing,
the $\ell_2$ certified robust radius $r$ is similar across datasets of different scales, resulting in the vanishing $\ell_\infty$ certified radius $r/{\sqrt d}$ when input dimension increases.
This limitation of existing certification methods of randomized smoothing is formally proved [Yang et al., 2020, Blum et al., 2020, Kumar et al., 2020, Hayes, 2020, Wu et al., 2021] and named “$\ell_\infty$ barrier” or “curse of dimensionality”.

Recent work tries to incorporate additional information besides true-class prediction probability to tighten the certification and bypass the barrier.
For $\ell_2$ and $\ell_\infty$ certification, to the best of our knowledge, gradient magnitude is the only exploited additional information [Mohapatra et al., 2020, Levine et al., 2020].
However, in practice, the improvement is relatively marginal and requires a large number of samples (see <Ref>).
Some other methods provide tighter certification given specific model structures [Kumar et al., 2020, Chiang et al., 2020, Lee et al., 2019, Awasthi et al., 2020].
instead focuses on leveraging model-structure-agnostic additional information.
More discussion on related work is in <Ref>.

§ PRELIMINARIES AND BACKGROUND

Define $[C] := \{1,\,\dots\,C\}$.
Let $\vDelta^C$ be the $C$-dimensional probability simplex.
We consider a multiclass classification model $F: \sR^d \to [C]$ as the base classifier, where $d$ is the input dimension, and the model outputs hard-label class prediction within $[C]$.
The original smoothing distribution $\gP$ and additional smoothing distribution $\gQ$ are both supported on $\sR^d$.
We let $p(\cdot)$ and $q(\cdot)$ be their density functions respectively.
We assume that both $p$ and $q$ are positive and differentiable almost everywhere, i.e., the set of singular points has zero measure under either $\gP$ or $\gQ$.
These assumptions hold for common smoothing distributions used in the literature such as Gaussian distribution [Lécuyer et al., 2019, Li et al., 2019, Cohen et al., 2019, Yang et al., 2020].

Randomized smoothing constructs a smoothed classifier from a given base classifier by adding input noise following original smoothing distribution $\gP$.
For input $\vx \in \sR^d$,
we define prediction probability under $\gP$ by function $f^\gP: \sR^d \to \vDelta^C$:

\begin{equation}
% \small
f^\gP(\vx)_c := \Pr_{\vepsilon\sim\gP} [F(\vx + \vepsilon) = c] \quad \text{where} \quad c\in [C].
\label{eq:f}
\vspace{-0.5em}
\end{equation}

The smoothed classifier $\wF^\gP: \sR \to [C]$ (or $\wF$ when $\gP$ is clear from the context) predicts the class with the highest confidence after smoothing with $\gP$:

\begin{equation}
% \small
\wF^\gP(\vx) := \argmax_{c\in [C]} f^\gP(\vx)_c.
\label{eq:F}
\end{equation}

We focus on robustness certification against $\ell_p$-bounded perturbations for smoothed classifier $\wF$, where the standard certification method is called Neyman-Pearson-based certification [Cohen et al., 2019] (details in <Ref>, certified radius from it denoted by $\rNP$).
Concretely, certification methods compute robust radius $r$ defined as below.

Under $\ell_p$ norm ($p \in \sR_+ \cup \{+\infty\}$), for given smoothed classifier $\wF^\gP$ and input $\vx_0 \in \sR^d$ with true label $y_0\in [C]$,
a radius $r \ge 0$ is called certified (robust) radius for $\wF^\gP$ if $\wF^\gP$ always predicts $y_0$ for any input within the $r$-radius ball centered at $\vx_0$:

\begin{equation}
\forall \vdelta \in \sR^d, \|\vdelta\|_p < r,
\wF^\gP(\vx_0 + \vdelta) = y_0.
\end{equation}

§ OVERVIEW

We propose (), which leverages the prediction probability from an additional smoothing distribution $\gQ$ (formally $Q_A := f^\gQ(\vx_0)_{y_0} = \Pr_{\vepsilon\sim\gQ}[F(\vx+\vepsilon)=y_0]$), along with the prediction probability from the original smoothing distribution $\gP$ (formally $P_A := f^\gP(\vx_0)_{y_0}$ as in <Ref>, also used in Neyman-Pearson-based certification), to provide robustness certification for $\gP$-smoothed classifier $\wF^\gP$.
Note that both $P_A$ and $Q_A$ can be obtained from Monte-Carlo sampling (see <Ref>).
Formally, we let $\rDSRS$ denote the tightest possible certified radius with prediction probability from $\gQ$, then $\rDSRS$ can be defined as below.

Given $P_A$ and $Q_A$,

\begin{equation}
% \small
\label{eq:r-dsrs}
\begin{aligned}
& \rDSRS := \max r \quad \mathrm{s.t.} \\
& \forall F: \sR \to [C], f^\gP(\vx_0)_{y_0} = P_A, f^\gQ(\vx_0)_{y_0} = Q_A \\
\Rightarrow & \forall \vx, \|\vx - \vx_0\|_p < r, \wF^\gP(\vx) = y_0.
\end{aligned}
\end{equation}

Intuitively, $\rDSRS$ is the maximum possible radius, such that any smoothed classifier constructed from base classifier satisfying $P_A$ and $Q_A$ constraints cannot predict other labels when the perturbation magnitude is within the radius.

In <Ref>, we will analyze the theoretical properties of , including comparing $\rDSRS$ and $\rNP$ under the concentration assumption.
Computing $\rDSRS$ is nontrivial, so in <Ref>, we will introduce a practical computational method that exactly solves $\rDSRS$ when $\gP$ and $\gQ$ are standard and generalized (truncated) Gaussian.
In <Ref>, we will show method variants to deal with other forms of $\gP$ and $\gQ$ distributions.
In <Ref>, we will further generalize the framework.

Smoothing Distributions.
Now we formally define the smoothing distributions used in .
We mainly consider standard Gaussian $\gN$ [Cohen et al., 2019, Yang et al., 2020] and generalized Gaussian $\gNg$ [Zhang et al., 2020].
Let $\gN(\sigma)$ to represent standard Gaussian distribution with covariance matrix $\sigma^2\mI_d$ that has density function $\propto \exp(-\|\vepsilon\|_2^2/(2\sigma^2))$.[In this paper, $\gN(\sigma)$ is a shorthand of $\gN(0,\sigma^2\mI_d)$.]
For $k\in\sN$, we let $\gNg(k,\sigma)$ to represent generalized Gaussian whose density function $\propto \|\vepsilon\|_2^{-2k} \exp(-\|\vepsilon\|_2^2/(2\sigma'^2))$ where $\sigma' = \sqrt{d/(d-2k)}\sigma$.
Here we use $\sigma'$ instead of $\sigma$ to ensure that the expected noise $\sqrt{\E \|\vepsilon\|_2^2}$ of $\gNg(k,\sigma)$ is the same as $\gN(\sigma)$.
The generalized Gaussian as the smoothing distribution overcomes the “thin shell” problem of standard Gaussian and improves certified robustness [Zhang et al., 2020]; and we will reveal more of its theoretical advantages in <Ref>.

As the additional smoothing distribution $\gQ$, we will mainly consider truncated distributions within a small $\ell_2$ radius ball.
Specially, truncated standard Gaussian is denoted by $\gNtrunc(T,\sigma)$ with density function $\propto \exp(-\|\vepsilon\|_2^2/(2\sigma^2)) \cdot \1[\|\vepsilon\|_2 \le T]$;
and truncated generalized Gaussian is denoted by $\gNgtrunc(k,T,\sigma)$ with density function $\propto \|\vepsilon\|_2^{-2k}\exp(-\|\vepsilon\|_2^2/(2\sigma'^2)) \cdot \1[\|\vepsilon\|_2 \le T]$.

In <Ref>, we summarize these distribution definitions.

Definitions of smoothing distributions in this paper. In the table, $k\in \sN$, $\sigma' = \sqrt{d/(d-2k)}\sigma$.

Name     Notation     Density Function

Standard Gaussian     $\gN(\sigma)$     $\propto \exp\left(-\frac{\|\vepsilon\|_2^2}{2\sigma^2}\right)$

Generalized Gaussian     $\gNg(k,\sigma)$     $\propto \|\vepsilon\|_2^{-2k} \exp\left(-\frac{\|\vepsilon\|_2^2}{2\sigma'^2}\right)$

Truncated Standard Gaussian     $\gNtrunc(T,\sigma)$     $\propto \exp\left(-\frac{\|\vepsilon\|_2^2}{2\sigma^2}\right)\cdot \1[\|\vepsilon\|_2\le T]$

Truncated Generalized Gaussian     $\gNgtrunc(k,T,\sigma)$     $\propto \|\vepsilon\|_2^{-2k} \exp\left(-\frac{\|\vepsilon\|_2^2}{2\sigma'^2}\right)\cdot \1[\|\vepsilon\|_2\le T]$

§ THEORETICAL ANALYSIS OF

In this section, we theoretically analyze to answer the following core question:
Does $Q_A$, the prediction probability under additional smoothing distribution, provide sufficient information for tightening the robustness certification?
We first show that if the support of $\gQ$ is the decision region of true class, can certify the smoothed classifier's maximum possible robust radius.
Then, under concentration assumption, we show the $\ell_2$ certified radius of can be $\Omega(\sqrt d)$ that is asymptotically optimal for bounded inputs.
Finally, under more general conditions, we conduct both numerical simulations and real-data experiments to verify that the certified radius of increases with data dimension $d$.
These analyses provide a positive answer to the above core question.

§.§ can certify the tightest possible robust radius.

Given an original smoothing distribution $\gP$ and a base classifier $F_0$.
At input point $\vx_0 \in \sR^d$ with true label $y_0$, we define the tightest possible certified robust radius $\rtight$ to be the largest $\ell_p$ ball that contains no adversarial example for smoothed classifier $\wF_0^\gP$:

\begin{equation*}
\rtight := \max r \, \mathrm{s.t.} \, \forall \vdelta\in\sR^d, \|\vdelta\|_p < r, \wF_0^\gP(\vx_0+\vdelta) = y_0.
\vspace{-2mm}
\end{equation*}

Then, for binary classification, if we choose an additional smoothing distribution $\gQ$ whose support is the decision region or its complement, then can certify robust radius $\rtight$.

Suppose the original smoothing distribution $\gP$ has non-zero density everywhere, i.e., $p(\cdot) > 0$.
For binary classification with base classifier $F_0$,
at point $\vx_0\in\sR^d$,
let $\gQ$ be an additional distribution that satisfies: (1) its support is the decision region of an arbitrary class $c\in [C]$ shifted by $\vx_0$: $\supp(\gQ) = \{\vx - \vx_0: F_0(\vx) = c\}$;
(2) for any $\vx \in \supp(\gQ)$, $0 < q(\vx)/p(\vx) < +\infty$.
Then, plugging $P_A = f_0^\gP(\vx_0)_c$ and $Q_A = f_0^\gQ(\vx_0)_c$ (see <Ref>) into <Ref>, we have $\rDSRS = \rtight$ under any $\ell_p$ ($p \ge 1$).

We defer the proof to <Ref>.
At a high level, with this type of $\gQ$, we have $Q_A = 1$ or $Q_A = 0$.
Then, from the mass of the $\gQ$'s support on $\gP$ and $P_A$, we can conclude that the $\gQ$'s support is exactly the decision region of label $c$ or its complement.
Thus, the DSRS constraints (in <Ref>) are satisfied iff $F$ differs from $F_0$ in a zero-measure set, and thus we exactly compute the smoothed classifier $\wF_0^\gP$'s maximum certified robust radius in .
An extension to multiclass setting is in <Ref>.

For any base classifier $F_0$, $\gQ$ that satisfies conditions in <Ref> exists, implying that with , certifying a strictly tight robust radius is possible.
In contrast, Neyman-Pearson-based is proved to certify tight robust radius for linear base classifiers <cit.>, but for arbitrary base classifiers, its tightness is not guaranteed.
This result suggests that, to certify a tight radius, just one additional smoothing distribution $\gQ$ is sufficient rather than multiple ones.

On the other hand, it is challenging to find $\gQ$ whose support (or its complement) exactly matches the decision region of an NN classifier.
In the following, we analyze the tightness of under weaker assumptions.

§.§ can certify $\Omega(\sqrt d)$ $\ell_2$ radius under concentration assumption.

We begin by defining the concentration property.

Given a base classifier $F_0$, at input $\vx_0 \in \sR$ with true label $y_0$, we call $F_0$ satisfies ($\sigma,\Pcon$)-concentration property, if for within $\Pcon$-percentile of small $\ell_2$ magnitude Gaussian $\gN(\sigma)$ noise, the adversarial example occupies zero measure.
Formally, ($\sigma,\Pcon$)-concentration means

\begin{align}
& \Pr_{\vepsilon \sim \gN(\sigma)} [F_0(\vx_0 + \vepsilon) = y_0 \,|\, \|\vepsilon_0\|_2 \le T] = 1
\label{eq:concentration}  \\
\mathrm{where} & \quad T \, \mathrm{satisfies} \, \Pr_{\vepsilon \sim \gN(\sigma)} [\|\vepsilon\|_2 \le T] = \Pcon.
\end{align}

Intuitively, ($\sigma,\Pcon$)-concentration implies that the base classifier has few adversarial examples for small magnitude noises during standard Gaussian smoothing.
In <Ref> (in <Ref>), we empirically verified that a well-trained base classifier on ImageNet may satisfy this property for a significant portion of inputs.
Furthermore, Salman et al., 2019 show that promoting this concentration property by adversarially training the smoothed classifier improves the certified robustness.
With this concentration property, certifies the radius $\Omega(\sqrt{d})$ under $\ell_2$ norm, as the following theorem shows.

Let $d$ be the input dimension and $F_0$ be the base classifier.
For an input point $\vx_0 \in \sR^d$ with true class $y_0$, suppose $F_0$ satisfies $(\sigma,\Pcon)$-Concentration property.
Then, for any sufficiently large $d$,
for the classifier $\wF_0^{\gP'}$ smoothed by generalized Gaussian $\gP'=\gNg(k,\sigma)$ with $d/2 - 15 \le k < d/2$, with additional smoothing distribution $\gQ = \gNgtrunc(k,T,\sigma)$ can certified $\ell_2$ radius

\begin{equation}
% \small
\rDSRS \ge 0.02\sigma\sqrt{d}
\vspace{-0.5em}
\end{equation}

where $T = \sigma\sqrt{2\GammaCDF_{d/2}^{-1}(\Pcon)}$ and $\GammaCDF_{d/2}$ is the CDF of gamma distribution $\Gamma(d/2,1)$.

We defer the proof to <Ref>.
At high level,
based on the standard Gaussian distribution's property (<Ref>), we find $Q_A = 1$ under concentration property (<Ref>).
With $Q_A = 1$, we derive a lower bound of $\rDSRS$ in <Ref>.
We then use: (1) the concentration of beta distribution $\Beta(\frac{d-1}{2},\frac{d-1}{2})$ (see <Ref>) for large $d$;
(2) the relative concentration of gamma $\Gamma(d/2,1)$ distribution around mean for large $d$ (see <Ref> and resulting <Ref>);
and (3) the misalignment of gamma distribution $\Gamma(d/2-k,1)$'s mean and median for small $(d/2-k)$ (see <Ref>) to lower bound the quantity in <Ref> and show it is large or equal to $0.5$.
Then, using the conclusion in <Ref> we conclude that $\rDSRS \ge 0.02\sigma\sqrt{d}$.

(1) For standard Neyman-Pearson based certificaton, $\rNP = \sigma\Phi^{-1}(f^\gP(\vx_0)_{y_0})$.
Along with the increase of input dimension $d$, to achieve growing $\ell_2$ certified radius, one needs the prediction probability of true class under $\gP$, namely $f^\gP(\vx_0)_{y_0}$, to grow simultaneously, which is challenging.
Indeed, across different datasets, $f^\gP(\vx_0)_{y_0}$ is almost a constant, which leads to a constant $\ell_2$ certified radius and shrinking $\ell_\infty$ radius for large $d$.
We further empirically illustrate this property in <Ref>.

(2) In contrast, as long as the model satisfies concentration property, which may be almost true on large datasets as reflected by <Ref>,
with our specific choices of $\gP$ and $\gQ$,
can achieve $\Omega(\sigma\sqrt d)$ $\ell_2$ radius on large datasets.
This rate translates to a constant $\Omega(\sigma)$ $\ell_\infty$ radius on large datasets and thus breaks the curse of dimensionality of randomized smoothing.
We remark that this $\sqrt d$ rate is optimal when dataset input is bounded such as images (otherwise, the $\omega(1)$ $\ell_\infty$ radius leads the radius to exceed the constant $\ell_\infty$ diameter for large $d$).
Therefore, under the assumption of concentration property, provides asymptotically optimal certification for randomized smoothing.

(3) Smoothing with generalized Gaussian distribution and choosing a parameter $k$ that is close to $d/2$ play an essential role in proving the $\Omega(\sigma\sqrt d)$ certified radius. Otherwise, in <Ref> we have <Ref> that shows any certification methods cannot certify an $\ell_2$ radius $c\sqrt d$ for any $c > 0$.
This adds another theoretical evidence for the superiority of generalized Gaussian that is cross-validated by Zhang et al., 2020.

[When holding probability in <Ref> is obtained from sampling $N$ times and confidence level $1-\alpha=99.9\%$, relation between certified radius ($y$-axis) and input data dimension $d$ ($x$-axis).
Different curves correspond to different $N$s.]

[Relation between certified radius ($x$-axis) and certified accuracy ($y$-axis) on ImageNet models.
Different curves correspond to Neyman-Pearson and certification with different $N$s. Sampling error considered, confidence level $=99.9\%$.]

Tendency of certified robust radius considering sampling error.
In both (a) and (b), certified radius grows along with the increase of sampling number $N$ but Neyman-Pearson radius is almost fixed.

§.§ certifies tighter radius under general scenarios.

When the concentration property does not absolutely hold, a rigorous theoretical analysis becomes challenging, since the impact of a noninfinite dual variable needs to be taken into account.
This dual variable is inside a Lambert $W$ function where typical approximation bounds are too loose to provide non-trivial convergence rates.
Thus, we leverage the numerical computational method introduced in <Ref> to provide numerical simulations and real-data experiments.
We generalize the concentration assumption by changing the holding probability in <Ref> from $1$ to $\alpha^{1/N}$, which corresponds to $(1-\alpha)$-confident lower bound of $Q_A$ given $N$ times of Monte-Carlo sampling, where we set $\alpha=0.1\%$ following the convention [Cohen et al., 2019].
In this scenario, we compare certification with Neyman-Pearson certification numerically in <Ref> (numerical simulations in <Ref> and ImageNet experiments in <Ref>).

In <Ref>, we assume $(\sigma,\Pcon)$-concentration with $\sigma=1$, $\Pcon=0.5$ and different sampling number $N$s.
We further assume $P_A = f^\gP(\vx_0)_{y_0} = 0.6$ as the true-class prediction probability under $\gP$.
In <Ref>, we take the model weights trained by Salman et al., 2019 on ImageNet and apply generalized Gaussian smoothing with $d/2 - k = 4$ and $\sigma = 0.50$.
We uniformly pick $100$ samples from the test set and compute $(1-\alpha)$-confident certified radius for each sample.
We report certified accuracy (under different $\ell_2$ radius $r$) that is the fraction of certifiably correctly classified samples by the smoothed classifier.

When the sampling error and confidence interval come into play, they quickly suppress the $\Omega(\sqrt d)$ growth rate of certified radius (blue curve)  as shown in <Ref>.
Nonetheless, still certifies a larger radius than the standard Neyman-Pearson method and increasing the sampling number further enlarges the gap.

We consider another relaxed version of concentration property in <Ref>, where still provides significantly tighter robustness certification than Neyman-Pearson.

§ COMPUTATIONAL METHOD

Overview of computational method.

The theoretical analysis in <Ref> implies that additional smoothing distribution $\gQ$ helps to tighten the robustness certification over standard Neyman-Pearson-based certification significantly.
In this section, we propose an efficient computational method to compute this tight certified robust radius $\rDSRS$ (see <Ref>) when $\gP$ is generalized Gaussian and $\gQ$ is truncated $\gP$ as suggested by <Ref>.

Compared with the classical certification for randomized smoothing or its variants (cf. [Kumar et al., 2020]), incorporating additional information raises a big challenge:
the Neyman-Pearson lemma ([Neyman & Pearson, 1933]) can no longer be served as the foundation of the certification algorithm due to its incapability to handle the additional information.

Thus, we propose a novel computational method by formalizing robustness certification as a constrained optimization problem and proving its strong duality ($\S$<ref>).
Then, we propose an efficient algorithm to solve this specific dual optimization problem considering sampling error.
The detailed algorithm can be found in <Ref> in <Ref>:
1) we first perform a binary search on the certified radius $r$ to determine the maximum radius that we can certify;
for current $r$, we determine the smoothed prediction confidence $P_A$ and $Q_A$ from the confidence intervals of predicting the true class ($\S$<ref>);
3) then, for current $r$ we solve the dual problem by quick binary search for dual variables $\lambda_1$ and $\lambda_2$ (see <Ref>) along with numerical integration ($\S$<ref>).
To guarantee the soundness of numerical-integration-based certification, we take the maximum possible error into account during the binary search.
We will discuss further extensions in $\S$<ref>.

§.§ as Constrained Optimization

We first formulate the robustness certification as a constrained optimization problem and then show several foundational properties of the problem.

Following the notation of <Ref>, from the given base classifier $F_0$, we can use Monte-Carlo sampling to obtain

\begin{equation}
% \small
\label{eq:first-appear-p-a-q-a}
P_A = f_0^\gP(\vx_0)_{y_0},\quad Q_A = f_0^\gQ(\vx_0)_{y_0}.
\end{equation}

In $\S$<ref> we will discuss how to handle confidence intervals of $P_A$ and $Q_A$.
For now, we assume $P_A$ and $Q_A$ are fixed.

Given perturbation vector $\vdelta\in \sR^d$, to test whether smoothed classifier $\wF_0^\gP$ still predicts true label $y_0$, we only need to check whether the prediction probability $f_0^\gP(\vx_0 + \vdelta)_{y_0} > 0.5$. This can be formulated as a constrained optimization problem $(\tC)$:

\begin{align}
\hspace{-0.9em} \underset{f}{\mathrm{minimize}} \quad & \E_{\vepsilon\sim\gP} [f(\vepsilon + \vdelta)] \label{eq:primal-opt} \\
\mathrm{s.t.}\quad & \E_{\vepsilon\sim\gP} [f(\vepsilon)] = P_A,  \quad % \\
% &
\E_{\vepsilon\sim\gQ} [f(\vepsilon)] = Q_A, \label{eq:primal-pa-qa} \\
& 0 \le f(\vepsilon) \le 1 \quad \forall \vepsilon\in\sR^d. \label{eq:primal-f-cons}
\end{align}

$(\tC)$ seeks for the minimum possible $f^\gP(\vx_0 + \vdelta)_{y_0}$ given <Ref>'s constraint.
Concretely, we let $f$ represent whether the base classifier predicts label $y_0$: $f(\cdot) = \1[F(\cdot\, + \vx_0) = y_0]$, and accordingly impose $f \in [0,\,1]$ in <Ref>.
Then, <Ref> unfold $f^\gP(\vx_0 + \vdelta)_{y_0}, f^\gP(\vx_0)_{y_0},$ and $f^\gQ(\vx_0)_{y_0}$ respectively and impose <Ref>'s constraint.

We let $\tC_\vdelta(P_A,\,Q_A)$ denote the optimal value of <Ref> when feasible.
Thus, under norm $p$,  to certify the robustness within radius $r$, we only need to check whether

\begin{equation}
% \small
\forall \vdelta, \|\vdelta\|_p < r \Rightarrow \tC_\vdelta(P_A,\,Q_A) > 0.5.
\label{eq:check-cond}
\end{equation}

This formulation yields the tightest robustness certification given information from $\gP$ and $\gQ$ under the binary setting.
Under the multiclass setting, there are efforts towards tighter certification by using “$>$ maximum over other classes” instead of “$> 0.5$” in <Ref> [Dvijotham et al., 2020].
For saving the sampling cost and also to follow the convention [Cohen et al., 2019, Yang et al., 2020, Jeong & Shin, 2020, Zhai et al., 2020], we mainly consider “$> 0.5$” for multiclass setting, and extension to the other form is straightforward.

Since our choices of $\gP$ and $\gQ$ (standard/generalized (truncated) Gaussian) are isotropic and centered around origin, when certifying radius $r$, for $\ell_2$ certification we only need to test $\tC_\vdelta(P_A,\,Q_A) > 0.5$ with $\vdelta = (r,\,0,\,\dots,\,0)^\T$;
for $\ell_\infty$ we only need to divide $\ell_2$ radius by $\sqrt{d}$.
This trick can also be extended for $\ell_1$ case [Zhang et al., 2020].

Directly solving $(\tC)$ is challenging. Thus, we construct the Lagrangian dual problem $(\tD)$:

\begin{equation}
\hspace{-0.2em} \underset{\lambda_1,\,\lambda_2 \in \sR}{\mathrm{maximize}} \, \Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda_1 p(\vepsilon+\vdelta) + \lambda_2 q(\vepsilon+\vdelta)]
\label{eq:dual-opt}
\end{equation}

\begin{equation}
\begin{aligned}
\mathrm{s.t.} & \Pr_{\vepsilon\sim\gP} [p(\vepsilon-\vdelta) < \lambda_1p(\vepsilon) + \lambda_2q(\vepsilon)] = P_A, \\
& \Pr_{\vepsilon\sim\gQ} [p(\vepsilon-\vdelta) < \lambda_1p(\vepsilon) + \lambda_2q(\vepsilon)] = Q_A.
\end{aligned}
\end{equation}

In <Ref>,  $p(\cdot)$ and $q(\cdot)$ are the density functions of distributions $\gP$ and $\gQ$ respectively.
We let $\tD_\vdelta(P_A,\, Q_A)$ denote the optimal objective value to <Ref> when it is feasible.

For given $\vdelta\in \sR^d$, $P_A$, and $Q_A$,
if $\tC$ and $\tD$ are both feasible, then $\tC_\vdelta(P_A,\, Q_A) = \tD_\vdelta(P_A,\, Q_A)$.

The theorem states the strong duality between $(\tC)$ and $(\tD)$.
We defer the proof to <Ref>. The proof is based on min-max inequality and feasibility condition of $(\tD)$.
Intuitively, we can view $(\tC)$, a functional optimization over $f$, as a linear programming (LP) problem over infinite number of variables $\{f(\vx): \vx\in\R^d\}$
so that the strong duality holds, which
guarantees the tightness of in the primal space.

§.§ Dealing with Confidence Intervals

It is practically intractable to know the exact $P_A$ and $Q_A$ in <Ref> by only querying the model's prediction for finite times.
The common practice is using Monte-Carlo sampling, which gives confidence intervals of $P_A$ and $Q_A$ with a predefined confidence level $1 - \alpha$.

Suppose we have confidence intervals $[\underline{P_A},\,\overline{P_A}]$
and $[\underline{Q_A},\,\overline{Q_A}]$.
To derive a sound certification, we need to certify that for any $P_A \in [\underline{P_A},\,\overline{P_A}]$ and any $Q_A \in [\underline{Q_A},\,\overline{Q_A}]$, $\tC_\vdelta(P_A,\,Q_A) > 0.5$.
Given the infinite number of possible $P_A$ and $Q_A$, the brute-force method is intractable. Here, without computing $\tC_\vdelta$, we show how to solve

\begin{equation}
% \small
\hspace{-0.215em}
(P_A,\,Q_A) = \argmin_{P_A' \in [\underline{P_A},\,\overline{P_A}],\\ Q_A' \in [\underline{Q_A},\,\overline{Q_A}]} \tC_\vdelta(P_A',\,Q_A').
\label{eq:preprocess-q-a-q-a}
\end{equation}

If solved $P_A$ and $Q_A$ satify $\tC_\vdelta(P_A,\,Q_A) > 0.5$,
then for any $P_A$ and $Q_A$ within the confidence intervals, we can certify the robustness against perturbation $\vdelta$.
We observe the following two properties of $\tC_\vdelta$.

$\tC_\vdelta(\cdot,\,\cdot)$ is convex in the feasible region.

With respect to $x \in [0,\,1]$,
functions $x \mapsto \min_{y} \tC_\vdelta(x,\,y)$ and $x \mapsto \argmin_{y} \tC_\vdelta(x,\,y)$ are monotonically non-decreasing.
Similarly, with respect to $y \in [0,\,1]$,
functions $y \mapsto \min_{x} \tC_\vdelta(x,\,y)$ and $y \mapsto \argmin_{x} \tC_\vdelta(x,\,y)$ are monotonically non-decreasing.

These two propositions characterize the landscape of $\tC_\vdelta(\cdot,\cdot)$—convex and monotonically non-decreasing along both $x$ and $y$ axes.
Thus, desired ($P_A$, $Q_A$) (location of minima within the bounded box) lies on the box boundary, and we only need to compute the location of boundary-line-sliced minima and compare it with box constraints to solve <Ref>.
Formally, we propose an efficient algorithm (<Ref>, omitted to <Ref>) to solve $(P_A,Q_A)$.

If <Ref> is feasible,
the $P_A$ and $Q_A$ returned by <Ref> solve <Ref>.

The above results are proved in <Ref>.
On a high level, we prove <Ref> by definition; we prove <Ref> via a reduction to classical Neyman-Pearson-based certification and analysis of this reduced problem; and we prove <Ref>
based on <Ref> along with exhaustive and nontrivial analyses of all possible cases.

§.§ Solving the Dual Problem

After the smoothed prediction confidences $P_A$ and $Q_A$ are determined from the confidence intervals, now we solve the dual problem $\tD_\vdelta(P_A,Q_A)$ as defined in <Ref>.
We solve the problem based on the following theorem:

In $\tD_\vdelta(P_A,Q_A)$, let $r = \|\vdelta\|_2$,
when $\gP = \gNg(k,\sigma)$ and $\gQ = \gNgtrunc(k,T,\sigma)$, let $\sigma' := \sqrt{d/(d-2k)}$ and let $\nu := \GammaCDF_{d/2-k}(T^2/(2\sigma'^2))$,

\begin{align*}
R(\lambda_1,\lambda_2) & := \Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda_1 p(\vepsilon + \vdelta) + \lambda_2 q(\vepsilon + \vdelta)] \\
= & \left\{ \begin{array}{lr}
\E_{t\sim\Gamma(d/2-k,1)} u_1(t), & \lambda_1 \le 0 \\
\E_{t\sim\Gamma(d/2-k,1)} u_1(t) + u_2(t), & \lambda_1 > 0
\end{array} \right. \mathrm{where} \\
u_1(t) = & \BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{\min\{T^2, 2\sigma'^2 k W(\frac{t}{k}e^{\frac t k} (\lambda_1 + \nu \lambda_2)^{\frac 1 k}) \}}{4r\sigma'\sqrt{2t}} \right. \\
& \hspace{3em} \left. -
\dfrac{(\sigma'\sqrt{2t} - r)^2}
\right), \\
u_2(t) = \hspace{1em} & \hspace{-1em} \max\left\{
\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{2\sigma'^2k W(\frac{t}{k} e^{\frac t k} \lambda_1^{\frac 1 k}) - (\sigma'\sqrt{2t} - r)^2}{4r\sigma'\sqrt{2t}}
\right)
\right. \\
& \hspace{3em} \left. -\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma'\sqrt{2t} - r)^2}{4r\sigma'\sqrt{2t}}
\right), 0 \right\}, \\
P(\lambda_1,\lambda_2) & := \Pr_{\vepsilon\sim\gP} [p(\vepsilon-\vdelta) < \lambda_1 p(\vepsilon) + \lambda_2 q(\vepsilon)] \\
& \hspace{-2em} = \E_{t\sim\Gamma(d/2-k,1)} \left\{ \begin{array}{lr}
u_3(t,\lambda_1), & t \ge T^2 / (2\sigma'^2) \\
u_3(t,\lambda_1 + \nu\lambda_2), & t < T^2 / (2\sigma'^2).
\end{array} \right. \mathrm{where} \\
u_3(t,\lambda) & = \BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{(r + \sigma'\sqrt{2t})^2}{4r\sigma'\sqrt{2t}}
- \dfrac{2k\sigma'^2 W(\frac{t}{k}e^{\frac t k}\lambda^{-\frac 1 k})}{4r\sigma'\sqrt{2t}} \right), \\
Q(\lambda_1,\lambda_2) & := \Pr_{\vepsilon\sim\gQ} [p(\vepsilon-\vdelta) < \lambda_1 p(\vepsilon) + \lambda_2 q(\vepsilon)] \\
& = \nu \E_{t\sim\Gamma(d/2-k,1)} u_3(t,\lambda_1 + \nu\lambda_2) \cdot \1[t \le T^2/(2\sigma'^2)].
\end{align*}

In above equations, $\Gamma(d/2-k,1)$ is gamma distribution and $\GammaCDF_{d/2-k}$ is its CDF,
$\BetaCDF_{\frac{d-1}{2}}$ is the CDF of distribution $\Beta(\frac{d-1}{2}, \frac{d-1}{2})$,
and $W$ is the principal branch of Lambert $W$ function.

When $\gP$ is standard Gaussian and $\gQ$ is truncated standard Gaussian, we derive similar expressions as detailed in <Ref>.
We prove <Ref> in <Ref>. The proof extends the level-set sliced integration and results from [Yang et al., 2020].
With the theorem, we can rewrite the dual problem $\tD_\vdelta(P_A,Q_A)$ as

\begin{equation}
\label{eq:dual-2} %\small
% \begin{aligned}
\underset{\lambda_1,\,\lambda_2 \in \sR}{\mathrm{max}}\,  R(\lambda_1,\,\lambda_2) %\\
\,
\mathrm{s.t.} \,  P(\lambda_1,\,\lambda_2) = P_A,\, Q(\lambda_1,\,\lambda_2) = Q_A,
% \end{aligned}
\end{equation}

Given concrete $\lambda_1$ and $\lambda_2$, from the theorem, these function values $P(\lambda_1,\lambda_2)$, $Q(\lambda_1,\lambda_2)$, and $R(\lambda_1,\lambda_2)$ can be easily computed with one-dimensional numerical integration using  package.

Now, solving $\tD_\vdelta(P_A,Q_A)$ reduces to finding dual variables $\lambda_1$ and $\lambda_2$ such that $P(\lambda_1,\lambda_2) = P_A$ and $Q(\lambda_1,\lambda_2) = Q_A$.
Generally, we find that there is only one unique feasible pair $(\lambda_1,\lambda_2)$ for <Ref>, so finding out such a pair is sufficient.
We prove the uniqueness and discuss how we deal with edge cases where multiple feasible pairs exist in <Ref>.

Normally, such solving process is expensive.
However, we find a particularly efficient method to solve $\lambda_1$ and $\lambda_2$ and the algorithm description is in <Ref> (in <Ref>).
At a high level, from <Ref>, we observe that $Q(\lambda_1,\lambda_2)$ is determined only by the sum $(\lambda_1 + \nu\lambda_2)$ and non-decreasing w.r.t. this sum.
Therefore, we apply binary search to find out $(\lambda_1 + \nu\lambda_2)$ that satisfies $Q(\lambda_1,\lambda_2) = Q_A$.
Then, we observe that

\begin{equation*}
\scriptsize
P(\lambda_1,\lambda_2) - \frac{Q(\lambda_1,\lambda_2)}{\nu} = \overbrace{\underset{\rt\sim\Gamma(d/2-k,1)}{\E} u_3(\rt,\lambda_1) \cdot \1\left[\rt \ge \frac{T^2}{2\sigma'^2}\right]}^{:=h(\lambda_1)}.
\vspace{-2mm}
\end{equation*}

Thus, we need to find $\lambda_1$ such that $h(\lambda_1) = P_A - Q_A / \nu$.
We observe that $h(\lambda_1)$ is non-decreasing w.r.t. $\lambda_1$, and we use binary search to solve $\lambda_1$.
Combining with the value of $(\lambda_1 + \nu\lambda_2)$, we also obtain $\lambda_2$.
We lastly leverage numerical integration to compute $R(\lambda_1,\lambda_2)$ following <Ref> to solve the dual problem $\tD_\vdelta(P_A,Q_A)$.

Practical Certification Soundness.
As a practical certification method, we need to guarantee the certification soundness in the presence of numerical error.
In , there are two sources of numerical error: numerical integration error when computing $P(\lambda_1,\lambda_2)$, $Q(\lambda_1,\lambda_2)$, and $R(\lambda_1,\lambda_2)$, and the finite precision of binary search on $\lambda_1$ and $\lambda_2$.
For numerical integration, we notice that typical numerical integration packages such as  support setting an absolute error threshold $\Delta$ and raising warnings when such threshold cannot be reached.
We set the absolute threshold $\Delta=1.5\times 10^{-8}$, and abstain when the threshold cannot be reached (which never happens in our experimental evaluation).
Then, when computing $P$, $Q$, and $R$, suppose the numerical value is $v$, we use the lower bound $(v-\Delta)$ and upper bound $(v+\Delta)$ in the corresponding context to guarantee the soundness.
For the finite precision in binary search, we use the left endpoint or the right endpoint of the final binary search interval to guarantee soundness.
For example, we use the left endpoint of $\lambda_1$ in $R$ computation, and use the left endpoint of $(\lambda_1 + \nu\lambda_2)$ minus right endpoint of $\lambda_1$ to get the lower bound of $\lambda_2$ to use in $R$ computation.
As a result, we always get an under-estimation of $R$ so the certification is sound.
Further discussion is in <Ref>.

To this point, we have introduced the computational method.
Complexity and efficiency analysis is omitted to <Ref>.
Implementation details are in <Ref>.

§.§ Extensions

We mainly discussed computational method for generalized Gaussian $\gP$ and truncated generalized Gaussian $\gQ$ under $\ell_2$ norm.
Can we extend it to other settings?
Indeed, is a general framework.
In appendices, we show following extensions:
(1) for generalized Gaussian with different variances as $\gP$ and $\gQ$ (in <Ref>);
(2) for other $\ell_p$ norms (in <Ref>);
and (3) that leverages other forms of additional information covering gradient magnitude information [Mohapatra et al., 2020, Levine et al., 2020] (in <Ref>).

§ EXPERIMENTAL EVALUATION

In this section, we systematically evaluate and demonstrate that it achieves tighter certification than the classical Neyman-Pearson-based certification against $\ell_2$ perturbations on MNIST, CIFAR-10, and ImageNet.
We focus on $\ell_2$ certification because additive randomized smoothing is not optimal for other norms (e.g., $\ell_1$ [Levine & Feizi, 2021]) or the certification can be directly translated from $\ell_2$ certification (e.g., $\ell_\infty$ [Yang et al., 2020] and semantic transformations [Li et al., 2021]).

§.§ Experimental Setup

Smoothing Distributions.
Following <Ref>, we use generalized Gaussian $\gNg(k,\sigma)$ as smoothing distribution $\gP$ where $d/2-15 \le k < d/2$.
we set $k$ to be $d/2-12$ on MNIST, $d/2-6$ on CIFAR-10, and $d/2-4$ on ImageNet.
We use three different $\sigma$'s: $0.25$, $0.50$, and $1.00$.

In terms of the additional smoothing distribution $\gQ$,
on MNIST and CIFAR-10, we empirically find that using generalized Gaussian with the same $k$ but different variance yields tighter robustness certification, and therefore we choose $\sigma_g$ to be $0.2$, $0.4$, and $0.8$ corresponding to $\gP$'s $\sigma$ being $0.25$, $0.50$, and $1.0$, respectively.
On ImageNet, the concentration property (see <Ref>) is more pronounced (detail study in <Ref>) and thus we use truncated generalized Gaussian $\gNgtrunc(k,T,\sigma)$ as $\gQ$.
We apply a simple but effective algorithm as explained in <Ref> to determine hyperparameter $T$ in $\gNgtrunc(k,T,\sigma)$.

Models and Training.
We consider three commonly-used or state-of-the-art training methods:
Gaussian augmentation [Cohen et al., 2019],
Consistency [Jeong & Shin, 2020],
and SmoothMix [Jeong et al., 2021].
We follow the default model architecture on each dataset respectively.
We train the models with augmentation noise sampled from the corresponding generalized Gaussian smoothing distribution $\gP$.
More training details can be found in <Ref>.

We consider the Neyman-Pearson-based certification method as the baseline.
This certification is widely used and is the tightest given only prediction probability under $\gP$ [Cohen et al., 2019, Yang et al., 2020, Jeong & Shin, 2020, Li et al., 2021].
We remark that although there are certification methods that leverage more information, to the best of our knowledge, they are not visibly better than the Neyman-Pearson-based method on $\ell_2$ certification under practical sampling number ($10^5$).
More comparisons in <Ref> show is also better than these baselines.

For both baseline and ,
following the convention, the certification confidence is $1 - \alpha = 99.9\%$, and we use $10^5$ samples for estimating $P_A$ and $Q_A$.
Neyman-Pearson certification does not use the information from additional distribution, and all $10^5$ samples are used to estimate the interval of $P_A$.
In , we use $5\times 10^4$ samples to estimate the interval of $P_A$ with confidence $1-\frac{\alpha}{2} = 99.95\%$ and the rest $5\times 10^4$ samples for $Q_A$ with the same confidence.
By union bound, the whole certification confidence is $99.9\%$.

We uniformly draw $1000$ samples from the test set and report certified robust accuracy (under each $\ell_2$ radius $r$) that is the fraction of samples that are both correctly classified and have certified robust radii larger than or equal to $r$.
Under each radius $r$, we report the highest certified robust accuracy among the three variances $\sigma \in \{0.25,0.50,1.00\}$ following [Cohen et al., 2019, Salman et al., 2019].
We also report evaluation results with ACR metric [Zhai et al., 2020] in <Ref>.

§.§ Evaluation Results

We show results in <Ref>.
The corresponding curves and separated tables for each variance $\sigma$ are in <Ref>.

For almost all models and radii $r$, yields significantly higher certified accuracy.
For example, for Gaussian augmented models, when $r=2.0$, on MNIST the robust accuracy increases from $25.5\%$ to $34.1\%$ ($+8.6\%$),
on CIFAR-10 from $4.1\%$ to $6.4\%$ ($+2.3\%$),
and on ImageNet from $13.7\%$ to $18.7\%$ ($+5.0\%$).
On average, on MNIST the improvements are around $6\%$ - $9\%$;
on CIFAR-10, the improvements are around $1.5\%$ - $3\%$;
and on ImageNet the improvements are around $2\%$ - $5\%$.
Thus, can be used in conjunction with different training approaches and provides consistently tighter robustness certification.

The improvements in the robust radius are not as substantial as those in <Ref> (which is around $2\times$).
We investigate the reason in <Ref>.
In summary, the model in <Ref> is trained with standard Gaussian smoothing augmentation and smoothed with generalized Gaussian.
The models in this section are trained with generalized Gaussian augmentation.
Such training gives higher certified robustness, but in the meantime, gives more advantage to Neyman-Pearson-based certification.

This finding implies that there may be a large space for exploring training approaches that favor certification since all existing training methods are designed for Neyman-Pearson-based certification.
Nevertheless, even with these “unsuitable” training methods, still achieves significantly tighter robustness certification than the baseline.

On the other hand, all the above results are restricted to generalized Gaussian smoothing.
We still observe that standard Gaussian smoothing combined with strong training methods [Salman et al., 2019, Jeong & Shin, 2020] and Neyman-Pearson certification (the SOTA setting) yields similar or slightly higher certified robust accuracy than generalized Gaussian smoothing even with certification.
Though with its suitable generalized Gaussian smoothing does not achieve SOTA certified robustness yet, given the theoretical advantages, we believe that with future tailored training approaches, with generalized Gaussian smoothing can bring strong certified robustness.
More discussion is in <Ref>.

Ablation Studies.
We present several ablation studies in the appendix. and verify:
(1) Effectiveness of our simple heuristic for selecting hyperparameter for $\gQ$:
We propose a simple heuristic to select the hyperparameter $T$ in smoothing distribution $\gQ = \gNgtrunc(k,T,\sigma)$.
In <Ref>, we propose a gradient-based optimization method to select such $\gQ$.
We find that our simple heuristic has similar performance compared to the more complex optimization method but is more efficient.
(2) Comparison of different types of $\gQ$:
by choosing different types of $\gQ$ distributions (truncated Gaussian or Gaussian with different variance), has different performance as mentioned in <Ref>.
In <Ref>, we investigate the reason.
In summary,
when concentration property (see <Ref>) is better satisfied, using truncated Gaussian as $\gQ$ is better; otherwise, using Gaussian with different variance is better.

§ CONCLUSION

We propose a general framework that exploits information based on an additional smoothing distribution to tighten the robustness certification.
We theoretically analyze and compare classical Neyman-Pearson and certification, showing that has the potential to break the curse of dimensionality of randomized smoothing.

§ ACKNOWLEDGEMENTS

We thank the anonymous reviewers for their constructive feedback.
This work is partially supported by
NSF grant No.1910100, NSF CNS No.2046726, C3 AI, and the Alfred P. Sloan Foundation.

[Albarghouthi, 2021]
Albarghouthi, A.
Introduction to neural network verification.
arXiv preprint arXiv:2109.10317, 2021.

[Alfarra et al., 2020]
Alfarra, M., Bibi, A., Torr, P. H., and Ghanem, B.
Data dependent randomized smoothing.
arXiv preprint arXiv:2012.04351, 2020.

[Athalye et al., 2018]
Athalye, A., Carlini, N., and Wagner, D. A.
Obfuscated gradients give a false sense of security: Circumventing
defenses to adversarial examples.
In Dy, J. G. and Krause, A. (eds.), Proceedings of the 35th
International Conference on Machine Learning, ICML 2018,
Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of
Proceedings of Machine Learning Research, pp.  274–283. PMLR,

[Awasthi et al., 2020]
Awasthi, P., Jain, H., Rawat, A. S., and Vijayaraghavan, A.
Adversarial robustness via robust low rank representations.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  11391–11403. Curran Associates, Inc.,

[Blum et al., 2020]
Blum, A., Dick, T., Manoj, N., and Zhang, H.
Random smoothing might be unable to certify $\ell_\infty$ robustness
for high-dimensional images.
Journal of Machine Learning Research, 21:0
211:1–211:21, 2020.

[Carlini & Wagner, 2017]
Carlini, N. and Wagner, D. A.
Towards evaluating the robustness of neural networks.
In 2017 IEEE Symposium on Security and Privacy, SP 2017,
San Jose, CA, USA, May 22-26, 2017, pp.  39–57. IEEE Computer Society,

[Carmon et al., 2019]
Carmon, Y., Raghunathan, A., Schmidt, L., Duchi, J. C., and Liang, P.
Unlabeled data improves adversarial robustness.
In Wallach, H. M., Larochelle, H., Beygelzimer, A.,
d'Alché-Buc, F., Fox, E. B., and Garnett, R. (eds.), Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019,
Vancouver, BC, Canada, pp.  11190–11201. Curran Associates, Inc., 2019.

[Chernoff & Scheffe, 1952]
Chernoff, H. and Scheffe, H.
A generalization of the neyman-pearson fundamental lemma.
The Annals of Mathematical Statistics, 230
(4):0 213–225, 1952.

[Chiang et al., 2020]
Chiang, P., Curry, M. J., Abdelkader, A., Kumar, A., Dickerson, J., and
Goldstein, T.
Detection as regression: Certified object detection with median
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  1275–1286. Curran Associates, Inc.,

[Cohen et al., 2019]
Cohen, J. M., Rosenfeld, E., and Kolter, J. Z.
Certified adversarial robustness via randomized smoothing.
In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of
the 36th International Conference on Machine Learning, ICML 2019, 9-15 June
2019, Long Beach, California, USA, volume 97 of Proceedings of
Machine Learning Research, pp.  1310–1320. PMLR, 2019.

[Dvijotham et al., 2020]
Dvijotham, K. D., Hayes, J., Balle, B., Kolter, J. Z., Qin, C., György,
A., Xiao, K., Gowal, S., and Kohli, P.
A framework for robustness certification of smoothed classifiers
using f-divergences.
In 8th International Conference on Learning Representations,
ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020. OpenReview.net, 2020.

[Eiras et al., 2021]
Eiras, F., Alfarra, M., Kumar, M. P., Torr, P. H. S., Dokania, P. K., Ghanem,
B., and Bibi, A.
ANCER: anisotropic certification via sample-wise volume
arXiv preprint arXiv:2107.04570, 2021.

[Eykholt et al., 2018]
Eykholt, K., Evtimov, I., Fernandes, E., Li, B., Rahmati, A., Xiao, C.,
Prakash, A., Kohno, T., and Song, D.
Robust physical-world attacks on deep learning visual classification.
In 2018 IEEE Conference on Computer Vision and Pattern
Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, pp.  1625–1634. Computer Vision Foundation / IEEE Computer Society, 2018.

[Goodfellow et al., 2015]
Goodfellow, I. J., Shlens, J., and Szegedy, C.
Explaining and harnessing adversarial examples.
In Bengio, Y. and LeCun, Y. (eds.), 3rd International
Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May
7-9, 2015, Conference Track Proceedings, 2015.

[Gowal et al., 2019]
Gowal, S., Dvijotham, K., Stanforth, R., Bunel, R., Qin, C., Uesato, J.,
Arandjelovic, R., Mann, T. A., and Kohli, P.
Scalable verified training for provably robust image classification.
In 2019 IEEE/CVF International Conference on Computer Vision,
ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019, pp.  4841–4850. IEEE, 2019.

[Hayes, 2020]
Hayes, J.
Extensions and limitations of randomized smoothing for robustness
In 2020 IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR Workshops 2020, Seattle, WA, USA, June 14-19, 2020, pp.   3413–3421. Computer Vision Foundation / IEEE, 2020.

[He et al., 2016]
He, K., Zhang, X., Ren, S., and Sun, J.
Deep residual learning for image recognition.
In 2016 IEEE Conference on Computer Vision and Pattern
Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, pp.  770–778. IEEE Computer Society, 2016.

[Jeong & Shin, 2020]
Jeong, J. and Shin, J.
Consistency regularization for certified robustness of smoothed
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  10558–10570. Curran Associates, Inc.,

[Jeong et al., 2021]
Jeong, J., Park, S., Kim, M., Lee, H., Kim, D., and Shin, J.
Smoothmix: Training confidence-calibrated smoothed classifiers for
certified robustness.
In Ranzato, M., Beygelzimer, A., Dauphin, Y. N., Liang, P., and
Vaughan, J. W. (eds.), Advances in Neural Information Processing
Systems 34: Annual Conference on Neural Information Processing Systems 2021,
NeurIPS 2021, December 6-14, 2021, virtual, pp.  30153–30168. Curran
Associates, Inc., 2021.

[Kumar et al., 2020]
Kumar, A., Levine, A., Feizi, S., and Goldstein, T.
Certifying confidence via randomized smoothing.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  5165–5177. Curran Associates, Inc.,

[Kumar et al., 2020]
Kumar, A., Levine, A., Goldstein, T., and Feizi, S.
Curse of dimensionality on randomized smoothing for certifiable
In Proceedings of the 37th International Conference on Machine
Learning, ICML 2020, 13-18 July 2020, Virtual Event, volume 119 of
Proceedings of Machine Learning Research, pp.  5458–5467. PMLR,

[Lécuyer et al., 2019]
Lécuyer, M., Atlidakis, V., Geambasu, R., Hsu, D., and Jana, S.
Certified robustness to adversarial examples with differential
In 2019 IEEE Symposium on Security and Privacy, SP 2019,
San Francisco, CA, USA, May 19-23, 2019, pp.  656–672. IEEE, 2019.

[Lee et al., 2019]
Lee, G., Yuan, Y., Chang, S., and Jaakkola, T. S.
Tight certificates of adversarial robustness for randomly smoothed
In Wallach, H. M., Larochelle, H., Beygelzimer, A.,
d'Alché-Buc, F., Fox, E. B., and Garnett, R. (eds.), Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019,
Vancouver, BC, Canada, pp.  4911–4922. Curran Associates, Inc., 2019.

[Levine & Feizi, 2021]
Levine, A. and Feizi, S.
Improved, deterministic smoothing for $\ell_1$ certified robustness.
In Meila, M. and Zhang, T. (eds.), Proceedings of the 38th
International Conference on Machine Learning, ICML 2021, 18-24 July 2021,
Virtual Event, volume 139 of Proceedings of Machine Learning
Research, pp.  6254–6264. PMLR, 2021.

[Levine et al., 2020]
Levine, A., Kumar, A., Goldstein, T. A., and Feizi, S.
Tight second-order certificates for randomized smoothing.
arXiv preprint arXiv:2010.10549, 2020.

[Li et al., 2019]
Li, B., Chen, C., Wang, W., and Carin, L.
Certified adversarial robustness with additive noise.
In Wallach, H. M., Larochelle, H., Beygelzimer, A.,
d'Alché-Buc, F., Fox, E. B., and Garnett, R. (eds.), Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019,
Vancouver, BC, Canada, pp.  9459–9469. Curran Associates, Inc.,

[Li et al., 2020]
Li, H., Xu, X., Zhang, X., Yang, S., and Li, B.
QEBA: query-efficient boundary-based blackbox attack.
In 2020 IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020, pp.  1218–1227. Computer Vision Foundation / IEEE, 2020a.

[Li et al., 2019]
Li, L., Zhong, Z., Li, B., and Xie, T.
Robustra: Training provable robust neural networks over reference
adversarial space.
In Kraus, S. (ed.), Proceedings of the Twenty-Eighth
International Joint Conference on Artificial Intelligence, IJCAI 2019,
Macao, China, August 10-16, 2019, pp.  4711–4717. ijcai.org,

[Li et al., 2020]
Li, L., Xie, T., and Li, B.
Sok: Certified robustness for deep neural networks.
arXiv preprint arXiv:2009.04131, 2020b.

[Li et al., 2021]
Li, L., Weber, M., Xu, X., Rimanic, L., Kailkhura, B., Xie, T., Zhang, C., and
Li, B.
TSS: transformation-specific smoothing for robustness
In Kim, Y., Kim, J., Vigna, G., and Shi, E. (eds.), 2021 ACM
SIGSAC Conference on Computer and Communications Security, CCS 2021,
Virtual Event, Republic of Korea, November 15 - 19, 2021, pp.  535–557.
ACM, 2021.

[Liu et al., 2021]
Liu, C., Arnon, T., Lazarus, C., Strong, C. A., Barrett, C. W., and
Kochenderfer, M. J.
Algorithms for verifying deep neural networks.
Foundations and Trends in Optimization, 40
(3-4):0 244–404, 2021.

[Lozier, 2003]
Lozier, D. W.
NIST digital library of mathematical functions.
Annals of Mathematics and Artificial Intelligence, 380
(1-3):0 105–119, 2003.

[Mirman et al., 2018]
Mirman, M., Gehr, T., and Vechev, M. T.
Differentiable abstract interpretation for provably robust neural
In Dy, J. G. and Krause, A. (eds.), Proceedings of the 35th
International Conference on Machine Learning, ICML 2018,
Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of
Proceedings of Machine Learning Research, pp.  3575–3583. PMLR,

[Mohapatra et al., 2020]
Mohapatra, J., Ko, C., Weng, T., Chen, P., Liu, S., and Daniel, L.
Higher-order certification for randomized smoothing.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  4501–4511. Curran Associates, Inc.,

[Neyman & Pearson, 1933]
Neyman, J. and Pearson, E. S.
On the problem of the most efficient tests of statistical hypotheses.
Philosophical Transactions of the Royal Society of London.
Series A, Containing Papers of a Mathematical or Physical Character,
231:0 289–337, 1933.

[Qiu et al., 2020]
Qiu, H., Xiao, C., Yang, L., Yan, X., Lee, H., and Li, B.
Semanticadv: Generating adversarial examples via
attribute-conditioned image editing.
In Vedaldi, A., Bischof, H., Brox, T., and Frahm, J. (eds.),
16th European Conference on Computer Vision, ECCV 2020, Glasgow, UK,
August 23-28, 2020, Proceedings, Part XIV, volume 12359 of Lecture
Notes in Computer Science, pp.  19–37. Springer, 2020.

[Raghunathan et al., 2018]
Raghunathan, A., Steinhardt, J., and Liang, P.
Certified defenses against adversarial examples.
In 6th International Conference on Learning Representations,
ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track
Proceedings. OpenReview.net, 2018.

[Salman et al., 2019]
Salman, H., Li, J., Razenshteyn, I. P., Zhang, P., Zhang, H., Bubeck, S., and
Yang, G.
Provably robust deep learning via adversarially trained smoothed
In Wallach, H. M., Larochelle, H., Beygelzimer, A.,
d'Alché-Buc, F., Fox, E. B., and Garnett, R. (eds.), Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019,
Vancouver, BC, Canada, pp.  11289–11300. Curran Associates, Inc., 2019.

[Schuchardt et al., 2022]
Schuchardt, J., Wollschläger, T., Bojchevski, A., and Günnemann, S.
Localized randomized smoothing for collective robustness
certification, 2022.
URL <https://openreview.net/forum?id=mF122BuAnnW>.

[Súkeník et al., 2021]
Súkeník, P., Kuvshinov, A., and Günnemann, S.
Intriguing properties of input-dependent randomized smoothing.
arXiv preprint arXiv:2110.05365, 2021.

[Szegedy et al., 2014]
Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow,
I. J., and Fergus, R.
Intriguing properties of neural networks.
In Bengio, Y. and LeCun, Y. (eds.), 2nd International
Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April
14-16, 2014, Conference Track Proceedings, 2014.

[Tramèr et al., 2020]
Tramèr, F., Carlini, N., Brendel, W., and Madry, A.
On adaptive attacks to adversarial example defenses.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  1633–1645. Curran Associates, Inc.,

[Wang et al., 2021]
Wang, B., Xu, C., Wang, S., Gan, Z., Cheng, Y., Gao, J., Awadallah, A. H., and
Li, B.
Adversarial GLUE: A multi-task benchmark for robustness
evaluation of language models.
In Vanschoren, J. and Yeung, S. (eds.), Proceedings of the
Neural Information Processing Systems Track on Datasets and Benchmarks 1,
NeurIPS Datasets and Benchmarks 2021, December 2021, virtual, volume 1,

[Wong & Kolter, 2018]
Wong, E. and Kolter, J. Z.
Provable defenses against adversarial examples via the convex outer
adversarial polytope.
In Dy, J. G. and Krause, A. (eds.), Proceedings of the 35th
International Conference on Machine Learning, ICML 2018,
Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of
Proceedings of Machine Learning Research, pp.  5283–5292. PMLR,

[Wu et al., 2021]
Wu, Y., Bojchevski, A., Kuvshinov, A., and Günnemann, S.
Completing the picture: Randomized smoothing suffers from the curse
of dimensionality for a large family of distributions.
In Banerjee, A. and Fukumizu, K. (eds.), The 24th International
Conference on Artificial Intelligence and Statistics, AISTATS 2021, April
13-15, 2021, Virtual Event, volume 130 of Proceedings of Machine
Learning Research, pp.  3763–3771. PMLR, 2021.

[Xu et al., 2020]
Xu, K., Shi, Z., Zhang, H., Wang, Y., Chang, K., Huang, M., Kailkhura, B., Lin,
X., and Hsieh, C.
Automatic perturbation analysis for scalable certified robustness and
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  1129–1141. Curran Associates, Inc.,

[Yang et al., 2020]
Yang, G., Duan, T., Hu, J. E., Salman, H., Razenshteyn, I. P., and Li, J.
Randomized smoothing of all shapes and sizes.
In Proceedings of the 37th International Conference on Machine
Learning, ICML 2020, 13-18 July 2020, Virtual Event, volume 119 of
Proceedings of Machine Learning Research, pp.  10693–10705. PMLR,

[Yang et al., 2022]
Yang, Z., Li, L., Xu, X., Kailkhura, B., Xie, T., and Li, B.
On the certified robustness for ensemble models and beyond.
In 10th International Conference on Learning Representations,
ICLR 2022, Virtual Event, April 25-29, 2022. OpenReview.net, 2022.

[Zhai et al., 2020]
Zhai, R., Dan, C., He, D., Zhang, H., Gong, B., Ravikumar, P., Hsieh, C., and
Wang, L.
MACER: attack-free and scalable robust training via maximizing
certified radius.
In 8th International Conference on Learning Representations,
ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020. OpenReview.net, 2020.

[Zhang et al., 2020]
Zhang, D., Ye, M., Gong, C., Zhu, Z., and Liu, Q.
Black-box certification with randomized smoothing: A functional
optimization based framework.
In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H.
(eds.), Advances in Neural Information Processing Systems 33: Annual
Conference on Neural Information Processing Systems 2020, NeurIPS 2020,
December 6-12, 2020, virtual, pp.  2316–2326. Curran Associates, Inc.,

[Zhang et al., 2020]
Zhang, H., Chen, H., Xiao, C., Gowal, S., Stanforth, R., Li, B., Boning, D. S.,
and Hsieh, C.
Towards stable and efficient training of verifiably robust neural
In 8th International Conference on Learning Representations,
ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020. OpenReview.net,

[Zhang et al., 2021]
Zhang, J., Li, L., Li, H., Zhang, X., Yang, S., and Li, B.
Progressive-scale boundary blackbox attack via projective gradient
In Meila, M. and Zhang, T. (eds.), Proceedings of the 38th
International Conference on Machine Learning, ICML 2021, 18-24 July 2021,
Virtual Event, volume 139 of Proceedings of Machine Learning
Research, pp.  12479–12490. PMLR, 2021.

§ NEYMAN-PEARSON CERTIFICATION

The Neyman-Pearson-based robustness certification is the tightest certification given only prediction probability under $\gP$ [Cohen et al., 2019].
This certification and its equivalent variants are widely used for randomized smoothing.

We use $\rNP$ to represent the certified radius from the Neyman-Pearson-based method.

If the smoothing distribution $\gP$ is standard Gaussian, the following proposition gives the closed-form certified robust radius derived from the Neyman-Pearson lemma [Neyman & Pearson, 1933].

Under $\ell_2$ norm, given input $\vx_0 \in \sR^d$ with true label $y_0$.
Let $\gP = \gN(\sigma)$ be the smoothing distribution, then Neyman-Pearson-based certification yields certified radius $\rNP = \sigma\Phi^{-1}(f^\gP(\vx_0)_{y_0})$, where $\Phi^{-1}$ is the inverse CDF of unit-variance Gaussian.

For other smoothing distributions, the concretization of the Neyman-Pearson certification method can be found in [Yang et al., 2020].

In practice, the routine is to use Monte-Carlo sampling to obtain a high-confidence interval of $f^\gP(\vx_0)_{y_0}$, which implies a high-confidence certification ($\rNP$) of robust radius.
A tighter radius can be obtained when the runner-up prediction probability is known: $\rNP' = \frac{\sigma}{2}\left(\Phi^{-1}(f^\gP(\vx_0)_{y_0}) - \max_{y\in [C]: y\neq y_0} \Phi^{-1}(f^\gP(\vx_0)_y)\right)$. However, due to efficiency concern (for $C$-way classification the sampling number needs to be more than $C$ times if using $\rNP'$ for certification instead of $\rNP$), the standard routine is to only use top-class probability and $\rNP$ <cit.>.
follows this routine.

§ ILLUSTRATION OF CONCENTRATION ASSUMPTION ON IMAGENET

Blue curves: Probability of true-prediction w.r.t. $\ell_2$ length of perturbations for a base classifier from [Salman et al., 2019] on ImageNet.
Each line corresponds to one of $100$ uniformly drawn samples from test set (detailed setup in <Ref>).
Green curve: Normalized density of $\ell_2$ noise magnitude for ImageNet standard Gaussian $\gN(\sigma)$ with $\sigma = 0.5$, which highly concentrates on $\sigma\sqrt{d}$.
Thus, for constant $\Pcon$, ($\sigma, \Pcon$)-concentration can be satisfied for a significant portion of input samples.

§ ILLUSTRATION OF UNCHANGED PA WITH INCREASING D

In the first remark of <Ref>, we mention that $P_A = f^\gP(\vx_0)_{y_0}$ does not grow simultaneously along with the increase of input dimension $d$.
In <Ref>, to illustrate this phenomenon, we plot $P_A$ histograms for $1,000$ test samples from TinyImageNet and ImageNet.
Note that TinyImageNet images are downscaled ImageNet images, so the data distribution only differs in the input dimension $d$.
As we can observe, though $d$ varies, the $P_A$ distribution is highly similar, so $P_A$ is roughly unchanged along with the increase of input dimension $d$.

$P_A$ histograms for models trained on TinyImageNet (left) and ImageNet (right) with same $\sigma=0.50$.

§ UNDER RELAXED CONCENTRATION ASSUMPTION

In main text (<Ref>), we generalize the concentration assumption by replacing the holding probability $1$ in <Ref> by probability considering sampling confidence.
In this appendix, we replace the holding probability in <Ref> by $\exp(-d^\alpha)$ for $\alpha\in \{0.1, 0.2, 0.3, 0.4, 0,5\}$.

With this relaxation,  we conduct numerical simulations using the same settings as in the main text, and the corresponding results are shown in <Ref>.
Note that some solid curves terminate when $d$ is large, which is due to the limitation of floating-point precision in numerical simulations, and we use dashed lines of the same color to plot the projected radius when $d$ is large.

When the concentration property holds with probability $\exp(-d^\alpha)$ ($0 < \alpha \le 0.5$) other than $1$, from <Ref>, we observe that $\rNP d^{\alpha / 1.18}$ predicts the certified radius of well where $\rNP$ is Neyman-Pearson certified radius.
Therefore, although the $\sqrt d$ growth rate of $\ell_2$ certified radius does not hold, the radius still increases along with the dimension $d$.
Interestingly, along with the increase of dimension $d$, the vanishing probability $\exp(-d^\alpha)$ still implies the increasing volume of adversarial examples, and smoothed classifier is still certifiably robust with increasing radius reflected by despite the increasing adversarial volume.

Tendency of certified robust radius with different input dimensions $d$ under relaxed concentration assumption:
when holding probability in <Ref> is $\exp(-d^\alpha)$ with $\alpha$ from $0.1$ to $0.5$;.
Blue line: when holding probability in <Ref> is $1$.
Dotted line: Neyman-Pearson certification.
Other solid lines: when holding probability in <Ref> is $\exp(-d^\alpha)$.
Other dashed lines: projected radius by $\rDSRS^{\text{proj}} = \rNP d^{\alpha / 1.18}$.
$\alpha \in \{0.1,0.2,0.3,0.4,0.5\}$.
Both $x$- and $y$-axes are logarithmic.

§ OMITTED DETAILS OF COMPUTATIONAL METHOD

§.§ Algorithm Description

Distributions $\gP$ and $\gQ$; $\vdelta$; $[\underline{P_A},\,\overline{P_A}]$ and $[\underline{Q_A},\,\overline{Q_A}]$
$P_A$ and $Q_A$ satisfying <Ref>
Compute $\underline{q} \gets \argmin_y \tC_\vdelta(\underline{P_A},\,y)$

$\underline{q} > \underline{Q_A}$

Compute $\underline{p} \gets \argmin_x \tC_\vdelta(x,\,\underline{Q_A})$

$(\max\{\min\{\underline{p},\,\overline{P_A}\}, \underline{P_A}\},\,\underline{Q_A})$

Determining $P_A$ and $Q_A$ from confidence intervals (see <ref>).

<Ref> is a subroutine (Line $7$) of <Ref>.
Note that Lines $1$ and $5$ of <Ref> solve the constrained optimization with only one constraint (either one of <Ref>), reducing to the well-studied and solvable Neyman-Pearson-based certification.

Note that we do not need to evaluate any value of $\tC_\vdelta$ in <Ref>.
Although $\underline{q}$ and $\underline{p}$ in the algorithm are “$\argmin_{x \text{ or } y}$” over $\tC_\vdelta$, the free choices of $x$ or $y$ leave $\tC_\vdelta$'s constrained optimization with only one constraint and then $\underline{q}$ and $\underline{p}$ can be solved by Neyman-Pearson instead of evaluating $\tC_\vdelta$ directly.

clean input $\vx_0$, base classifier $F_0$; distributions $\gP$ and $\gQ$; norm type $p$; confidence level $\alpha$; numerical integration error bound $\Delta$
Certified radius $r$

Query prediction $y_0 \gets \wF_0^\gP(\vx_0)$

Sample and estimate the intervals of smoothed confidence $[\underline{P_A}, \overline{P_A}]$ under $\gP$ and $[\underline{Q_A}, \overline{Q_A}]$ under $\gQ$ with confidence $(1-\alpha)$ following [Cohen et al., 2019]

Initialize: $r_l \gets 0, r_u \gets r_{\max}$

(*[f]Binary search on radius $r$)$r_u - r_l > \mathsf{eps}$
$r_m \gets (r_l + r_u) / 2$

$\vdelta \gets (r_m,\,0,\,\dots,\,0)^\T$
*[r]for $\ell_2$ certification with $\ell_2$ symmetric $\gP$ and $\gQ$; for $\ell_\infty$ or $\ell_1$, can be adjusted following [Zhang et al., 2020]
Determine $P_A \in [\underline{P_A}, \overline{P_A}]$ and $Q_A \in [\underline{Q_A}, \overline{Q_A}]$
*[r]See <Ref> and <Ref>
$(\lambda_1,\,\lambda_2) \gets$ DualBinarySearch($P_A,\,Q_A$)
*[r]See <Ref> and <Ref>

$v \gets R(\lambda_1,\,\lambda_2) - \Delta$
*[r]Using <Ref>
$v> 0.5$
$r_l\gets r_m$

$r_u\gets r_m$

computational method.

<Ref> is the pseudocode of the whole computational method as introduced in <Ref>.

Query access to $P(\cdot,\,\cdot)$ and $Q(\cdot,\,\cdot)$; $P_A$; $Q_A$; $\nu$; precision parameter $\epsilon$; numerical integration error bound $\Delta$
$\lambda_1$ and $\lambda_2$ satisfying constraints $P(\lambda_1,\lambda_2) = P_A, Q(\lambda_1,\lambda_2) = Q_A$ (see <Ref>)

$a^L \gets 0$, $a^U \gets M$
*[r]search for $a = \lambda_1 + \nu\lambda_2$, $M$ is a large positive number

$a^U - a^L > \epsilon$
$a^m \gets (a^L + a^U) / 2$

$Q(a^m,0) < Q_A$

$a^L \gets a^m$

$a^U \gets a^m$

Following while-loop enlarges $a^L$ and $a^U$ until $[a^L, a^U]$ covers $a^*$ such that $Q(a^*,0) = Q_A$ under numerical integration error
$(Q(a^L,0) + \Delta > Q_A)$  or  $(Q(a^U,0) - \Delta < Q_A)$
$t \gets a^U - a^L$

$a^L \gets a^L - t/2$

$a^U \gets a^U + t/2$

$\lambda_1^L \gets 0$, $\lambda_1^U \gets M$
*[r]search for $\lambda_1$, $M$ is a large positive number

$\lambda_1^U - \lambda_1^L > \epsilon$
$\lambda_1^m \gets (\lambda_1^L + \lambda_1^U) / 2$

$h(\lambda_1^m) - \Delta < P_A - Q_A / \nu$

$\lambda_1^L \gets \lambda_1^m$

$\lambda_1^U \gets \lambda_1^m$

Following while-loop enlarges $\lambda_1^L$ and $\lambda_1^U$ until $[\lambda_1^L, \lambda_1^U]$ covers $\lambda_1^*$ such that $h(\lambda_1^*) = P_A-Q_A/\nu$ under numerical integration error
$(h(\lambda_1^L) + \Delta > P_A - Q_A/\nu)$  or  $(h(\lambda_1^U) - \Delta < P_A - Q_A/\nu)$
$t \gets \lambda_1^U - \lambda_1^L$

$\lambda_1^L \gets \lambda_1^L - t/2$

$\lambda_1^U \gets \lambda_1^U + t/2$

$(\lambda_1^L,\,(a^L - \lambda_1^U)/\nu)$ *[r]for soundness, choose the left endpoint of $\lambda_1$ and $\lambda_2$ range
DualBinarySearch for $\lambda_1$ and $\lambda_2$.

<Ref> is the dual variable search algorithm described in <Ref>.
From Line $1$ to $8$, we conduct binary search for quantity $\lambda_1 + \nu\lambda_2$;
from Line $14$ to $21$, we conduct binary search for quantity $\lambda_1$.
Notice that our binary search interval is initialized to be the non-negative interval.
This is because $Q(a^m,0) =0$ and $h(\lambda_1^m) = 0$ if $a^m$ and $\lambda_1^m$ are non-positive observed from <Ref>.

§.§ Guaranteeing Numerical Soundness

In computational method, to compute a practically sound robustness guarantee, we take the binary search error and numerical integration error into consideration.

Specifically, in <Ref>, we return a pair $(\lambda_1',\lambda_2')$ whose $R(\lambda_1',\lambda_2')$ lower bounds $R(\lambda_1, \lambda_2)$ where $(\lambda_1,\lambda_2)$ is the precise feasible pair.
We achieve so by returning $(\lambda_1^L,\,(a^L - \lambda_1^U)/\nu)$.
as we can see $a^L$ is an underestimation of actual $(\lambda_1 + \nu\lambda_2)$ in the presence of binary search error and numerical integration error.
We bound the numerical integration error by setting absolute error bound $\Delta = 1.5\times 10^{-8}$ in  function.
Then, $\lambda_1^L$ and $\lambda_1^U$ are underestimation and overestimation of the actual $\lambda_1$ in the presence of errors respectively.
As a result, $(a^L - \lambda_1^U)/\nu$ is an underestimation of $\lambda_2$.
Therefore, both $\lambda_1$ and $\lambda_2$ are underestimated and by the monotonicity of $R(\cdot,\cdot)$, the actual $R(\lambda_1,\lambda_2)$ would be underestimated to guarantee the soundness.

Then, since the evaluation of $R$ involves numerical integration, we compare the lower bound of $R$: $R(\lambda_1,\lambda_2) - \Delta$ in <Ref> (line 9) with $0.5$ to determine whether current robustness radius can be certified or not.

§.§ Complexity and Efficiency Analysis

In this appendix, we briefly analyze the computational complexity of computational method introduced in <Ref>.
Suppose the binary search precision is $\epsilon$, and each numerical integration costs $C$ time.
First, the search of certified robust radius costs $O(\log (\sqrt{d}/\epsilon))$.
For each searched radius, we first determine $P_A$ and $Q_A$ by running Neyman-Pearson-based certification, which has cost $O(\log(1/\epsilon)C)$.
Then, solving dual variables takes two binary search rounds, which has cost $O(\log(1/\epsilon)C)$.
The final one-time integration of $R(\lambda_1,\lambda_2)$ has cost $O(C)$.
Thus, overall time complexity is $O(\log (\sqrt{d}/\epsilon) \log(1/\epsilon) C)$, which is the same as classical Neyman-Pearson certification and grows slowly (in logarithmic factor) w.r.t. input dimension $d$.

In practice, the certification time is on average $\SI{5}{s}$ to $\SI{10}{s}$ per sample across different datasets.
For example, with $\sigma=0.50$ as the smoothing variance parameter, the certification time, as an overhead over Neyman-Pearson-based certification, is 10.53s, 4.53s, and 3.21s on MNIST, CIFAR-10, and ImageNet respectively.
This overhead is almost negligible compared with the sampling time for estimating $P_A$ and $Q_A$ which is around $\SI{200}{s}$ on ImageNet and is the shared cost of all randomized smoothing certification methods.
In summary, compared with standard Neyman-Pearson-based certification, the running time of is roughly the same.

§ EXTENSIONS AND PROOFS IN TIGHTNESS ANALYSIS

In this appendix, we provide formal proofs and theoretical extensions for the results in <Ref>.

§.§ Proof of Theorem 1

We let $D_c$ denote the decision region of $F_0$ for class $c$, i.e., $D_c := \{\vx:\,F_0(\vx) = c\}$.
Since $\gQ$ is supported on the decision region shifted by $\vx_0$, $f_0^\gQ(\vx_0)_c = 1$.
Thus, from $f_0^\gQ(\vx_0)_c$, we know $(\supp(\gQ) + \vx_0) \setminus S \subseteq D_c$, where $S$ is some set with zero measure under $\gQ + \vx_0$.
Since $0<q(\vx)/p(\vx)<+\infty$, $S$ also has zero measure under $\gP + \vx_0$.
On the other hand, by $0 < q(\vx)/p(\vx) < +\infty$, we can determine the probability mass of $\supp(\gQ)$ on $\gP$, i.e., $\Pr_{\vepsilon\sim\gP} [\vepsilon \in \supp(\gQ)]$.
Then, we observe that

\begin{align*}
f_0^\gP(\vx_0)_c =
% m_{\smooth}^M(\vx_0)_c =
\Pr_{\vepsilon\sim\gP} [F_0(\vx_0 + \vepsilon) = c] \\
= & \Pr_{\vepsilon\sim\gP} [\vepsilon \in \sR^d \setminus \supp(\gQ)]
\\
& \hspace{2em}
\cdot \Pr_{\vepsilon\sim\gP} [F_0(\vx_0+\vepsilon)=c \,|\, \vepsilon\in \sR^d\setminus \supp(\gQ)] \\
& \vspace{2em} + \Pr_{\vepsilon\sim\gP} [\vepsilon \in \supp(\gQ)] \cdot \Pr_{\vepsilon\sim\gP} [F_0(\vx_0 + \vepsilon)=c \,|\, \vepsilon \in \supp(\gQ)] \\
= & \Pr_{\vepsilon\sim\gP} [\vepsilon \in \sR^d \setminus \supp(\gQ)]
\\
& \hspace{2em}
\cdot \Pr_{\vepsilon\sim\gP} [F_0(\vx_0+\vepsilon)=c \,|\, \vepsilon\in \sR^d\setminus \supp(\gQ)] \\
&  + \Pr_{\vepsilon\sim\gP} [\vepsilon \in \supp(\gQ)].
\end{align*}

By the definition of $\supp(\gQ)$, we observe that $\Pr_{\vepsilon\sim\gP} [F_0(\vx_0+\vepsilon)=c \,|\, \vepsilon\in \sR^d\setminus \supp(\gQ)] = 0$.
As a result, we will find that $f_0^\gP(\vx_0)_c = \Pr_{\vepsilon\sim\gP} [\vepsilon \in \supp(\gQ)]$.
Then the certification method can know $\left(\left(\sR^d \setminus \supp(\gQ)\right) + \vx_0\right) \cap D_c$ has zero measure under $\gP + \vx_0$, i.e.,
In summary, the certification method can determine that $\supp(\gQ) + \vx_0$ differs from $D_c$ on some set $\Delta$ with zero measure under $\gP+\vx_0$.
Because $\gP$ has positive density everywhere, $\Delta$ also has zero measure under $\gP + \vx_0 + \vdelta$ for arbitrary $\vdelta\in\sR^d$.
Thus, for arbitrary $\vdelta\in\sR^d$, the certification method can compute out

\begin{align}
f_0^\gP(\vx_0+\vdelta)_c & = \Pr_{\vepsilon\sim\gP} [F_0(\vx_0+\vdelta+\vepsilon)=c] \\ & \hspace{-2em} = \Pr_{\vepsilon\sim\gP+\vdelta} [\vx_0 + \vepsilon\in D_c] = \Pr_{\vepsilon\sim\gP+\vdelta}[\supp(\gQ)]. \nonumber
\end{align}

Under the binary classification setting, it implies that for any $\vdelta\in\sR^d$, the $f_0^\gP(\vx_0+\vdelta)$ can be uniquely determined by certification method.
Since the smoothed classifier's decision at any $\vx_0+\vdelta$ is uniquely determiend by $f_0^\gP(\vx_0+\vdelta)$ (see <Ref>), the certification method can exactly know $\wF_0^\gP(\vx+\vdelta)$ for any $\vdelta$ and thus determine tightest possible certified robust radius $\rtight$.

§.§ Extending Theorem 1 to Multiclass Setting

For the multiclass setting, we define a variant of as follows.

Given $P_A \in [0,\,1]$ and $Q_A^\multi \in \sR^{C-1}$,

\begin{equation}
\small \label{eq:r-dsrs-multi}
\begin{aligned}
& \rDSRSmulti := \max r \quad \mathrm{s.t.} \\
& \forall F: \sR \to [C], f^\gP(\vx_0)_{y_0} = P_A, \\
& f^{\gQ_c}(\vx_0)_{c} = (Q_A^\multi)_c, c\in [C-1] \\
% f^\gQ(\vx_0)_{y_0} = Q_A \\
\Rightarrow & \forall \vx, \|\vx - \vx_0\|_p < r, \wF^\gP(\vx) = y_0.
\end{aligned}
\end{equation}

In the above definition, $\rDSRSmulti$ is the tightest possible certified radius with prediction probability $Q_A^\multi$, where each component of $Q_A^\multi$, namely $(Q_A^\multi)_c$, corresponds to the prediction probability for label $c$ under additional smoothing distribution $\gQ_c$.
Note that there are $(C-1)$ additional smoothing distributions $\{Q_c\}_{c=1}^{C-1}$ in this generalization for multiclass setting.

With this generation, the following corollary extends the tightness analysis in <Ref> from binary to the multiclass setting.

Suppose the original smoothing distribution $\gP$ has positive density everywhere, i.e., $p(\cdot) > 0$.
For multiclass classification with base classifier $F_0$, at point $\vx_0\in\sR^d$, for each class $c \in [C-1]$, let $\gQ_c$ be a distribution that satisfies: (1) its support is the decision region of $c$ shifted by $\vx_0$: $\supp(\gQ_c) = \{\vx - \vx_0: F_0(\vx) = c\}$; (2) for any $\vx\in\supp(\gQ)$, $0 < q_c(\vx)/p(\vx) < +\infty$.
Then, plugging $P_A = f_0^\gP(\vx_0)_c$ and $Q_A^\multi$ where $(Q_A^\multi)_c = f_0^{\gQ_c}(\vx_0)_c$ for $c\in [C-1]$ into <Ref>, we have $\rDSRSmulti = \rtight$ under any $\ell_p$ ($p \ge 1$).

Similar as the proof of <Ref>, the certification method can observe that for any $c \in [C-1]$, $f_0^{\gQ_c}(\vx_0)_c = 1$.
Thus, the method knows $\supp(\gQ_c) + \vx_0 \approx D_c$ for arbitrary $c \in [C-1]$.
Here, the “$\approx$” means that the difference between the two sets has zero measure under $\gP+\vx_0$.
Thus, for arbitrary $\vdelta\in\sR^d$ the certification method can precisely compute out $f^\gP(\vx_0)_c$ for any $c \in [C-1]$.
Since $f^\gP(\vx_0) \in \Delta^C$, we also know $f^\gP(\vx_0)_C$ and the smoothed classifier's prediction on $\vx_0+\vdelta$ can be uniquely determined.
Then, following the same argument in <Ref>'s proof, we can determine the tightest possible certified robust radius $\rDSRSmulti$.

To achieve the tightest possible certified radius $\rtight$, for binary classification, we only need one extra scalar as the additional information ($Q_A$), while for multiclass classification, we need $(C-1)$ extra scalars as the additional information ($Q_A^\multi \in \sR^{C-1}$).
Following the convention as discussed in <Ref>, we are interested in tight certification under the binary classification for sampling efficiency concerns.
Therefore, we focus on using only one extra scalar ($Q_A$) to additional information in .
In both <Ref> and <Ref>, we only need finite quantities to achieve tight certification for any smoothed classifier.
In contrast, other existing work requires infinite quantities to achieve such optimal tightness <cit.>.

§.§ Proof of Theorem 2

The proof of <Ref> is a bit complicated, which relies on several propositions and lemmas along with theoretical results in <Ref>.
At high level,
based on the standard Gaussian distribution's property (<Ref>), we find $Q_A = 1$ under concentration property (<Ref>).
With $Q_A = 1$, we derive a lower bound of $\rDSRS$ in <Ref>.
We then use: (1) the concentration of beta distribution $\Beta(\frac{d-1}{2},\frac{d-1}{2})$ (see <Ref>) for large $d$;
(2) the relative concentration of gamma $\Gamma(d/2,1)$ distribution around mean for large $d$ (see <Ref> and resulting <Ref>);
and (3) the misalignment of gamma distribution $\Gamma(d/2-k,1)$'s mean and median for small $(d/2-k)$ (see <Ref>) to lower bound the quantity in <Ref> and show it is large or equal to $0.5$.
Then, using the conclusion in <Ref> we conclude that $\rDSRS \ge 0.02\sigma\sqrt{d}$.

If random vector $\vepsilon \in \sR^d$ follows standard Gaussian distribution $\gN(\sigma)$,

\begin{equation}
\Pr [\|\vepsilon\|_2 \le T] =  \GammaCDF_{d/2}\left(\frac{T^2}{2\sigma^2}\right),
\end{equation}

where $\GammaCDF_{d/2}$ is the CDF of gamma distribution $\Gamma(d/2,1)$.

According to <cit.>, he volume of a $d$-dimensional ball, i.e., $d$-ball, with radius $r$ is
V_d(r) = \frac{\pi^{d/2}}{\Gamma(\frac{d}{2}+1)} r^d

\begin{equation}
\Vol(\{\vepsilon: \|\vepsilon\|_2=r\}) = V_d(r)' = \dfrac{d\pi^{d/2}}{\Gamma(\frac{d}{2}+1)} r^{d-1}.
\end{equation}

\begin{align*}
& \Pr[\|\vepsilon\|_2 \le T] \\
% = & \int_0^T \dfrac{1}{(2\pi\sigma^2)^{d/2}} \cdot \exp\left(-\frac{r^2}{2\sigma^2}\right) \cdot \Vol(\{\vepsilon: \|\vepsilon\|_2=T\}) \dif r \\
= & \int_0^T
\dfrac{1}{(2\pi\sigma^2)^{d/2}} \cdot \exp\left(-\frac{r^2}{2\sigma^2}\right) \cdot \mathrm{Vol}(\{\vepsilon: \|\vepsilon\|_2=r\}) \dif r \\
= & \int_0^T
\dfrac{1}{(2\pi\sigma^2)^{d/2}} \cdot \exp\left(-\frac{r^2}{2\sigma^2}\right) \cdot
\dfrac{d\pi^{d/2}}{\Gamma(\frac d 2 + 1)} r^{d-1} \dif r \\
= & \int_0^{T^2} \dfrac{1}{(2\sigma^2)^{d/2} \Gamma(\frac d 2)} \exp\left(-\frac{r}{2\sigma^2}\right) r^{d/2-1} \dif r \\
= & \dfrac{1}{\Gamma(\frac d 2)} \int_0^{\frac{T^2}{2\sigma^2}} \exp(-r) r^{d/2-1} \dif r = \GammaCDF_{d/2} \left( \frac{T^2}{2\sigma^2} \right).
\end{align*}

With <Ref>, now we can show that $Q_A = 1$ under the condition of <Ref> as stated in the following lemma.

Suppose $F_0$ satisfies $(\sigma,\Pcon)$-concentration property at input point $\vx_0 \in \sR^d$,
with additional smoothing distribution $\gQ = \gNgtrunc(k,T,\sigma)$ where $T^2 = 2\sigma^2\GammaCDF_{d/2}^{-1}(\Pcon)$ and $d/2 - 15 \le k < d/2$, we have

\begin{equation}
Q_A = \Pr_{\vepsilon\sim\gQ} [F_0(\vx_0+\vepsilon)=y_0] = 1.
\end{equation}

According to <Ref>, for $T'$ that satisfies

\begin{equation}
\Pr_{\vepsilon\sim\gN(\sigma)} [\|\vepsilon\|_2 \le T'] = \Pcon
\label{eq:pf-lem-thm-2-2-1}
\end{equation}

we have

\begin{equation}
\Pr_{\vepsilon\sim\gN(\sigma)} [F_0(\vx_0 + \vepsilon) = y_0 \,|\, \|\vepsilon\|_2 \le T'] = 1.
\label{eq:pf-lem-thm-2-2-2}
\end{equation}

With <Ref>,
from <Ref>, we have

\begin{equation}
T'^2 = 2\sigma^2 \GammaCDF_{d/2}^{-1}(\Pcon).
\end{equation}

Thus, <Ref> implies

\begin{equation}
\small
\Pr_{\vepsilon\sim\gN(\sigma)} [F_0(\vx_0 + \vepsilon) = y_0 \,|\, \|\vepsilon\|_2 \le \sqrt{2\sigma^2\GammaCDF_{d/2}^{-1}(\Pcon)}] = 1.
\end{equation}

Notice that $\gN(\sigma)$ has finite and positive density anywhere within $\{\vepsilon: \|\vepsilon\|_2 \le T'\}$.
Thus, $F_0(\vx_0+\vepsilon)=y_0$ for any $\vepsilon$ with $\|\vepsilon\|_2 \le T'$ unless a zero-measure set.

Now, we consider $\gQ$.
$\gQ = \gNgtrunc(k,T,\sigma)$ where $T = T'$, and $\gNtrunc$ has finite and positive density anywhere within $\{\vepsilon: \|\vepsilon\|_2 \le T'\} \setminus \{\mathbf{0}\}$.

\begin{align*}
Q_A = & \Pr_{\vepsilon\sim\gNgtrunc(k,T,\sigma)} [F_0(\vx_0 + \vepsilon) = y_0] \\
= & \Pr_{\vepsilon\sim \gNg(k,\sigma)} [F_0(\vx_0+\vepsilon)=y_0 \,|\, \|\vepsilon\|_2 \le T] \\
\overset{(*)}{=} & \Pr_{\vepsilon\sim \gNg(k,\sigma)}
[F_0(\vx_0+\vepsilon)=y_0 \,|\, \|\vepsilon\|_2 \le T, \vepsilon \neq \mathbf{0}] \\
= & 1.
\end{align*}

In the above equations, $(*)$ is because

\begin{align*}
& \Pr_{\vepsilon\sim\gNg(k,\sigma)} [\vepsilon = \mathbf{0} \,|\, \|\vepsilon\|_2 \le T] \\
\le & \lim_{r\to 0} \Pr_{\vepsilon\sim\gNg(k,\sigma)} [\|\vepsilon\|_2 \le r \,|\, \|\vepsilon\|_2 \le T] \\
= & \lim_{r\to 0}
C \int_0^r
x^{-2k} \exp\left(-\frac{x^2}{2\sigma'^2}\right) d x^{d-1} \dif x \\
= & C\lim_{r\to 0} \int_0^r d x^{d-2k-1} \exp\left(-\frac{x^2}{2\sigma'^2}\right) \dif x
= 0
\end{align*}

where $C$ is a constant, and the last equality is due to $d /2 > k \Rightarrow d - 2k - 1 \ge 0$.

With $Q_A = 1$, we can have a lower bound of $\rDSRS$ as stated in the following lemma.

Under the same condition as in <Ref>,
we let $\BetaCDF_{\frac{d-1}{2}}$ be the CDF of distributon $\Beta(\frac{d-1}{2},\frac{d-1}{2})$,

\begin{equation}
\begin{aligned}
& r_0 = \max u \\
\mathrm{s.t.} & \E_{t\sim\Gamma(\frac{d}{2}-k)}
\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma' \sqrt{2t} - u)^2}{4u\sigma' \sqrt{2t}} \right) \ge 0.5,
\end{aligned}
\label{eq:lem-thm-2-3-main}
\end{equation}

and let $\rDSRS$ be the tightest possible certified radius in under $\ell_2$ when smoothing distribution $\gP = \gNg(k,\sigma)$, then

\begin{equation}
\rDSRS \ge r_0.
\end{equation}

The proof shares the same core methodology as computational method introduced in <Ref>.
Basically, according to <Ref>, for any radius $r$,
let $\vdelta = (r,0,\dots,0)^\T$, if $\tC_{\vdelta}(P_A, Q_A) > 0.5$, then $\rDSRS \ge r$, where $Q_A = 1$ according to <Ref>, and by the $(\sigma,\Pcon)$-concentration property

\begin{equation}
P_A \ge \Pr_{\vepsilon\sim\gNg(k,\sigma)} [\|\vepsilon\|_2 \le T].
\end{equation}

Therefore, to prove the lemma, we only need to show that when

\begin{equation}
\E_{t\sim\Gamma(\frac{d}{2}-k)}
\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma' \sqrt{2t} - u)^2}{4u\sigma' \sqrt{2t}} \right) \ge 0.5,
\end{equation}

for any $\vdelta = (u,0,\dots,0)^\T$, $\tC_{\vdelta}(P_A, Q_A) > 0.5$.

By definition (<Ref>),

\begin{align*}
& \tC_{\vdelta} (P_A,Q_A) \\
= & \min_{f} \E_{\vepsilon\sim\gP} [f(\vepsilon+\vdelta)] \quad\mathrm{s.t.} \\
& \hspace{1em} \E_{\vepsilon\sim\gP} [f(\vepsilon)] = P_A,
\E_{\vepsilon\sim\gQ} [f(\vepsilon)] = Q_A, \\
& \hspace{1em}
0 \le f(\vepsilon) \le 1 \quad \forall \epsilon\in\sR^d \\
\ge & \min_{f} \E_{\vepsilon\sim\gP}
[f(\vepsilon+\vdelta)] \quad\mathrm{s.t.} \\
& \hspace{1em}
\E_{\vepsilon\sim\gQ} [f(\vepsilon)] = 1, \\
& \hspace{1em}
0 \le f(\vepsilon) \le 1 \quad \forall \epsilon\in\sR^d \\
= & \E_{\vepsilon\sim\gP} [f(\vepsilon+\vdelta)] \quad\mathrm{where}\quad  %\\
% & \hspace{1em}
f(\vepsilon) = \left\{
\begin{array}{lr}
1, & \|\vepsilon\|_2 \le T \\
0. & \|\vepsilon\|_2 > T
\end{array}
\right. \\
=: & V.
\end{align*}

We now compute $V$:

\begin{align*}
V = & \E_{\vepsilon\sim\gP} [\|\vepsilon + \vdelta\|_2 \le T] \\
= & \int_{\sR^d} p(\vx) \1[\|\vx + \vdelta\|_2 \le T] \dif \vx \\
\overset{(1)}{=} & \int_0^{\infty} y \dif y \int_{\substack{p(\vx) = y\\ \1[\|\vx + \vdelta\|_2 \le T}}
\dfrac{\dif \vx}{\|\nabla p(\vx)\|_2} \\
\overset{(2)}{=} & \int_0^\infty y\dif y \int_{\substack{p(\vx) = y\\ \1[\|\vx + \vdelta\|_2 \le T}}
- \dfrac{\dif \vx}{r_p'(r_p^{-1}(y))} \\
\overset{(3)}{=} & \int_0^\infty y \dif y
\dfrac{2\pi^{d/2}}{\Gamma(\frac d 2)} r_p^{-1}(y)^{d-1}
\cdot
\left(- \dfrac{1}{r_p'(r_p^{-1}(y))} \right)
\cdot \\
& \hspace{2em}
\Pr[\|\vx + \vdelta\|_2 \le T \,|\, p(\vx) = y] \\
\overset{(4)}{=} & \int_0^\infty
r_p(t)\dif t \dfrac{2\pi^{d/2}}{\Gamma(\frac d 2)} t^{d-1}
% \cdot \\
% & \hspace{2em}
\Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = t] \\
\overset{(5)}{=} & \int_0^\infty
\dfrac{1}{(2\sigma'^2)^{d/2-k}\pi^{d/2}}\cdot
\dfrac{\Gamma(d/2)}{\Gamma(d/2-k)}
t^{-2k} \exp\left(-\frac{t^2}{2\sigma'^2}\right) \cdot\\
& \hspace{2em} \dfrac{2\pi^{d/2}}{\Gamma(\frac d 2)} t^{d-1} \Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = t] \dif t \\
= & \dfrac{1}{\Gamma(\frac d 2 - k)} \int_0^\infty
t^{d/2-k-1} \exp(-t) \cdot \\
& \hspace{2em} \Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}] \dif t \\
= & \E_{t\sim\Gamma(\frac d 2 -k)} \Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}].
\end{align*}

In above equations, $(1)$ follows from level-set sliced integration extended from [Yang et al., 2020] and $p(\vx)$ is the density of distribution $\gP=\gNg(k,\sigma)$ at point $\vx$;
in $(2)$ we define $r_p(\|\vx\|_2) := p(\vx)$ noting that $\gP$ is $\ell_2$ symmetric and all $\vx$ with same $\ell_2$ length having the same $p(\vx)$, and we have $\|\nabla p(\vx)\|_2 = -r'_p(r_p^{-1}(y))$ since $y = r_p(\|\vx\|_2)$ and $r_p$ is monotonically decreasing;
$(3)$ uses

\begin{equation}
\begin{aligned}
\mathrm{Vol}(\{\vx: p(\vx)=y\}) & = \mathrm{Vol}(\{\vx: \|\vx\|_2 = r_p^{-1}(y)\}) \\
& = \dfrac{2\pi^{d/2}}{\Gamma(\frac d 2)} r_p^{-1}(y)^{d-1};
\end{aligned}
\end{equation}

$(4)$ changes the integration variable from $y$ to $t = r_p^{-1}(y)$;
and $(5)$ injects the concrete expression of $r_p$.

Now, we inspect $\Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}]$:
When $\|\vx\|_2 = \sigma'\sqrt{2t}$,
$\sum_{i=1}^d x_i^2 = 2t\sigma'^2$.
Meanwhile, $\|\vx + \vdelta\|_2 \le T$ means $(x_1+u)^2 + \sum_{i=2}^d x_i^2 \le T^2$.
Thus, when $\|\vx\|_2 = \sigma'\sqrt{2t}$,

\begin{equation}
\|\vx + \vdelta\|_2 \le T
\iff
\dfrac{x_1}{\sigma'\sqrt{2t}} \le \dfrac{T^2-u^2-2t\sigma'^2}{2u\sigma'\sqrt{2t}}.
\end{equation}

According to <cit.>, for $\vx$ uniformly sampled from sphere with radius $\sigma'\sqrt{2t}$, the component coordinate $\dfrac{1 + \frac{x_1}{\sigma'\sqrt{2t}}}{2} \sim \Beta(\frac{d-1}{2}, \frac{d-1}{2})$.

\begin{align}
& \Pr[\|\vx + \vdelta\|_2 \le T \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}] \nonumber \\
= & \BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{1 + \frac{T^2-u^2-2t\sigma'^2}{2u\sigma'\sqrt{2t}}}{2}
\right) \nonumber \\
= & \BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma'\sqrt{2t} - u)^2}{4u\sigma'\sqrt{2t}}
\right).
\end{align}

Finally, we get

\begin{equation}
V = \E_{t\sim\Gamma(\frac{d}{2}-k,1)}
\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma' \sqrt{2t} - u)^2}{4u\sigma' \sqrt{2t}} \right).
\label{eq:V}
\end{equation}

In other words,

\begin{equation}
\E_{t\sim\Gamma(\frac{d}{2}-k,1)}
\BetaCDF_{\frac{d-1}{2}} \left(
\dfrac{T^2 - (\sigma' \sqrt{2t} - u)^2}{4u\sigma' \sqrt{2t}} \right) \ge 0.5,
\end{equation}

we have $V \ge 0.5$, and thus $\tC_{\delta}(P_A,Q_A) \ge V \ge 0.5$, $\rDSRS \ge u$, which concludes the proof.

We require the following property of $\BetaCDF$.

There exists $d_0 \in \sN_+$, for any $d \ge d_0$,

\begin{equation}
\BetaCDF_{\frac{d-1}{2}}(0.6) \ge 0.999.
\end{equation}

We let $\rv \sim \Beta(\frac{d-1}{2},\frac{d-2}{2})$, then

\begin{align}
& \E [\rv] = 1/2, \\
& \Var [\rv] = \dfrac{(\frac{d-1}{2})^2}{(\frac{d-1}{2} + \frac{d-1}{2})^2(\frac{d-1}{2} + \frac{d-1}{2} + 1)} = \frac{1}{4d}.
\end{align}

Now, applying Chebyshev's inequality, we have

\begin{equation}
\Pr[|\rv-0.5| \ge 0.1] \le \dfrac{1}{0.04d}.
\end{equation}

\begin{align}
& \BetaCDF_{\frac{d-1}{2}}(0.6) \\
= & \Pr[\rv < 0.6] \\
= & 1 - \Pr[\rv \ge 0.6] \\
\ge & 1 - \Pr[|\rv - 0.5| \ge 0.1]
\ge 1 - \dfrac{1}{0.04d}.
\end{align}

Thus, when $d \ge 25000$, $\BetaCDF_{\frac{d-1}{2}}(0.6) \ge 0.999$.

We also require the following properties of the gamma distribution.

For any $\Pcon \in (0,1)$, there exists $d_0 \in \sN_+$, for any $d \ge d_0$,

\begin{equation}
\GammaCDF^{-1}_{d/2}(\Pcon) \ge 0.99 \cdot \dfrac{d}{2}.
\end{equation}

We let $\rv \sim \Gamma(d/2)$, then

\begin{align}
& \E[\rv] = d/2, \\
& \Var[\rv] = d/2.
\end{align}

We now apply Chebyshev's inequality and get

\begin{align*}
& \Pr[\rv < 0.99 \cdot d/2] \\
\le & \Pr[|\rv - d/2| > 0.01 \cdot d/2] \\
\le & \dfrac{20000}{d}.
\end{align*}

Thus, for any $\Pcon \in (0,1)$, when $d \ge \dfrac{20000}{\Pcon}$,

\begin{equation}
\GammaCDF_{d/2}\left(0.99\cdot \frac d 2\right) \le \dfrac{20000}{d} \le \Pcon,
\end{equation}

\GammaCDF^{-1}_{d/2}(\Pcon) \ge 0.99\cdot \frac d 2.

When $d/2 - 15 \le k < d/2$, $k, d \in \sN_+$,

\begin{equation}
\Pr_{t\sim\Gamma(d/2-k)} \left[t \le 0.98 \left(\frac d 2 - k\right)\right] \ge \frac{0.5}{0.999}.
\end{equation}

We prove the proposition by enumeration.
Notice that $d/2 - k \in \{0.5, 1.0, \cdots, 14.5, 15.0\}$, we enumerate $\GammaCDF_{d/2-k}(0.98(d/2-k))$ for each $(d/2-k)$ and get the following table.

$\frac{d}{2}-k$     $\GammaCDF_{d/2-k}(0.98(d/2-k))$
$\frac{d}{2}-k$     $\GammaCDF_{d/2-k}(0.98(d/2-k))$

0.5     0.6778
8.0     0.5245

1.0     0.6247
8.5     0.5224

1.5     0.5990
9.0     0.5204

2.0     0.5831
9.5     0.5186

2.5     0.5718
10.0     0.5168

3.0     0.5632
10.5     0.5152

3.5     0.5564
11.0     0.5136

4.0     0.5507
11.5     0.5121

4.5     0.5459
12.0     0.5107

5.0     0.5418
12.5     0.5093

5.5     0.5381
13.0     0.5080

6.0     0.5349
13.5     0.5068

6.5     0.5319
14.0     0.5056

7.0     0.5292
14.5     0.5044

7.5     0.5268
15.0     0.5033

On the other hand, $\frac{0.5}{0.999} \le 0.5001$, which concludes the proof.

Now we are ready to prove the main theorem.

According to <Ref>,
we only need to show that for $u = 0.02\sigma\sqrt d$,
<Ref> holds.
For sufficiently large $d$, indeed,

\begin{align*}
& \E_{t\sim\Gamma(d/2-k)} \BetaCDF_{(d-1)/2} \left(
\dfrac{T^2 - (\sigma'\sqrt{2t} - u)^2}{4u\sigma'\sqrt{2t}} \right) \\
& \overset{\text{\Cref{lem:thm-2-4}}}{\ge} 0.999 \E_{t\sim\Gamma(d/2-k)} \1\left[
\dfrac{T^2 - (\sigma'\sqrt{2t} - u)^2}{4u\sigma'\sqrt{2t}} \ge 0.6
\right] \\
& \overset{(*)}{\ge} 0.999 \E_{t\sim\Gamma(d/2-k)} \1\left[
t \le 0.98 \left(\frac d 2 - k\right)
\right] \\
& \overset{\text{\Cref{prop:thm-2-6}}}{\ge} 0.999 \cdot \dfrac{0.5}{0.999} = 0.5.
\end{align*}

Thus, from <Ref> we have $\rDSRS \ge u = 0.02\sqrt d$.

The inequality $(*)$ follows from <Ref>.

Under the condition of <Ref>,
for sufficiently large $d$,

\begin{equation}
t \le 0.98 \left(\frac d 2 - k\right) \Rightarrow
\dfrac{T^2 - (\sigma'\sqrt{2t} - u)^2}{4u\sigma'\sqrt{2t}} \ge 0.6.
\end{equation}

\begin{align*}
& \dfrac{T^2 - (\sigma'\sqrt{2t}-u)^2}{4u\sigma'\sqrt{2t}} \ge 0.6 \\
\iff & T^2 - (\sigma'\sqrt{2t} - u)^2 \ge 2.4u\sigma'\sqrt{2t} \\
\overset{\text{$x:=\sigma'\sqrt{2t}$}}{\iff} &
T^2 - (x-u)^2 \ge 2.4 ux \\
\iff & x^2 + 0.4ux + u^2 - T^2 \le 0.
\end{align*}

From <Ref>, we have

\begin{equation}
\begin{aligned}
& u^2 - T^2 \\
= & 0.0004 d\sigma^2 - 2\sigma^2 \GammaCDF^{-1}_{d/2}(\Pcon) \\
\le & 0.0004 d\sigma^2 - 0.99d\sigma^2 < 0.
\end{aligned}
\end{equation}

\begin{align*}
& x^2 + 0.4ux + u^2 - T^2 \le 0 \\
\iff & x \le \dfrac{-0.4u + \sqrt{0.16u^2 - 4(u^2-T^2)}}{2} \\
& \hspace{3em} = -0.2u + \sqrt{T^2 - 0.96u^2} \\
\iff & t \le \dfrac{(-0.2u + \sqrt{T^2 - 0.96u^2})^2}{2\sigma'^2}.
\end{align*}

Again, from <Ref>,

\begin{equation}
\begin{aligned}
& T^2 - 0.96u^2 \\
= & 2\sigma^2\GammaCDF_{d/2}^{-1}(\Pcon) - 0.96 \times 0.0004 d\sigma^2 \\
\ge & (0.99 - 0.96 \times 0.0004) d\sigma^2,
\end{aligned}
\end{equation}

and therefore

\begin{equation}
\begin{aligned}
& (-0.2u + \sqrt{T^2 - 0.96u^2})^2 \\
\ge & d\sigma^2 \left(-0.004 + \sqrt{0.99 - 0.96 \times 0.0004}\right)^2 \approx 0.9816 d\sigma^2  \\
\ge & 0.98 d\sigma^2.
\end{aligned}
\end{equation}

\begin{align*}
& t \le \dfrac{(-0.2u + \sqrt{T^2 - 0.96u^2})^2}{2\sigma'^2} \\
\Longleftarrow & t \le \dfrac{0.98d\sigma^2}{2\sigma^2} \cdot \dfrac{d-2k}{d} \\
\iff & t \le 0.98 \left( \frac d 2 - k \right).
\end{align*}

§.§ Theorem 6

Let $d$ be the input dimension and $F_0$ be the base classifier.
For an input point $\vx_0\in \sR^d$ with true class $y_0$, suppose $F_0$ satisfies $(\sigma,\Pcon)$-Concentration property and $\Pr_{\vepsilon \sim \gN(\sigma)} [F_0(\vx_0 + \vepsilon) = y_0] = \Pcon$ where $\Pcon < 1$.
The smoothed classifier $\wF_0^{\gP'}$ is constructed from $F_0$ and smoothed by generalized Gaussian $\gP'=\gNg(k_0,\sigma)$ where $k_0$ is a constant independent of input dimension $d$.
for any constant $c > 0$, there exists $d_0$, such that when input dimension $d \ge d_0$, any method cannot certify $\ell_2$ radius $c\sqrt{d}$,
where $T = \sigma\sqrt{2\GammaCDF_{d/2}^{-1}(\Pcon)}$ and $\GammaCDF_{d/2}$ is the CDF of gamma distribution $\Gamma(d/2,1)$.

We defer the proof to <Ref>.
This theorem suggests that, if we use generalized Gaussian whose $k$ is a constant with respect to input dimension $d$ or use standard Gaussian (whose $k=0$ is a constant) for smoothing, we cannot achieve $\Omega(\sqrt{d})$ certified radius rate from and any other certification method.

§.§ Proof of Theorem 3

The proof of <Ref> is based on three lemmas listed below.

Given $k_0 \in \sN$, for any $\epsilon > 0$, there exists $d_0$, such that when $d > d_0$,

\begin{equation}
\Pr_{\rt\sim\Gamma\left(\frac d 2 - k_0,1\right)} \left[\rt\le (1-\epsilon)\left(\frac d 2 - k_0\right)\right] \le \frac{0.48}{0.99}.
\label{eq:thm-3-pf-1}
\end{equation}

Given $\Pcon \ge 0$, for any $\epsilon > 0$, there exists $d_0$, such that when $d > d_0$,

\begin{equation}
T := \sigma\sqrt{2\GammaCDF^{-1}_{d/2}(\Pcon)} \le \sigma\sqrt{(1+\epsilon)d}.
\label{eq:thm-3-pf-2}
\end{equation}

For any $\epsilon > 0$, there exists $d_0$, such that when $d > d_0$,

\begin{equation}
\BetaCDF_{\frac{d-1}{2}} (0.5-\epsilon) \le 0.01.
\label{eq:thm-3-pf-3}
\end{equation}

Proofs of these lemmas are based on Chebyshev's inequality.

For $\rt\sim\Gamma(d/2 - k_0,1)$, we have

\begin{equation}
\E[\rt] = d/2 - k_0, \, \Var[\rt] = d/2 - k_0.
\end{equation}

By Chebyshev's inequality,

\begin{align}
% & \Pr \left[ \rt \le (1 - \epsilon)\left(\frac d 2 - k_0\right) \right] \nonumber \\
\hspace{-2em}\Pr \left[ \rt \le (1 - \epsilon)\left(\frac d 2 - k_0\right) \right] \le & \Pr \left[ |\rt - \E[\rt]| \ge \epsilon\left(\frac d 2 - k_0\right) \right] \nonumber \\
\le & \dfrac{1}{\epsilon^2 (\frac d 2 - k_0)}.
\end{align}

Picking $d_0 = 2\left( \frac{0.99}{0.48\epsilon^2} + k_0 \right)$ concludes the proof.

We define random variable $\rv \sim \Gamma(d/2,1)$, so $\E[\rv]=d/2, \Var[\rv]=d/2$.
By Chebyshev's inequality,

\begin{equation}
\begin{aligned}
& \Pr[\rv \le (1+\epsilon)d/2] \\
\ge & 1 - \Pr[\rv \ge (1+\epsilon)d/2] \\
\ge & 1 - \Pr[|\rv - \E[\rv]| \ge \epsilon d/2] \\
\ge & 1 - \dfrac{2}{d\epsilon^2}.
\end{aligned}
\end{equation}

Let $d_0 = \frac{2}{\epsilon^2 (1 - \Pcon)}$.
Thus, when $d > d_0$,
$\Pr[\rv \le (1+\epsilon)d/2] \ge 1 - \dfrac{2}{d_0\epsilon^2} = \Pcon$, which implies that $\GammaCDF_{d/2}((1+\epsilon)d/2) \ge \Pcon$ and $\GammaCDF_{d/2}^{-1}(\Pcon) \le (1+\epsilon)d/2$ and concludes the proof.

We define random variable $\rv\sim\Beta(\frac{d-1}{2},\frac{d-1}{2})$, and we have $\E[\rv] = 1/2, \Var[\rv] = \frac{1}{4d}$.
By Chebyshev's inequality,

\begin{equation}
\Pr[\rv\le 0.5-\epsilon] \le \Pr[|\rv-\E[\rv]|\ge \epsilon]
\le \dfrac{1}{4d\epsilon^2}.
\end{equation}

Let $d_0 = \frac{25}{d\epsilon^2}$, when $d > d_0$, $\Pr[\rv \le 0.5-\epsilon] \le 0.01$ and hence $\BetaCDF_{\frac{d-1}{2}} (0.5-\epsilon) \le 0.01$.

Now we are ready to prove the main theorem.

According to the above three lemmas, we pick $d_0$, such that when $d > d_0$, the followings hold simultaneously.

\begin{align}
& \Pr_{\rt\sim\Gamma\left(\frac d 2 - k_0,1\right)} \left[\rt\le \left(1-\frac{c^2}{8\sigma^2}\right)\left(\frac d 2 - k_0\right)\right] \le \frac{0.48}{0.99}, \label{eq:thm-3-pf-4} \\
& T = \sigma\sqrt{2\GammaCDF^{-1}_{d/2}(\Pcon)} \le \sigma\sqrt{\left(1+\frac{c^2}{8\sigma^2}\right)d}, \label{eq:thm-3-pf-5} \\
& \BetaCDF_{\frac{d-1}{2}} \left(0.5-\frac{c}{8\sigma}\right) \le 0.01. \label{eq:thm-3-pf-6}
\end{align}

We define vector $\vdelta = (c\sqrt d, 0, 0, \dots, 0)^\T$.
Since $F_0$ satisfies $(\sigma,\Pcon)$-concentration property and $\Pr_{\vepsilon\sim\gN(\sigma)} [F_0(\vx_0+\vepsilon)=y_0] = \Pcon$, up to a set of zero measure, the region $\{\vepsilon: F_0(\vx_0+\vepsilon)=y_0\}$ and region $\{\vepsilon: \|\vepsilon\|_2\le T\}$ coincide.

We now show that $\E_{\vepsilon\sim\gP'=\gNg(k_0,\sigma)}[F_0(\vx_0+\vdelta+\vepsilon)=y_0] < 0.5$ when $c \le \sigma \sqrt{8/7}$.

\begin{align}
& \E_{\vepsilon\sim\gP'} [F_0(\vx_0+\vdelta+\vepsilon)=y_0] \nonumber \\
= & \Pr_{\vepsilon\sim\gP'} [\|\vdelta + \vepsilon\|_2 \le T] \nonumber \\
= & \E_{\rt\sim\Gamma(\frac d 2 - k_0,1)} \BetaCDF_{\frac{d-1}{2}} \left( \dfrac{T^2 - (\sigma'\sqrt{2\rt} - c\sqrt{d})^2}{4\sigma'\sqrt{2\rt}\cdot c\sqrt{d}} \right) \nonumber \\
& \hspace{5em} \text{(from \Cref{eq:V})} \nonumber \\
\le & 0.99 \E_{\rt\sim\Gamma(\frac d 2 - k_0,1)} \1\left[ \dfrac{T^2 - (\sigma'\sqrt{2\rt} - c\sqrt{d})^2}{4\sigma'\sqrt{2\rt}\cdot c\sqrt{d}} \ge 0.5 - \frac{c}{8\sigma} \right] \nonumber \\
& + 0.01 \hspace{3em} \text{(by \Cref{eq:thm-3-pf-6} and $\BetaCDF_{\cdot}(\cdot) \le 1$)} \nonumber \\
= & 0.99 \Pr_{\rt\sim\Gamma(\frac d 2 - k_0,1)} \left[ \dfrac{T^2 - (\sigma'\sqrt{2\rt} - c\sqrt{d})^2}{4\sigma'\sqrt{2\rt}\cdot c\sqrt{d}} \ge 0.5 - \frac{c}{8\sigma} \right] \nonumber \\
& + 0.01. \label{eq:thm-3-pf-7}
\end{align}

\begin{align}
& \dfrac{T^2 - (\sigma'\sqrt{2\rt} - c\sqrt{d})^2}{4\sigma'\sqrt{2\rt}\cdot c\sqrt{d}} \ge 0.5 - \frac{c}{8\sigma} \nonumber \\
\iff & T^2 - 2\rt \sigma^2 \frac{d}{d-2k_0} - dc^2 + \frac{c^2d}{2}\sqrt{\frac{2\rt}{d-2k_0}} \ge 0 \\
\overset{\text{\Cref{eq:thm-3-pf-5}}}{\Longrightarrow} & \left(1+\frac{c^2}{8\sigma^2}\right)d\sigma^2 - 2\rt \sigma^2 \frac{d}{d-2k_0} - dc^2 \label{eq:thm-3-pf-8} \\
& \hspace{3em} + \frac{c^2d}{2}\sqrt{\frac{2\rt}{d-2k_0}} \ge 0.  \nonumber
\end{align}

We now inject $\rt=0$ and $\rt=\left(1-\frac{c^2}{8\sigma^2}\right)\left(\frac{d}{2} - k_0\right)$ to the LHS of <Ref>.

* When $\rt=0$,

\begin{align*}
\text{LHS of \Cref{eq:thm-3-pf-8}} & = \left(1+\frac{c^2}{8\sigma^2}\right) d\sigma^2 - dc^2 \\
& = d\left(\sigma^2 - \frac{7}{8}c^2\right) \ge 0.
\end{align*}

* When $\rt=\left(1-\frac{c^2}{8\sigma^2}\right)\left(\frac{d}{2} - k_0\right)$,

\begin{align*}
& \text{LHS of \Cref{eq:thm-3-pf-8}} \\
= & \left( 1 + \frac{c^2}{8\sigma^2} \right) d\sigma^2 - \frac{2d\sigma^2}{d-2k_0} \left(1-\frac{c^2}{8\sigma^2}\right) \left(\frac{d}{2} - k_0\right) \\
& - dc^2 + \frac{dc^2}{2} \sqrt{ 1 - \frac{c^2}{8\sigma^2} } \\
= & d\sigma^2 + \frac{dc^2}{8} - d\sigma^2 + \dfrac{dc^2}{8} - dc^2 + \frac{dc^2}{2} \sqrt{ 1 - \frac{c^2}{8\sigma^2} } \\
\le & \frac{dc^2}{4} - dc^2 + \dfrac{dc^2}{2} < 0.
\end{align*}

Notice that the LHS of <Ref> is a parabola with negative second-order coefficient.

\begin{equation}
\text{\Cref{eq:thm-3-pf-8}} \Longrightarrow
\rt \in \left[0, \left(1-\frac{c^2}{8\sigma^2}\right)\left(\frac{d}{2} - k_0\right)\right]
\end{equation}

and hence

\begin{equation}
\begin{aligned}
& \Pr_{\rt\sim\Gamma(\frac d 2 - k_0,1)} \left[ \dfrac{T^2 - (\sigma'\sqrt{2\rt} - c\sqrt{d})^2}{4\sigma'\sqrt{2\rt}\cdot c\sqrt{d}} \ge 0.5 - \frac{c}{8\sigma} \right] \\
\le & \Pr_{\rt\sim\Gamma(\frac d 2 - k_0,1)} \left[\rt \le  \left(1-\frac{c^2}{8\sigma^2}\right)\left(\frac{d}{2} - k_0\right) \right] \\
\le & \dfrac{0.48}{0.99}. \hspace{5em} \text{(by \Cref{eq:thm-3-pf-4})}
\end{aligned}
\end{equation}

Plugging this inequality to <Ref>, we get

\begin{equation}
\E_{\vepsilon\sim\gP'} [F_0(\vx_0+\vdelta+\vepsilon) = y_0] \le 0.99 \cdot \frac{0.48}{0.99} + 0.01 = 0.49.
\end{equation}

As a result, when $c \le \sigma\sqrt{8/7}$, the smoothed classifier $\wF_0^{\gP'}$ is not robust given the perturbation $\vdelta=(c\sqrt{d},0,0,\cdots,0)^\T$, since there may exist another $y'\neq y_0$ with $\E_{\vepsilon\sim\gP'} [F_0(\vx_0+\vdelta+\vepsilon) = y'] \ge 0.51$ so $\wF_0^{\gP'}(\vx_0+\vdelta) =y' \neq y_0$.

When $c > \sigma\sqrt{8/7}$, within the $\ell_2$ radius ball $c\sqrt{d}$, there exists perturbation vector $\vdelta=(c'\sqrt{d},0,0,\cdots,0)^\T$ fooling smoothed classifier $\wF_0^{\gP'}$ where $c'=\sigma\sqrt{8/7}$.
Hence, for any $c > 0$, there exists a perturbation within $\ell_2$ ball with radius $c\sqrt{d}$, such that smoothed classifier $\wF_0^{\gP'}$ can be fooled, and then any robustness certification method cannot certify $\ell_2$ radius $c\sqrt{d}$ since the smoothed classifier itself is not robust.

§ PROOFS OF COMPUTATIONAL METHOD

§.§ Proof of Strong Duality (Theorem 3)

We write down the Lagrangian dual function of  <Ref>:

\begin{equation}
\begin{aligned}
\Lambda(f,\,\lambda_1,\,\lambda_2)
:= & \E_{\vepsilon\sim\gP} [f(\vdelta + \vepsilon)]
- \lambda_1 \left( \E_{\vepsilon\sim\gP} [f(\vepsilon)] - P_A\right) \\
- \lambda_2 \left( \E_{\vepsilon\sim\gQ} [f(\vepsilon)] - Q_A\right).
\end{aligned}
\label{eq:new-proof-thm-3-lambda-def}
\end{equation}

Then, from $\tC$'s expression (<Ref>), we have

\begin{align}
& \tC_{\vdelta}(P_A, Q_A) \nonumber \\
= & \min_f \E_{\vepsilon\sim\gP} [f(\vepsilon + \vdelta)] \,\mathrm{s.t.}\, 0\le f(\vepsilon) \le 1 \, \forall \vepsilon\in\sR^d, \nonumber \\
& \hspace{3em}  \E_{\vepsilon\sim\gP} [f(\vepsilon)] = P_A, \E_{\vepsilon\sim\gQ} [f(\vepsilon)] = Q_A  \nonumber\\
= & \min_{f: \sR^d \to [0,1]} \max_{\lambda_1,\lambda_2\in \sR} \Lambda(f,\,\lambda_1,\,\lambda_2) \nonumber \\
\overset{(i)}{\ge} & \max_{\lambda_1,\lambda_2\in \sR} \min_{f: \sR^d \to [0,1]}
\Lambda(f,\,\lambda_1,\,\lambda_2) \nonumber \\
\overset{(ii)}{=} & \max_{\lambda_1,\lambda_2\in\sR}
\Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda_1 p(\vepsilon+\vdelta) + \lambda_2 q(\vepsilon+\vdelta)] \nonumber \\
& \hspace{2em} - \lambda_1 \Pr_{\vepsilon\sim\gP} [p(\vepsilon-\vdelta) < \lambda_1 p(\vepsilon) + \lambda_2 q(\vepsilon)] \nonumber  \\
& \hspace{2em} - \lambda_2 \Pr_{\vepsilon\sim\gQ} [p(\vepsilon-\vdelta) < \lambda_1 p(\vepsilon) + \lambda_2 q(\vepsilon)] \nonumber  \\
& \hspace{2em} + \lambda_1 P_A + \lambda_2 Q_A. \label{eq:new-proof-thm-3-ineq}
\end{align}

In the above equation, $(i)$ is from the min-max inequality.
For completeness, we provide the proof as such:
Define $g: \sR^2 \to \sR$ such that $g(\lambda_1,\lambda_2) := \min_{f:\sR^d \to [0,1]} \Lambda(f,\,\lambda_1,\lambda_2)$.
As a result, for any $\lambda_1,\lambda_2\in\sR$ and any $f:\sR^d \to [0,1]$,
$g(\lambda_1,\lambda_2) \le \Lambda(f,\lambda_1,\lambda_2)$.
So for any $f: \sR^d \to [0,1]$, $\max_{\lambda_1,\lambda_2\in\sR} g(\lambda_1,\lambda_2) \le \max_{\lambda_1,\lambda_2\in \sR} \Lambda(f,\lambda_1,\lambda_2)$, which implies

\begin{equation}
\max_{\lambda_1,\lambda_2\in\sR} g(\lambda_1,\lambda_2) \le \min_{f:\sR\to[0,1]} \max_{\lambda_1,\lambda_2\in\sR} \Lambda(f,\lambda_1,\lambda_2),
\label{eq:new-proof-thm-3-ineq-i}
\end{equation}

where LHS is the RHS of $(i)$ and RHS is the LHS of $(i)$.

In above equation, $(ii)$ comes from a closed-form solution of $f$ for $\Lambda(f,\lambda_1,\lambda_2)$ given $(\lambda_1,\lambda_2)\in \sR^2$.
Notice that we can rewrite $\Lambda(f,\lambda_1,\lambda_2)$ as an integral over $\sR^d$:

\begin{align}
& \Lambda(f,\,\lambda_1,\,\lambda_2) \nonumber \\
= & \E_{\vepsilon\sim\gP} [f(\vdelta+\vepsilon)] - \lambda_1 \E_{\vepsilon\sim\gP} [f(\vepsilon)] - \lambda_2 \E_{\vepsilon\sim\gQ} [f(\vepsilon)] \nonumber \\
& \hspace{2em}
+ \lambda_1 P_A + \lambda_2 Q_A \nonumber \\
= & \int_{\sR^d} f(\vx) \cdot \left( p(\vx - \vdelta) - \lambda_1 p(\vx) - \lambda_2 q(\vx) \right) \dif\vx \nonumber \\
& \hspace{2em}
+ \lambda_1 P_A + \lambda_2 Q_A.
\label{eq:new-proof-them-3-ineq-ii}
\end{align}

We would like to minimize over $f: \sR^d \to [0,\,1]$ in <Ref> and simple greedy solution reveals that we should choose

\begin{equation}
f(\vx) = \left\{
\begin{aligned}
1, \quad & p(\vx - \vdelta) - \lambda_1p(\vx) - \lambda_2q(\vx) < 0 \\
0. \quad & p(\vx - \vdelta) - \lambda_1p(\vx) - \lambda_2q(\vx) \ge 0
\end{aligned}
\right.
\label{eq:proof-thm-3-new-optimal-f}
\end{equation}

We inject this $f$ into <Ref> and get

\begin{equation}
\begin{aligned}
& \min_{f:\sR^d\to[0,1]} \Lambda(f,\,\lambda_1,\,\lambda_2) \\
= & \Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda_1 p(\vepsilon + \vdelta) + \lambda_2 q(\vepsilon + \vdelta)] \\
& \hspace{1em} + \lambda_1 \left( P_A - \Pr_{\vepsilon\sim\gP} [p(\vepsilon - \vdelta) < \lambda_1p(\vepsilon) + \lambda_2q(\vepsilon)] \right) \\
& \hspace{1em} + \lambda_2 \left( Q_A - \Pr_{\vepsilon\sim\gQ} [p(\vepsilon - \vdelta) < \lambda_1p(\vepsilon) + \lambda_2q(\vepsilon)] \right). \\
\end{aligned}
\label{eq:new-proof-thm-3-ineq-ii-2}
\end{equation}

Hence $(ii)$ holds.

On the other hand, we know that $\tD_\vdelta(P_A,Q_A)$ (defined by <Ref>) is feasible by theorem statement.
Denote $(\lambda_1^*, \lambda_2^*)\in \sR^2$ to a feasible solution to $\tD_\vdelta(P_A,Q_A)$ and $d^*$ to the objective value, then from the constraints of $(\tD)$ we know

\begin{equation}
\begin{aligned}
\Pr_{\vepsilon\sim\gP} [p(\vepsilon-\vdelta) < \lambda_1^* p(\vepsilon) + \lambda_2^* q(\vepsilon)] = P_A, \\
\Pr_{\vepsilon\sim\gQ} [p(\vepsilon-\vdelta) < \lambda_1^* p(\vepsilon) + \lambda_2^* q(\vepsilon)] = Q_A, \\
\Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda_1^* p(\vepsilon+\vdelta) + \lambda_2^* q(\vepsilon+\vdelta)] = d^*.
\end{aligned}
\label{eq:new-proof-thm-3-equalities}
\end{equation}

Plugging in these equalities into <Ref>, we have

\begin{equation}
\tC_\vdelta(P_A,Q_A)
\ge d^* - \lambda_1P_A - \lambda_2Q_A + \lambda_1P_A + \lambda_2 Q_A
= d^*.
\label{eq:new-proof-thm-3-oneside}
\end{equation}

At the same time, we define function $f^*: \sR^d \to [0,1]$ such that

\begin{equation}
f^*(\vx) = \1[p(\vx-\vdelta) - \lambda_1^* p(\vx) - \lambda_2^* q(\vx) < 0].
\end{equation}

From <Ref>, $f^*$ satisfies the constraints of ($\tC$) (<Ref>).
Since $(\tC)$ minimizes over all possible functions $f: \sR^d\to[0,1]$,
we have

\begin{equation}
\tC_\vdelta(P_A,Q_A) \le \E_{\vepsilon\sim\gP} [f^*(\vepsilon+\vdelta)] = d^*.
\label{eq:new-proof-thm-3-anotherside}
\end{equation}

Combining  <Ref>, we get $\tC_{\vdelta}(P_A,Q_A) = d^*$ and hence the strong duality holds.

§.§ Proofs of Results in Section 5.2

Suppose $f_1$ is the optimal solution to $\tC_\vdelta(P_A^1,\,Q_A^1)$ and $f_2$ is the optimal solution to $\tC_\vdelta(P_A^2,\,Q_A^2)$.
Due to the linearity of expectation, $(f_1 + f_2)/2$ satisfies all the constraints of <Ref> for $P_A = (P_A^1 + P_A^2) / 2$ and $Q_A = (Q_A^1 + Q_A^2) / 2$, i.e., $(f_1 + f_2)/2$ is feasible for $P_A = (P_A^1 + P_A^2) / 2$ and $Q_A = (Q_A^1 + Q_A^2) / 2$ with objective value $\left( \tC_\vdelta(P_A^1,\,Q_A^1) + \tC_\vdelta(P_A^2,\,Q_A^2) \right) / 2$.
Thus, we have

\begin{equation}
\begin{aligned}
& \tC_\vdelta\left( \frac{P_A^1 + P_A^2}{2},\, \frac{Q_A^1 + Q_A^2}{2} \right) \le \\
& \hspace{2em}
\dfrac{1}{2} \left( \tC_\vdelta(P_A^1,\,Q_A^1) + \tC_\vdelta(P_A^2,\,Q_A^2) \right)
\end{aligned}
\end{equation}

since $\tC$ is a minimization problem.
By definition, $\tC_\vdelta (P_A,\,Q_A)$ is convex.

Since $\tC_\vdelta(P_A,\,Q_A)$ is defined on a compact $\sR^2$ subspace, the convexity implies continuity.
The continuity property is used in the following proof of <Ref>.

Here, we only prove the monotonicity for functions $x \mapsto \min_y \tC_\vdelta (x,\,y)$ and $x \mapsto \argmin_{y} \tC_\vdelta (x,\,y)$.
The same statement for $y \mapsto \min_x \tC_\vdelta(x,\,y)$ and $y \mapsto \min_x \tC_\vdelta(x,\,y)$ is then straightforward due to the symmetry.

For simplification, we define $\tC'_\vdelta: x \mapsto \min_y \tC_\vdelta(x,\,y)$ and let $\tilde\tC_\vdelta: x\mapsto \argmin_y \tC_\vdelta(x,\,y)$.
We notice that both functions can be exactly mapped to the constrained optimization problem ($\tC'$) which removes the second constraint in <Ref> in ($\tC$):

\begin{align}
\underset{f}{\mathrm{minimize}} \quad & \E_{\vepsilon\sim\gP} [f(\vdelta + \vepsilon)] \\
\mathrm{s.t.} \quad & \E_{\vepsilon\sim\gP}[f(\vepsilon)] = x,\,\\
& 0 \le f(\vepsilon) \le 1 \quad \forall \vepsilon\sim\sR^d.
\end{align}

$\tC'_\vdelta(x)$ is the optimal objective to ($\tC'$) and $\bar\tC_\vdelta(x)$ is $\E_{\vepsilon\sim\gQ}[f^*(\vepsilon)]$ where $f^*$ is the optimal solution.

Either based on Neyman-Pearson lemma [[Neyman & Pearson, 1933]] or strong duality, $(\tC')$ is equivalent to $(\tD')$ defined as such:

\begin{align}
& \Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda p(\vepsilon + \vdelta)] \\
\mathrm{s.t.} \quad & \Pr_{\vepsilon\sim\gP} [p(\vepsilon - \vdelta) < \lambda p(\vepsilon)] = x. \label{eq:proof-prop-2-d-prime-cons}
\end{align}

For a given $x$, we only need to find $\lambda$ satisfying <Ref>.

\begin{align}
\tC'_\vdelta(x) & = \Pr_{\vepsilon\sim\gP} [p(\vepsilon) < \lambda p(\vepsilon + \vdelta)], \label{eq:proof-prop-2-alter-C-prime} \\
\bar\tC_\vdelta(x) & = \Pr_{\vepsilon\sim\gQ} [p(\vepsilon - \vdelta) < \lambda p(\vepsilon)]. \label{eq:proof-prop-2-alter-C-bar}
\end{align}

Now the monotonicity (what we would like to prove) is apparent.
For $x_1 < x_2$, from <Ref>, we have $\lambda_1 < \lambda_2$, since the probability density function $p$ is non-negative.
Thus, we inject $\lambda_1$ and $\lambda_2$ into <Ref> and <Ref>, and yield

\begin{equation}
\tC'_\vdelta(x_1) \le \tC'_\vdelta(x_2),\,
\bar \tC_\vdelta(x_1) \le \bar \tC_\vdelta(x_2),
\end{equation}

which concludes the proof.

We discuss the cases according to the branching statement in the algorithm (<Ref>).

If $\underline{q} > \underline{Q_A}$,

* if $\underline{q} \le \overline{Q_A}$, by definition we have $\tC_\vdelta(\underline{P_A},\,\underline{q}) \le \tC_\vdelta(\underline{P_A},\,y)$ for arbitrary $y$.
According to <Ref>, we also have $\tC_\vdelta(\underline{P_A},\,\underline{q}) \le \tC_\vdelta(x,\,y)$ for arbitrary $x \ge \underline{P_A}$ and arbitrary $y$.
Given that $(\underline{P_A},\,\underline{q}) \in [\underline{P_A},\,\overline{P_A}] \times [\underline{Q_A},\,\overline{Q_A}]$,
$(\underline{P_A},\,\underline{q})$ solves <Ref>;

* if $\underline{q} > \overline{Q_A}$, by convexity, $\tC_\vdelta(\underline{P_A},\,\overline{Q_A}) \le \tC_\vdelta(\underline{P_A},\,y)$ for $y \in [\underline{Q_A},\,\overline{Q_A}]$.

We further show that $\tC_\vdelta(\underline{P_A},\,\overline{Q_A}) \le \tC_\vdelta(x,\,\overline{Q_A})$ for $x \in [\underline{P_A},\,\overline{P_A}]$:
assume that this is not true, by <Ref>, the function $y \mapsto \argmin_{x} \tC_\vdelta(x,\,y)$ has function value larger than $\underline{P_A}$ at $y = \overline{Q_A}$.
Since $\tC_\vdelta(0,\,0) = 0$ is the global minimum of $\tC_\vdelta$, the function value at $y = 0$ is $x = 0$.
By <Ref>, there exists $y_0 \in [0,\,\overline{Q_A}]$ such that $\underline{P_A} = \argmin_x \tC_\vdelta(x,\,y_0)$.
Then, we get

\begin{equation}
\begin{aligned}
& \tC_\vdelta(\underline{P_A},\,y_0) \\
\overset{(i.)}{\le} & \tC_\vdelta(\argmin_x \tC_\vdelta(x,\,\overline{Q_A}),\,\overline{Q_A}) \\
\overset{(ii.)}{\le} & \tC_\vdelta(\underline{P_A},\,\overline{Q_A}),
\end{aligned}
\label{eq:proof-thm-4-1}
\end{equation}

where $(i.)$ follows from <Ref> for $y \mapsto \argmin_x \tC_\vdelta(x,\,y)$; $(ii.)$ is implied in the meaning of $\argmin_x \tC_\vdelta(x,\,\overline{Q_A})$.
Since $y_0 \in [0,\,\overline{Q_A}]$, <Ref> implies that $\underline{q}$ should be in $[0,\,\overline{Q_A}]$ as well, which violates the branching condition.
Thus, $\tC_\vdelta(\underline{P_A},\,\overline{Q_A}) \le \tC_\vdelta(x,\,\overline{Q_A})$ for $x\in [\underline{P_A},\,\overline{P_A}]$.

Using <Ref> for function $x \mapsto \argmin_y \tC_\vdelta(x,\,y)$ in interval $[\underline{P_A},\,\overline{P_A}]$ together with <Ref>, we get $\tC_\vdelta(x,\,\overline{Q_A}) \le \tC_\vdelta(x,\,y)$ for $x \in [\underline{P_A},\,\overline{P_A}]$ and $y \in [\underline{Q_A},\,\overline{Q_A}]$.
Thus, $(\underline{P_A},\,\overline{Q_A})$ solves <Ref>.

If $\underline{q} \le \underline{Q_A}$,

* if $\underline{p} \le \overline{P_A}$,
by definition we have $\tC_\vdelta(\max\{\underline{p},\underline{P_A}\},\,\underline{Q_A}) \le \tC_\vdelta(x,\,\underline{Q_A})$ for $x\in[\underline{P_A},\,\overline{P_A}]$.
According to <Ref> and condition $\underline{q} \le \underline{Q_A}$, we further have $\tC_\vdelta(\max\{\underline{p},\underline{P_A}\},\,\underline{Q_A}) \le \tC_\vdelta(\max\{\underline{p},\underline{P_A}\},y) \le \tC_\vdelta(x,\,y)$ for arbitrary $x \in [\underline{P_A}, \overline{P_A}]$ and $y \in [\underline{Q_A},\overline{Q_A}]$.
Given that $(\max\{\underline{p},\underline{P_A}\},\,\underline{Q_A}) \in [\underline{P_A},\,\overline{P_A}] \times [\underline{Q_A},\,\overline{Q_A}]$,
$(\underline{p},\,\underline{Q_A})$ solves <Ref>;

* if $\underline{p} > \overline{P_A}$,
according to <Ref>, $\tC_\vdelta(\overline{P_A},\,\underline{Q_A}) \le \tC_\vdelta(x,\,\underline{Q_A})$ for $x\in [\underline{P_A},\,\overline{P_A}]$.

We further show that $\tC_\vdelta(\overline{P_A},\, \underline{Q_A}) \le \tC_\vdelta(\overline{P_A},\, y)$ for $y \in [\underline{Q_A},\,\overline{Q_A}]$:
assume that this is not true, by <Ref>, the function $x \mapsto \argmin_y \tC_\vdelta(x,\,y)$ has function value larger than $\underline{Q_A}$ at $x = \overline{P_A}$.
Since $\tC_\vdelta(0,\,0)$ is the global minimum, by <Ref> on $x \mapsto \argmin_y \tC_\vdelta(x,\,y)$, there exists $x_0 \in [0,\,\overline{P_A}]$ such that $\underline{Q_A} = \argmin_y \tC_\vdelta(x_0,\,y)$.
Then, we get

\begin{equation}
\begin{aligned}
\tC_\vdelta(x_0,\,\underline{Q_A}) \\
\le & \tC_\vdelta(\overline{P_A},\,\argmin_y \tC_\vdelta(\overline{P_A},\,y)) \\
\le & \tC_\vdelta(\overline{P_A},\,\underline{Q_A})
\end{aligned}
\label{eq:proof-thm-4-2}
\end{equation}

following the similar deduction as in <Ref>.
Since $x_0 \in [0,\,\overline{P_A}]$, <Ref> implies that $\underline{p}$ should be in $[0,\,\overline{P_A}]$ as well, which violates the branching condition.
Thus, $\tC_\vdelta(\overline{P_A},\,\underline{Q_A}) \le \tC_\vdelta(\overline{P_A},\,y)$ for $y \in [\underline{Q_A},\,\overline{Q_A}]$.

Using <Ref> for function $y \mapsto \argmin_x \tC_\vdelta(x,\,y)$ in interval $[\underline{Q_A},\,\overline{Q_A}]$ together with <Ref>, we get $\tC_\vdelta(\overline{P_A},\,y) \le \tC_\vdelta(x,\,y)$ for $y\in [\underline{Q_A},\,\overline{Q_A}]$ and $x\in [\underline{P_A},\,\overline{P_A}]$.
Thus, $(\overline{P_A},\,\underline{Q_A})$ solves <Ref>.

§.§ Proof of Theorem 5

We first define $r_p(\|\vepsilon\|_2) = p(\vepsilon)$ and $r_q(\|\vepsilon\|_2) = q(\vepsilon)$, then easily seen the concrete expressions of $r_p$ and $r_q$ are:

\begin{align}
r_p(t) & = \dfrac{1}{(2\sigma'^2)^{d/2-k} \pi^{d/2}} \cdot \dfrac{\Gamma(d/2)}{\Gamma(d/2-k)}, \label{eq:pf-thm5-1} \\
r_q(t) & = \dfrac{\nu}{(2\sigma'^2)^{d/2-k} \pi^{d/2}} \cdot \dfrac{\Gamma(d/2)}{\Gamma(d/2-k)},
\label{eq:pf-thm5-2}
\end{align}

\begin{equation}
\nu := \dfrac{\Gamma(d/2-k)}{\gamma(d/2-k, \frac{T^2}{2\sigma'^2})} > 1
\end{equation}

and $\gamma$ is the lower incomplete Gamma function.

Now we use level-set integration similar as the proof in <Ref> to get the expressions of $P$, $Q$, and $R$ respectively.
Since $\gP$ and $\gQ$ are $\ell_2$-symmetric, without loss of generality, we let $\vdelta = (r,0,\dots,0)^\T$.

Suppose $P_T = r_p(T)$.

\begin{align*}
& P(\lambda_1,\lambda_2) \\
= & \Pr_{\vepsilon\sim\gP = \gNg(k,\sigma)} [p(\vepsilon-\vdelta) < \lambda_1 p(\vepsilon) + \lambda_2 q(\vepsilon)] \\
= & \int_{\sR^d} \1[p(\vx - \vdelta) < \lambda_1 p(\vx) + \lambda_2 q(\vx)] p(\vx) \dif \vx \\
= & \int_0^{P_T} y \dif y \int_{\substack{p(\vx)=y\\ p(\vx-\vdelta) < \lambda_1p(\vx)}} \dfrac{\dif \vx}{\|\nabla p(\vx)\|_2}
+ \\
& \hspace{1em}
\int_{P_T}^\infty y \dif y \int_{\substack{p(\vx)=y\\ p(\vx-\vdelta) < (\lambda_1 + \lambda_2\nu) p(\vx)}} \dfrac{\dif \vx}{\|\nabla p(\vx)\|_2} \\
= & \int_0^{P_T} y\dif y \dfrac{2\pi^{d/2}}{\Gamma(d/2)} r_p^{-1}(y)^{d-1} \left( - \dfrac{1}{r_p'(r_p^{-1}(y))} \right) \cdot \\
& \hspace{1em} \Pr[ p(\vx-\vdelta) \le \lambda_1 p(\vx) \,|\, p(\vx)=y ] + \\
& \hspace{1em} \int_{P_T}^{\infty} y\dif y \dfrac{2\pi^{d/2}}{\Gamma(d/2)} r_p^{-1}(y)^{d-1} \left( - \dfrac{1}{r_p'(r_p^{-1}(y))} \right) \cdot \\
& \hspace{1em} \Pr[p(\vx-\vdelta) < (\lambda_1+\lambda_2\nu) p(\vx) \,|\, p(\vx) = y] \\
\overset{y=r_p(t)}{=} &
\int_T^\infty r_p(t)\dif t \dfrac{2\pi^{d/2}}{\Gamma(d/2)} t^{d-1} \cdot \\
& \hspace{1em} \Pr[p(\vx-\vdelta) < \lambda_1p(\vx) \,|\, \|\vx\|_2 = t] + \\
& \hspace{1em}
\int_0^T r_p(t)\dif t \dfrac{2\pi^{d/2}}{\Gamma(d/2)} t^{d-1} \cdot \\
& \hspace{1em} \Pr[p(\vx-\vdelta) < (\lambda_1 + \lambda_2\nu) p(\vx) \,|\, \|\vx\|_2 = t] \\
= & \dfrac{1}{\Gamma(d/2-k)} \int_{T^2/(2\sigma'^2)}^\infty t^{d/2-k-1} \exp(-t) \dif t \cdot \\
& \hspace{1em} \Pr[p(\vx-\vdelta) < \lambda_1p(\vx) \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}] + \\
& \hspace{1em} \dfrac{1}{\Gamma(d/2-k)} \int_{0}^{T^2/(2\sigma'^2)} t^{d/2-k-1} \exp(-t) \dif t \cdot \\
& \hspace{1em} \Pr[p(\vx-\vdelta) < (\lambda_1 + \lambda_2\nu) p(\vx) \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}] \\
= & \E_{t\sim\Gamma(d/2-k,1)}
\left\{
\begin{array}{lr}
u_3(t,\lambda_1), & t \ge T^2 / (2\sigma'^2) \\
u_3(t,\lambda_1+\lambda_2\nu). & t < T^2 / (2\sigma'^2)
\end{array}
\right.
\end{align*}

Here, $u_3(t,\lambda) = \Pr[p(\vx-\vdelta) < \lambda p(\vx) \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}]$.

\begin{align*}
& Q(\lambda_1,\lambda_2) \\
= & \Pr_{\vepsilon\sim\gQ = \gNgtrunc(k,T,\sigma)} [p(\vepsilon-\vdelta) < \lambda_1p(\vepsilon) + \lambda_2q(\vepsilon)] \\
= & \int_{\|\vx\|_2\le T} \1[p(\vx-\vdelta) < \lambda_1p(\vx) + \lambda_2 q(\vx)] q(\vx) \dif \vx \\
= & \int_{\nu P_T}^\infty y \dif y \int_{\substack{q(\vx) = y\\ p(\vx-\vdelta) < (\lambda_1 + \lambda_2\nu) p(\vx)}} \dfrac{\dif \vx}{\|\nabla q(\vx)\|_2}  \\
\overset{y = r_q(t)}{=} & \int_0^T
r_q(t) \dif t \dfrac{2\pi^{d/2}}{\Gamma(d/2)} t^{d-1} \cdot \\
& \hspace{1em}
\Pr[p(\vx-\vdelta) < (\lambda_1 + \lambda_2 \nu) p(\vx) \,|\, \|\vx\|_2=t] \\
= & \nu \E_{t\sim \Gamma(d/2-k,1)} u_3(t,\lambda_1 + \lambda_2\nu) \cdot \1\left[ t \le \frac{T^2}{2\sigma'^2} \right].
\end{align*}

Now, for $R$:

\begin{align*}
& R(\lambda_1,\lambda_2) \\
= & \Pr_{\vepsilon\sim\gP=\gNg(k,\sigma)} [p(\vepsilon) < \lambda_1 p(\vepsilon + \vdelta) + \lambda_2 q(\vepsilon + \vdelta)] \\
= & \int_{\sR^d} \1[p(\vx) < \lambda_1 p(\vx + \vdelta) + \lambda_2 q(\vx + \vdelta)] p(\vx) \dif \vx \\
= & \E_{t\sim\Gamma(d/2-k,1)} u_4(t,\lambda_1,\lambda_2).
\end{align*}

Here, $u_4(t,\lambda_1,\lambda_2)=\Pr[\lambda_1 p(\vx + \vdelta) + \lambda_2 q(\vx + \vdelta) > r_p(\sigma'\sqrt{2t}) \,|\, \|\vx\|_2 = \sigma'\sqrt{2t}]$.

Plugging <Ref> into $P(\lambda_1,\lambda_2)$ and $Q(\lambda_1,\lambda_2)$, and then plugging <Ref> into $R(\lambda_1,\lambda_2)$, we yield the desired expressions in theorem statement.

Under the condition of <Ref>,
let $\vdelta = (r,0,\dots,0)^\T$,

\begin{equation}
\begin{aligned}