Additional experiments and resultsapp:results In the following we provide
additional results for our experiments with hierarchical randomized smoothing
for the tasks of node and image classification. For node classification we
experimentally observe that smoothed GAT models are more robust on average
(add:archs upper row). We therefore decide to show more results for smoothed
GAT models in the following. Recall that $r$ is the number of adversarial
nodes, $r_{d}$ is the number of adversarial deletions ($1\rightarrow 0$), and
$r_{a}$ the number of adversarial insertions ($0\rightarrow 1$). In add:1 we
show results for varying $r_{d}$ with fixed $r=1$ and $r_{a}=0$. In add:2 we
further vary $r_{a}$ with fixed $r=1$ and $r_{d}=0$. In add:3 we show results
for $r=1$, $r_{d}>0$ and $r_{a}>0$ and separately results for $r=2$. For image
classification, recall that $r$ is the number of pixels controlled by the
adversary and $\epsilon$ the perturbation strength in $\ell_{2}$-distance
across channels. In add:4 and add:5 we provide additional results for
different radii $r$ and $\epsilon$. Clearly, hierarchical randomized smoothing
is significantly expanding the Pareto-front w.r.t. the robustness-accuracy
trade-off in all of these settings. Since our certificates generalize ablation
and additive noise certificates into a new common framework, our method
performs either better or on-par with the existing methods under our threat
model. Note that in special threat models where e.g. adversaries have access
to all nodes ($r=N$), our certificates correspond to standard randomized
smoothing and consequently perform on-par with existing methods. In add:archs
(upper row) we compare different GNN architectures for the task of node
classification and generally observe that smoothed GAT models are more robust
on average. In add:archs (lower row) we provide further results for different
smoothing parameters when certifying the robustness of smoothed ResNet50
models on CIFAR10. Note the different effects on the certifiable robustness at
different radii $r$ when increasing $\sigma$ and $k$ (increasing $k$
corresponds to decreasing $p$). In add:6 we study the setting where
adversaries do not have access to target nodes in the graph but only to nodes
in the first- and second-hop neighborhoods (this setting has been previously
studied by [scholten2022randomized]). Such threat models allow us to introduce
skip-connections: We can forward the (unsmoothed) features of target nodes
separately through the GNN (without graph structure) and add the final
representations on the representations obtained with smoothing. In add:6 we
experimentally observe that such skip-connections can significantly expand the
Pareto-front w.r.t. robustness and accuracy of sparse and hierarchical
smoothing. In add:7 we perform experiments with and without the additional
additional column for smoothed GAT models on CoraML. Notably, training the
models on the extended matrices and using the additional column at inference
can result in better classifiers: As we experimentally observe for graph-
structured data, using the additional column significantly extends the Pareto-
front w.r.t. robustness and accuracy under specific threat models. In add:8 we
study the effect of different numbers of Monte-Carlo samples on the
certificate strength. We observe that for our setting using $\alpha=0.01$ and
$n_{1}=10{,}000$ is sufficient to compare our method to the baselines. Note
that in practice one would have to fix $\alpha$ and $n_{0},n_{1}$ prior to
computing certificates to ensure their statistical validity. Error bars. For
the node classification task we run every experiment 5 times with different
dataset splits and model initializations and report averaged results. The
average standard deviation is (3%, 3%, 4%, 4%) in clean accuracy and (3%, 3%,
2%, 3%) in certified accuracy (across all node classification experiments,
hierarchical smoothing only, sparse smoothing only, and ablation smoothing
only, respectively).

##  1 Full experimental setup (sec:exp)

We provide detailed information on the experimental setup as described in
sec:exp to ensure reproducibility of our results. Note that we also upload our
implementation.111Project page: https://www.cs.cit.tum.de/daml/hierarchical-
smoothing

###  Experimental setup for the task of node classification

Datasets and models. We implement smoothed classifiers for four architectures
with two message-passing layers: Graph convolutional networks (GCN)
kipf2017semi, graph attention networks (GAT and GATv2)
velivckovic2018graph,brody2021how, and APPNP klicpera2020predict. We train
models on the citation dataset Cora-ML
bojchevski2018deep,mccallum2000automating consisting of 2,810 nodes with 2,879
binary features, 7,981 edges and 7 classes. We follow the standard procedure
in graph machine learning shchur2018pitfalls and preprocess graphs into
undirected graphs containing only their largest connected component. Model
details. We implement all GNN models for two message-passing layers. We use 64
hidden channels for GCN kipf2017semi and APPNP klicpera2020predict, and 8
attention heads for 8 hidden channels for GAT and GATv2
velivckovic2018graph,brody2021how. For APPNP we further set the hyperparameter
k_hops to $10$ and the teleport probability to $0.15$. Node classification
task. For GNNs we follow the setup in scholten2022randomized. We evaluate our
certificates for the task of node classification in semi-supervised inductive
learning settings: We label 20 randomly drawn nodes per class for training and
validation, and 10% of the nodes for testing. We use the labelled training
nodes and all remaining unlabeled nodes as training graph, and insert the
validation and test nodes into the graph at validation and test time,
respectively. We train our models on training nodes, tune them on validation
nodes and compute certificates for all test nodes. Note that transductive
settings come with shortcomings when evaluating robustness certificates
scholten2022randomized, or in adversarial robustness on graph-structured data
in general gosch2023adversarial. Training details. We train models full-batch
using Adam (learning rate = 0.001, $\beta_{1}\\!=\\!0.9$,
$\beta_{2}\\!=\\!0.999$, $\epsilon\\!=\\!e{-08}$, weight decay = $5e{-04}$)
for a maximum of $1{,}000$ epochs with cross-entropy loss. Note that we
implement early stopping after 50 epochs. We also use a dropout of $0.5$ on
the hidden node representations after the first graph convolution for GCN and
GAT, and additionally a dropout of $0.5$ on the attention coefficients for GAT
and GATv2. Training-time smoothing parameters. During training we sample one
smoothed feature matrix (i.e. a matrix with (partially) added noise) in each
epoch and then we normally proceed training with the noised feature matrix.
Unless stated differently we also train our models on the higher-dimensional
matrices, i.e. we append the additional column to the matrices. Note that we
use the same smoothing parameters during training as for certification.
Experimental evaluation. We compare hierarchical smoothing to the additive
noise baseline by bojchevski2020efficient and the ablation baseline
scholten2022randomized as follows: We run $1{,}000$ experiments for each
method by (1) randomly drawing specific smoothing parameters ($p$ and
$p_{+},p_{-}$, depending on the method), (2) training $5$ smoothed classifiers
with different seeds for the datasplit and model initialization, (3) computing
robustness certificates under the corresponding smoothing distribution, and
(4) averaging resulting clean and certified accuracies. Exploring the Pareto-
front. In our exhaustive experiments we randomly explore the entire space of
possible parameter combinations. The reason for this is to demonstrate Pareto-
optimality, i.e. to find all ranges of smoothing parameters for which we can
offer both better robustness and accuracy. Note that in practical settings one
would have to conduct significantly less experiments to find suitable
parameters, and the task of efficiently finding the optimal parameters is a
research direction orthogonal to deriving robustness certificates. Smoothing
parameters for node classification. We randomly sample from the following
space of possible smoothing parameters (we draw from a uniform distribution
over these intervals):

  * [noitemsep]

  * •

Discrete hierarchical smoothing: $p\in[0.51,0.993]$, $p_{+}\in[0,0.05]$,
$p_{-}\in[0.5,1]$,

  * •

Sparse smoothing (additive noise baseline): $p_{+}\in[0,0.05]$,
$p_{-}\in[0.5,1]$

  * •

Interception smoothing (ablation smoothing baseline): $p\in[0.51,0.993]$

Note that the search space of hierarchical smoothing is the cross product of
the search spaces of the two baselines (sparse smoothing
bojchevski2020efficient and interception smoothing scholten2022randomized).
Notably, we still find that hierarchical smoothing is superior despite running
1,000 experiments equally for all three methods (see sec:exp). Note that we
run 1,000 experiments equally for all methods to establish an equal chance for
random outliers in accuracy. Also note that for clarity of the plots we only
include points with certified accuracies of larger than 0%.

###  Experimental setup for the task of image classification

Dataset and model architecture. For the task of image classification we follow
the experimental setup in levine2020robustness and train smoothed ResNet50
he2016resnet on CIFAR10 krizhevsky2009learning. CIFAR10 consists of 50,000
training images and 10,000 test images with 10 classes (6,000 images per
class). The colour images are of size 32x32 and have 3 channels. Note that for
runtime reasons we only certify a subset of the entire test set: For all
experiments we report clean and certified accuracy for the same subset
consisting of 1,000 randomly sampled test images. Training details. We train
the ResNet50 model with stochastic gradient descent (learning rate 0.01,
momentum 0.9, weight decay 5e-4) for 400 epochs using a batch-size of 128. We
also deploy a cosine learning rate scheduler loshchilov2017sgdr. At inference
we use a batch size of 300. We append the additional indicator channel (see
fig1) during training and certification. We train the smoothed classifiers
with the same noise level as during certification. Image preprocessing and
smoothing. We augment the training set with random horizontal flips and random
crops (with a padding of 4). We normalize all images using the channel-wise
dataset mean [0.4914, 0.4822, 0.4465] and standard deviation [0.2023, 0.1994,
0.2010]. For the ablation smoothing baseline we follow the implementation of
levine2020robustness and first normalize the images before appending three
more channels to the images. Specifically, for each channel $x$ we append the
channel $1-x$ (resulting in 6 channels in total). Finally, we set all channels
of ablated pixels to 0. We train and certify the extended images. See
levine2020robustness for a detailed setup of their method. For hierarchical
smoothing and the Gaussian smoothing baseline cohen2019certified we first add
the Gaussian noise and then normalize the images. Note that for hierarchical
smoothing we extend the images with the additional indicator channel (fig1),
resulting in 4 channels in total. Smoothing parameters for image
classification. Note that the ablation baseline levine2020robustness draws $k$
pixels to ablate from a uniform distribution $\mathcal{U}(k,N)$ over all
possible subsets with $k$ pixels of an image with $N$ total pixels. This is in
contrast to our approach, which selects pixels to add noise to with a
probability of $p$. For the comparison we choose the pixel selection
probability $p=1-\frac{k}{N}$ for every $k$, where $N$ is the number of pixels
in an image (i.e. we select $k$ pixels in expectation). This way, increasing
$k$ corresponds to decreasing $p$ and vice versa. We run experiments for every
10th $k$ from 0 to 500. Note that for $k$ larger that 512 we cannot obtain
certificates with either method since then $\Delta>0.5$. For Gaussian
smoothing we run experiments for every 100th $\sigma$ in the interval $[0,1]$,
and for hierarchical smoothing for every 20th $\sigma$ in $[0,1]$. Overall, we
run 50 experiments for the ablation baseline levine2020robustness, 100
experiments for the Gaussian noise baseline cohen2019certified, and 1,000
experiments for hierarchical smoothing.

###  Further experimental details

Certification parameters. All presented plots for graphs show results for the
tighter multi-class certificates, for images we compute the binary-class
certificates. Our certificates are probabilistic cohen2019certified, and we
use Monte-Carlo samples to estimate the smoothed classifiers with Clopper-
Pearson confidence intervals for $\alpha=0.01$ (thus our certificates hold
with high probability of $99\%$). We also apply Bonferroni correction to
ensure that the bounds hold simultaneously with confidence level $1-\alpha$.
We follow the procedure of cohen2019certified and, if not stated differently,
use $n_{0}=1{,}000$ Monte-Carlo samples for estimating the majority class at
test time and $n_{1}=10{,}000$ samples for the certification at test time.
Reproducibility. To ensure reproducibility we further use random seeds for all
random processes including for model training, for drawing smoothing
parameters and for drawing Monte-Carlo samples. We publish the source code
including reproducibility instructions and all random seeds. Computational
resources. All experiments were performed on a Xeon E5-2630 v4 CPU with a
NVIDIA GTX 1080TI GPU. In the following we report experiment runtimes (note
that for node classification we compute certificates for 5 different smoothed
models in each experiment): Regarding node classification experiments, the
average runtime of hierarchical smoothing for the $1{,}000$ experiments is
$1.6$ hours. For the additive noise baseline $1.5$ hours, and for the ablation
baseline $0.92$ hours. For image classification we only train a single
classifier. The overall training and certification process takes (on average)
$9.9$ hours for hierarchical smoothing, $8.6$ hours for Gaussian smoothing,
and $13.3$ hours for ablation smoothing. Third-party assets. For the sparse
smoothing baseline we use the implementation of
bojchevski2020efficient.222https://github.com/abojchevski/sparse_smoothing For
the ablation smoothing baseline for GNNs we use the implementation of
scholten2022randomized.333https://github.com/yascho/interception_smoothing We
implement all GNNs using PyTorch Geometric fey2019fast.444https://pytorch-
geometric.readthedocs.io The graph datasets are publicly available and can be
downloaded also e.g. using PyTorch Geometric. The CIFAR10 dataset is also
publicly available, for example via the torchvision
library.555https://pytorch.org/vision/stable/index.html

##  2 Robustness certificates for hierarchical smoothing (Proofs of
sec:certificates)

In the following we show the statements of sec:certificates about hierarchical
smoothing certificates and correctness of the following certification
procedure: Comment/* */ ruled KwDataInput KwResultReturn [H] $f$, $\bm{X}$,
$p$, Lower-level smoothing distribution $\mu$, Radius ($r$, $\epsilon$),
$n_{0}$, $n_{1}$, $1-\alpha$ counts0 $\leftarrow$ SampleUnderNoise($f$,
$\bm{X}$, $p$, $\mu$, $n_{0}$) counts1 $\leftarrow$ SampleUnderNoise($f$,
$\bm{X}$, $p$, $\mu$, $n_{1}$) $y_{A}$ $\leftarrow$ top index in counts0
$\underline{p_{A}}$ $\leftarrow$ ConfidenceBound(counts1[$y_{A}$], $n_{1}$,
$1-\alpha$) $\Delta\leftarrow 1-p^{r}$ $\underline{\tilde{p}_{A}}$
$\leftarrow$ LowerLevelWorstCaseLowerBound($\mu$, $\epsilon$,
$\frac{\underline{p_{A}}-\Delta}{1-\Delta}$)
$\underline{\tilde{p}_{A}}\cdot(1-\Delta)>\frac{1}{2}$ $y_{A}$ certified
ABSTAIN Binary-class Robustness Certificates for Hierarchical Randomized
Smoothing Here, LowerLevelWorstCaseLowerBound($\mu$,$\epsilon$,
$\underline{p_{A}}$) computes the existing certificate, specifically the
smallest $\underline{p_{\tilde{\bm{X}},y_{A}}}$ over the entire threat model,
e.g. $\Phi\left(\Phi^{-1}(\underline{p_{A}})-\frac{\epsilon}{\sigma}\right)$
for Gaussian smoothing. The remaining certification procedure is analogous to
the certification in cohen2019certified.

###  Proofs of sec:certificates

First, define the space
$\mathcal{Z}^{N\times(D+1)}\triangleq\mathcal{X}^{N\times
D}\times\\{0,1\\}^{N}$ and recall the definition of the reachable set
$\mathcal{A}\triangleq\mathcal{A}_{1}\cup\mathcal{A}_{2}$ with

$\mathcal{A}_{1}\triangleq\\{\bm{Z}\in\mathcal{Z}^{N\times(D+1)}\mid\forall
i:(\boldsymbol{\tau}_{i}=0)\Rightarrow(\bm{z}_{i}=\bm{x}_{i}|0)\\},$

and

$\mathcal{A}_{2}\triangleq\\{\bm{Z}\in\mathcal{Z}^{N\times(D+1)}\mid\forall
i:(\boldsymbol{\tau}_{i}=0)\Rightarrow(\bm{z}_{i}=\tilde{\bm{x}}_{i}|0)\\}.$

Intuitively, $\mathcal{A}$ contains only those $\bm{Z}$ that can be sampled
from either $\Psi_{\bm{X}}$ or $\Psi_{\tilde{\bm{X}}}$ (given the lower-level
smoothing distribution has infinite support). Note that these two sets are not
necessarily disjoint, so we cannot apply the Neyman-Pearson Lemma for discrete
random variables tocher1950extension on the upper-level smoothing distribution
$\phi$ yet. The following Proposition therefore provides a partitioning of the
space $\mathcal{Z}^{N\times(D+1)}$: 1prop:partitioning4 We want to show that
the four sets above (1) are pairwise disjoint, and (2) partition
$\mathcal{Z}^{N\times(D+1)}$, that is
$\mathcal{R}_{0}\uplus\mathcal{R}_{1}\uplus\mathcal{R}_{2}\uplus\mathcal{R}_{3}=\mathcal{Z}^{N\times(D+1)}$.
Clearly we have $\mathcal{Z}^{N\times(D+1)}=\mathcal{R}_{0}\uplus\mathcal{A}$
due to the definition of the region $\mathcal{R}_{0}$. Thus we only have to
show $\mathcal{R}_{1}\uplus\mathcal{R}_{2}\uplus\mathcal{R}_{3}=\mathcal{A}$:
(1) We show $\mathcal{R}_{1}$,$\mathcal{R}_{2}$,$\mathcal{R}_{3}$ are pairwise
disjoint:  We have $\mathcal{R}_{2}\cap\mathcal{R}_{1}=\emptyset$ and
$\mathcal{R}_{2}\cap\mathcal{R}_{3}=\emptyset$ since in region
$\mathcal{R}_{2}$ we have $\boldsymbol{\tau}_{i}=1$ for all $i\in\mathcal{C}$,
which does not hold for $\mathcal{R}_{2}$ or $\mathcal{R}_{3}$. Moreover,
$\mathcal{R}_{1}\cap\mathcal{R}_{3}=\emptyset$ because if
$\bm{Z}\in\mathcal{R}_{1}$ then $\bm{Z}\notin\mathcal{R}_{3}$ (and vice versa)
due the to definition of $\mathcal{A}$ and
$\mathcal{R}_{1}\subseteq\mathcal{A}_{1}$,
$\mathcal{R}_{3}\subseteq\mathcal{A}_{3}$. (2) We show
$\mathcal{R}_{1}\cup\mathcal{R}_{2}\cup\mathcal{R}_{3}=\mathcal{A}$:
“$\subseteq$”: Clearly
$\mathcal{R}_{1}\subseteq\mathcal{A}$,$\mathcal{R}_{2}\subseteq\mathcal{A}$,
and $\mathcal{R}_{3}\subseteq\mathcal{A}$ by definition of the regions.
“$\supseteq$”: Consider any
$\bm{Z}\in\mathcal{A}=\mathcal{A}_{1}\cup\mathcal{A}_{2}$. Either (1) all
perturbed rows are selected and $\bm{Z}\in\mathcal{R}_{2}$, or (2) at least
one perturbed row is not selected ($\exists i:\boldsymbol{\tau}_{i}=0$), in
which case we either have $\bm{Z}\in\mathcal{R}_{1}$ or
$\bm{Z}\in\mathcal{R}_{3}$, depending on whether $\bm{Z}\in\mathcal{A}_{1}$,
$\bm{Z}\in\mathcal{A}_{3}$, respectively. We further have to derive the
probability to sample $\bm{Z}$ of each region from both $\Psi_{\bm{X}}$ and
$\Psi_{\tilde{\bm{X}}}$ for applying the Neyman-Pearson Lemma
tocher1950extension for discrete random variables on the upper-level smoothing
distribution $\phi$ that selects rows before adding noise.
2prop:partitioningprobabilities2 The probability to select all $|\mathcal{C}|$
perturbed rows is $p^{r}$ since we sample each $\boldsymbol{\tau}_{i}$
independently from a Bernoulli: $\boldsymbol{\tau}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\text{Ber}(p)$. Consequently the probability to not select at
least one perturbed row is $1-p^{r}=\Delta$. The remaining follows from the
fact that matrices from region $\mathcal{R}_{1}$ cannot be sampled from
$\Psi_{\tilde{\bm{X}}}$ and matrices from $\mathcal{R}_{3}$ cannot be sampled
from $\Psi_{\bm{X}}$. Note that we need region $\mathcal{R}_{0}$ just for a
complete partitioning of $\mathcal{Z}^{N\times(D+1)}$. This region contains
essentially all matrices we cannot sample neither from $\Psi_{\bm{X}}$ nor
$\Psi_{\tilde{\bm{X}}}$. For example, matrices with a row that differs from
the same row in $\bm{X}$ and $\tilde{\bm{X}}$ without a previous selection of
that row cannot be sampled from the two distributions – simply because we
first select rows and then we only add noise to the rows that we selected.
Thus we have $\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[\bm{Z}\in\mathcal{R}_{0}]=0$ and
$\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[\bm{Z}\in\mathcal{R}_{0}]=0$. Before we
derive the lower bound on $p_{\bm{\tilde{X}},y}$ we first restate the Neyman-
Pearson Lemma as presented in sec:pre: 12 There are different views on
lemma:nplemma. In the context of statistical hypothesis testing, lemma:nplemma
states that the likelihood ratio test is the uniformly most powerful test when
testing the simple hypothesis $\mu_{\bm{X}}(\bm{W})$ against
$\mu_{\tilde{\bm{X}}}(\bm{W})$. For such views we refer to
[cohen2019certified]. Now we prove our main theorem: 1thm:mainthm1 We are
solving the following optimization problem:

$\underline{p_{\tilde{\bm{X}},y}}\quad\triangleq\quad\min_{h\in\mathbb{H}}\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y]\quad
s.t.\quad\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y]=p_{\bm{X},y}$

We can use the partitioning of prop:partitioning and apply the law of total
probability:

$\split\underline{p_{\bm{\tilde{X}},y}}=\min_{h\in\mathbb{H}}\quad&\sum_{i=0}^{3}\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{i}]\cdot\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[\bm{Z}\in\mathcal{R}_{i}]\\\
\text{s.t.}\quad\sum_{i=0}^{3}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{i}]\cdot\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[\bm{Z}\in\mathcal{R}_{i}]=p_{\bm{X},y}\\\
$

We apply prop:partitioningprobabilities and cancel out all zero probabilities:

$\split\underline{p_{\bm{\tilde{X}},y}}=\min_{h\in\mathbb{H}}\quad&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)+\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{3}]\cdot\Delta\\\
\text{s.t.\quad}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)+\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{1}]\cdot\Delta=p_{\bm{X},y}\\\
$

This is a minimization problem and to minimize we choose
$\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{3}]=0$
and
$\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{1}]=1$. In
other words, in the worst-case all matrices in $\mathcal{R}_{1}$ that we never
observe around $\tilde{\bm{X}}$ are correctly classified and the worst-case
classifier first uses the budget $p_{\bm{X},y}$ on region $\mathcal{R}_{1}$.
We obtain:

$\split\underline{p_{\bm{\tilde{X}},y}}=\min_{h\in\mathbb{H}}\quad&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)\\\
\text{s.t.\quad}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)=p_{\bm{X},y}-\Delta\\\
$

Here we can see that for $p=0$ (i.e. $\Delta=1$) we have just
$\underline{p_{\bm{\tilde{X}},y}}=0$. We can assume $p>0$ in the following. We
can further rearrange the terms:

$\split\underline{p_{\bm{\tilde{X}},y}}=(1-\Delta)\cdot\min_{h\in\mathbb{H}}\quad&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\\\
\text{s.t.\quad}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]=\frac{p_{\bm{X},y}-\Delta}{1-\Delta}\\\
$

Now we can apply lemma:nplemma:

$\split\underline{p_{\bm{\tilde{X}},y}}=&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[\bm{Z}\in
S_{\kappa}\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)\\\
\text{with}\kappa\in\mathbb{R}_{+}\text{s.t.}\\\
\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[\bm{Z}\in
S_{\kappa}\mid\bm{Z}\in\mathcal{R}_{2}]=\frac{p_{\bm{X},y}-\Delta}{(1-\Delta)}\\\
$ (1)

for
$S_{\kappa}\triangleq\left\\{\bm{Z}\in\mathcal{Z}^{N\times(D+1)}:\Psi_{\tilde{\bm{X}}}(\bm{Z})\leq\kappa\cdot\Psi_{\bm{X}}(\bm{Z})\right\\}$.
Note the difference between $S_{\kappa}$ and $\hat{S}_{\kappa}$. We need to
realize that due to independence we have:

$\aligned&\left\\{\bm{Z}\in\mathcal{R}_{2}:\Psi_{\tilde{\bm{X}}}(\bm{Z})\leq\kappa\cdot\Psi_{\bm{X}}(\bm{Z})\right\\}\\\
=\left\\{\bm{Z}\in\mathcal{R}_{2}:\Psi_{\tilde{\bm{X}}}(\bm{W},\boldsymbol{\tau})\leq\kappa\cdot\Psi_{\bm{X}}(\bm{W},\boldsymbol{\tau})\right\\}\\\
=\left\\{\bm{Z}\in\mathcal{R}_{2}:\mu_{\tilde{\bm{X}}}(\bm{W}|\boldsymbol{\tau})\phi(\boldsymbol{\tau})\leq\kappa\cdot\mu_{\bm{X}}(\bm{W}|\boldsymbol{\tau})\phi(\boldsymbol{\tau})\right\\}\\\
=\left\\{\bm{Z}\in\mathcal{R}_{2}:\prod_{i=1}^{N}\mu_{\tilde{\bm{x}}_{i}}(\bm{w}_{i}|\boldsymbol{\tau}_{i})\leq\kappa\cdot\prod_{i=1}^{N}\mu_{\bm{x}_{i}}(\bm{w}_{i}|\boldsymbol{\tau}_{i})\right\\}\\\
\stackrel{{\scriptstyle(*)}}{{=}}\left\\{\bm{Z}\in\mathcal{R}_{2}:\prod_{i\in\mathcal{C}}\mu_{\tilde{\bm{x}}_{i}}(\bm{w}_{i})\leq\kappa\cdot\prod_{i\in\mathcal{C}}\mu_{\bm{x}_{i}}(\bm{w}_{i})\right\\}\\\
=\left\\{\bm{Z}\in\mathcal{R}_{2}:\mu_{\tilde{\bm{X}}_{\mathcal{C}}}(\bm{W}_{\mathcal{C}})\leq\kappa\cdot\mu_{{\bm{X}}_{\mathcal{C}}}(\bm{W}_{\mathcal{C}})\right\\}\\\
$ (2)

where (*) is due to $\bm{x}_{i}=\tilde{\bm{x}}_{i}$ for all
$i\notin\mathcal{C}$ and $\boldsymbol{\tau}_{i}=1$ for $i\in\mathcal{C}$. In
other words this means that if $\bm{Z}\in S_{\kappa}\cap\mathcal{R}_{2}$ then
all $\bm{Z}^{\prime}\in\mathcal{R}_{2}$ that differ only in rows
$i\notin\mathcal{C}$ from $\bm{Z}$ are also in $S_{\kappa}\cap\mathcal{R}_{2}$
and integrate out (the “likelihood ratio” is only determined by rows
$i\in\mathcal{C}$). It follows that for fixed $\kappa\in\mathbb{R}_{+}$ we
have

$\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[\bm{Z}\in
S_{\kappa}\mid\bm{Z}\in\mathcal{R}_{2}]=\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]$

and

$\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[\bm{Z}\in
S_{\kappa}\mid\bm{Z}\in\mathcal{R}_{2}]=\Pr_{\bm{W}\sim\mu_{\bm{X}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]$

for
$\hat{S}_{\kappa}\triangleq\left\\{\bm{W}\in\mathcal{X}^{|\mathcal{C}|\times
D}:\mu_{\tilde{\bm{X}}_{\mathcal{C}}}(\bm{W})\leq\kappa\cdot\mu_{\bm{X}_{\mathcal{C}}}(\bm{W})\right\\}$.
All together we obtain:

$\underline{p_{\bm{\tilde{X}},y}}=\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]\cdot(1-\Delta)\quad\text{with}\kappa\in\mathbb{R}_{+}\text{s.t.}\quad\Pr_{\bm{W}\sim\mu_{\bm{X}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]=\frac{p_{\bm{X},y}-\Delta}{1-\Delta}$

##  3 Multi-class certificates for hierarchical randomized smoothing

In the main section of our paper we only derive so-called binary-class
certificates: If $\underline{p_{\tilde{\bm{X}},y^{*}}}>0.5$, then the smoothed
classifier $g$ is certifiably robust. We can derive a tighter certificate by
guaranteeing $\underline{p_{\tilde{\bm{X}},y^{*}}}>\max_{y\neq
y^{*}}\overline{p_{\tilde{\bm{X}},y}}$. To this end, we derive an upper bound
$p_{\tilde{\bm{X}},y}\leq\overline{p_{\tilde{\bm{X}},y}}$ in the following:

$\overline{p_{\tilde{\bm{X}},y}}\quad\triangleq\quad\max_{h\in\mathbb{H}}\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}}}[h(\bm{W})=y]\quad
s.t.\quad\Pr_{\bm{W}\sim\mu_{\bm{X}}}[h(\bm{W})=y]=p_{\bm{X},y}$ (3)

[cohen2019certified] show that the maximum of eq:maximize can be obtained by
using following lemma: [Neyman-Pearson upper bound] Given
$\bm{X},\tilde{\bm{X}}\in\mathcal{X}^{N\times D}$, distributions
$\mu_{\bm{X}},\mu_{\tilde{\bm{X}}}$, class label $y\in\mathcal{Y}$,
probability $p_{\bm{X},y}$ and the set
$S_{\kappa}\triangleq\left\\{\bm{W}\in\mathcal{X}^{N\times
D}:\mu_{\tilde{\bm{X}}}(\bm{W})\geq\kappa\cdot\mu_{\bm{X}}(\bm{W})\right\\}$
we have

$\overline{p_{\bm{\tilde{X}},y}}=\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}}}[\bm{W}\in
S_{\kappa}]\quad\text{with}\kappa\in\mathbb{R}_{+}\text{s.t.}\quad\Pr_{\bm{W}\sim\mu_{\bm{X}}}[\bm{W}\in
S_{\kappa}]=p_{\bm{X},y}$

See neyman1933ix,cohen2019certified. We can now use *lemma:nplemma2 to derive
an upper bound under hierarchical smoothing: [Neyman-Pearson upper bound for
hierarchical smoothing] Given fixed
$\bm{X},\tilde{\bm{X}}\in\mathcal{X}^{N\times D}$, let $\mu_{\bm{X}}$ denote a
smoothing distribution that operates independently on matrix elements, and
$\Psi_{\bm{X}}$ the induced hierarchical distribution over
$\mathcal{Z}^{N\times(D+1)}$. Given label $y\in\mathcal{Y}$ and the
probability $p_{\bm{X},y}$ to classify $\bm{X}$ as $y$ under $\Psi$. Define
$\hat{S}_{\kappa}\triangleq\left\\{\bm{W}\in\mathcal{X}^{r\times
D}:\mu_{\tilde{\bm{X}}_{\mathcal{C}}}(\bm{W})\geq\kappa\cdot\mu_{\bm{X}_{\mathcal{C}}}(\bm{W})\right\\}$.
We have:

$\overline{p_{\bm{\tilde{X}},y}}=\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]\cdot(1-\Delta)+\Delta\quad\text{with}\kappa\in\mathbb{R}_{+}\text{s.t.}\quad\Pr_{\bm{W}\sim\mu_{\bm{X}_{\mathcal{C}}}}[\bm{W}\in\hat{S}_{\kappa}]=\frac{p_{\bm{X},y}}{1-\Delta}$

where $\bm{X}_{\mathcal{C}}$ and $\tilde{\bm{X}}_{\mathcal{C}}$ denote those
rows $\bm{x}_{i}$ of $\bm{X}$ and $\tilde{\bm{x}}_{i}$ of $\tilde{\bm{X}}$
with $i\in\mathcal{C}$, that is $\bm{x}_{i}\neq\tilde{\bm{x}}_{i}$. We are
solving the following maximization problem:

$\overline{p_{\tilde{\bm{X}},y}}\quad\triangleq\quad\max_{h\in\mathbb{H}}\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y]\quad
s.t.\quad\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y]=p_{\bm{X},y}$

We can use the partitioning of prop:partitioning and apply the law of total
probability:

$\split\overline{p_{\bm{\tilde{X}},y}}=\max_{h\in\mathbb{H}}\quad&\sum_{i=0}^{3}\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{i}]\cdot\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[\bm{Z}\in\mathcal{R}_{i}]\\\
\text{s.t.}\quad\sum_{i=0}^{3}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{i}]\cdot\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[\bm{Z}\in\mathcal{R}_{i}]=p_{\bm{X},y}\\\
$

We apply prop:partitioningprobabilities and cancel out all zero probabilities:

$\split\overline{p_{\bm{\tilde{X}},y}}=\max_{h\in\mathbb{H}}\quad&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)+\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{3}]\cdot\Delta\\\
\text{s.t.\quad}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)+\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{1}]\cdot\Delta=p_{\bm{X},y}\\\
$

This is a maximization problem and to maximize we choose
$\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{3}]=1$
and
$\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{1}]=0$. In
other words, all matrices in $\mathcal{R}_{3}$ that we never observe around
${\bm{X}}$ are correctly classified and the worst-case classifier does not use
any budget $p_{\bm{X},y}$ on this region $\mathcal{R}_{3}$. We obtain:

$\split\overline{p_{\bm{\tilde{X}},y}}=\max_{h\in\mathbb{H}}\quad&\Pr_{\bm{Z}\sim\Psi_{\tilde{\bm{X}}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)+\Delta\\\
\text{s.t.\quad}\Pr_{\bm{Z}\sim\Psi_{\bm{X}}}[h(\bm{Z})=y\mid\bm{Z}\in\mathcal{R}_{2}]\cdot(1-\Delta)=p_{\bm{X},y}\\\
$

The remaining statement follows analogous to the proof in appA. In this
section we derived multi-class robustness certificates for randomized
smoothing and showed the correctness of the following certification procedure:
[H] $f$, $\bm{X}$, $p$, Lower-level smoothing distribution $\mu$, Radius ($r$,
$\epsilon$), $n_{0}$, $n_{1}$, $1-\alpha$ counts0 $\leftarrow$
SampleUnderNoise($f$, $\bm{X}$, $p$, $\mu$, $n_{0}$) counts1 $\leftarrow$
SampleUnderNoise($f$, $\bm{X}$, $p$, $\mu$, $n_{1}$) $y_{A}$, $y_{B}$
$\leftarrow$ top two indices in counts0 $\underline{p_{A}}$,
$\overline{p_{B}}$ $\leftarrow$
ConfidenceBounds(counts1[$y_{A}$],counts1[$y_{B}$], $n_{1}$, $1-\alpha$)
$\Delta\leftarrow 1-p^{r}$ $\overline{\tilde{p}_{B}}$ $\leftarrow$
LowerLevelWorstCaseUpperBound($\mu$, $\epsilon$,
$\frac{\underline{p_{B}}}{1-\Delta}$) $\underline{\tilde{p}_{A}}$ $\leftarrow$
LowerLevelWorstCaseLowerBound($\mu$, $\epsilon$,
$\frac{\underline{p_{A}}-\Delta}{1-\Delta}$)
$\underline{\tilde{p}_{A}}\cdot(1-\Delta)>\overline{\tilde{p}_{B}}\cdot(1-\Delta)+\Delta$
$y_{A}$ certified ABSTAIN Multi-class Robustness Certificates for Hierarchical
Randomized Smoothing Here, LowerLevelWorstCaseUpperBound($\mu$,$\epsilon$,
$\overline{p_{B}}$) computes the existing certificate, specifically the
largest $\overline{p_{\tilde{\bm{X}},y_{B}}}$ over the entire threat model,
e.g. $\Phi\left(\Phi^{-1}(\overline{p_{B}})+\frac{\epsilon}{\sigma}\right)$
for Gaussian smoothing. The remaining certification procedure is analogous to
the certification in cohen2019certified.

##  4 Hierarchical randomized smoothing using Gaussian isotropic smoothing
(Proofs of sec:initialization)

Having derived our general certification framework for hierarchical randomized
smoothing, we can instantiate it with specific smoothing distributions to
certify specific threat models. In the following we show the robustness
guarantees for initializing hierarchical randomized smoothing with Gaussian
isotropic smoothing, as already discussed in sec:initialization. We will also
derive the corresponding multi-class robustness certificates. 3prop:gaussian3
As outlined in subsection:background, [cohen2019certified] show that the
optimal value of eq:minimize under Gaussian smoothing is
$\Phi\left(\Phi^{-1}(p_{\bm{x},y})-\frac{||\bm{\delta}||_{2}}{\sigma}\right)$,
where $\Phi$ denotes the CDF of the normal distribution and
$\bm{\delta}=\bm{x}-\tilde{\bm{x}}$. Plugging this optimal value into
thm:mainthm yields for any
$\tilde{\bm{X}}\in\mathcal{B}_{2,\epsilon}^{r}(\bm{X})$:

$\underline{p_{\tilde{\bm{X}},y}}=\Phi\left(\Phi^{-1}\left(\frac{p_{\bm{X},y}-\Delta}{1-\Delta}\right)-\frac{||\bm{\delta}||_{2}}{\sigma}\right)(1-\Delta)$

where $\Delta=1-p^{|\mathcal{C}|}$ and
$\bm{\delta}=vec(\bm{X}-\tilde{\bm{X}})$. The remaining statements about the
minimum over the entire threat model follows from the fact that
$\underline{p_{\tilde{\bm{X}},y}}$ is monotonously decreasing in
$||\bm{\delta}||_{2}$, and monotonously decreasing in $r$ (see also prop:ball
in sec:certificates). Note that for fixed $\tilde{\bm{X}}$ we write
$\Delta=1-p^{|\mathcal{C}|}$ but for the entire threat model we have to use
the largest $\Delta=1-p^{r}$. 1cor:gaussian3 For the binary-class certificate,
the smoothed classifier $g$ is certifiably robust for fixed input $\bm{X}$ if
$\min_{\tilde{\bm{X}}\in\mathcal{B}_{2,\epsilon}^{r}(\bm{X})}\underline{p_{\tilde{\bm{X}},y^{*}}}>0.5$.
By using prop:gaussian we have

$\min_{\tilde{\bm{X}}\in\mathcal{B}_{2,\epsilon}^{r}(\bm{X})}\underline{p_{\tilde{\bm{X}},y^{*}}}>\frac{1}{2}$
$\Leftrightarrow\Phi\left(\Phi^{-1}\left(\frac{p_{\bm{X},y^{*}}-\Delta}{1-\Delta}\right)-\frac{\epsilon}{\sigma}\right)(1-\Delta)>\frac{1}{2}$
$\Leftrightarrow\Phi^{-1}\left(\frac{p_{\bm{X},y^{*}}-\Delta}{1-\Delta}\right)-\frac{\epsilon}{\sigma}>\Phi^{-1}\left(\frac{1}{2(1-\Delta)}\right)$
$\Leftrightarrow\epsilon<\sigma\left(\Phi^{-1}\left(\frac{p_{\bm{X},y^{*}}-\Delta}{1-\Delta}\right)-\Phi^{-1}\left(\frac{1}{2(1-\Delta)}\right)\right)$

In the following we additionally derive the corresponding multi-class
robustness certificates: Given continuous $\mathcal{X}\triangleq\mathbb{R}$,
consider the threat model $\mathcal{B}_{2,\epsilon}^{r}(\bm{X})$ for fixed
$\bm{X}\in\mathcal{X}^{N\times D}$. Initialize hierarchical smoothing with the
isotropic Gaussian smoothing distribution defined as
$\mu_{\bm{X}}(\bm{W})\triangleq\prod_{i=1}^{N}\mathcal{N}(\bm{w}_{i}|\bm{x}_{i},\sigma^{2}\bm{I})$.
Let $y\in\mathcal{Y}$ denote a label and $p_{\bm{X},y}$ the probability to
classify $\bm{X}$ as $y$ under $\Psi$. Define $\Delta\triangleq 1-p^{r}$. Then
we have:

$\max_{\tilde{\bm{X}}\in\mathcal{B}_{2,\epsilon}^{r}(\bm{X})}\overline{p_{\tilde{\bm{X}},y}}=\Phi\left(\Phi^{-1}\left(\frac{p_{\bm{X},y}}{1-\Delta}\right)+\frac{\epsilon}{\sigma}\right)(1-\Delta)+\Delta$

[cohen2019certified] show that the optimal value of eq:maximize under Gaussian
smoothing is
$\Phi\left(\Phi^{-1}(p_{\bm{x},y})+\frac{\epsilon}{\sigma}\right)$, where
$\Phi$ denotes the CDF of the normal distribution. The statement follows
analogously to prop:gaussian by directly plugging this optimal value into
nplemma:upper. [Multi-class certificates] Given continuous
$\mathcal{X}\triangleq\mathbb{R}$, consider the threat model
$\mathcal{B}_{2,\epsilon}^{r}(\bm{X})$ for fixed
$\bm{X}\in\mathcal{X}^{N\times D}$. Initialize the hierarchical smoothing
distribution $\Psi$ using the isotropic Gaussian distribution
$\mu_{\bm{X}}(\bm{W})\triangleq\prod_{i=1}^{N}\mathcal{N}(\bm{w}_{i}|\bm{x}_{i},\sigma^{2}\bm{I})$
that applies noise independently on each matrix element. Given majority class
$y^{*}\in\mathcal{Y}$ and the probability $p_{\bm{X},y^{*}}$ to classify
$\bm{X}$ as $y^{*}$ under $\Psi$. Then the smoothed classifier $g$ is
certifiably robust $g(\bm{X})=g(\tilde{\bm{X}})$ for any
$\tilde{\bm{X}}\in\mathcal{B}_{2,\epsilon}^{r}(\bm{X})$ if

$\Phi\left(\Phi^{-1}\left(\frac{p_{\bm{X},y^{*}}-\Delta}{1-\Delta}\right)-\frac{\epsilon}{\sigma}\right)>\max_{y\neq
y^{*}}\Phi\left(\Phi^{-1}\left(\frac{p_{\bm{X},y}}{1-\Delta}\right)+\frac{\epsilon}{\sigma}\right)+\frac{\Delta}{1-\Delta}$

where $\Phi^{-1}$ denotes the inverse CDF of the normal distribution and
$\Delta\triangleq 1-p^{r}$. For the tighter multi-class certificate we have to
guarantee $\underline{p_{\tilde{\bm{X}},y^{*}}}>\max_{y\neq
y^{*}}\overline{p_{\tilde{\bm{X}},y}}$. The statement follows directly from
prop:gaussian and prop:gaussian2.

##  5 Hierarchical smoothing for discrete lower-level smoothing distributions

Threat model details for discrete hierarchical smoothing experiments. In our
discrete experiments in sec:exp we model adversaries that perturb multiple
nodes in the graph by deleting $r_{d}$ ones (flip $1\rightarrow 0$) and adding
$r_{a}$ ones (flip $0\rightarrow 1$) in the feature matrix $\bm{X}$, that is
we consider the ball $\mathcal{B}_{r_{a},r_{d}}^{r}(\bm{X})$. Here we define
this threat model more formally:

$\aligned\mathcal{B}_{r_{a},r_{d}}^{r}(\bm{X})\triangleq\\{\tilde{\mathbf{X}}\in\mathcal{X}^{N\times
D}\mid&\sum_{i=1}^{N}\mathds{1}[{\mathbf{x}_{i}\neq\tilde{\mathbf{x}}_{i}}]\leq
r\\\ \sum_{i}\sum_{j}\mathds{1}(\tilde{\bm{X}}_{ij}=\bm{X}_{ij}+1)\leq
r_{a}\\\ \sum_{i}\sum_{j}\mathds{1}(\tilde{\bm{X}}_{ij}=\bm{X}_{ij}-1)\leq
r_{d}\\}$ (4)

where $\mathds{1}$ is an indicator function, $r$ the number of controlled
nodes, $r_{a}$ the number of adversarial insertions ($0\rightarrow 1$), and
$r_{d}$ the number of adversarial deletions ($1\rightarrow 0$). We expand the
result of thm1 to certificates for discrete lower-level smoothing
distributions $\mu$: 3thm23 This follows from applying the Neyman-Pearson-
Lemma for discrete random variables twice (see lemma:discreteNP). While
$\mathcal{R}_{1},\mathcal{R}_{2},\mathcal{R}_{3}$ represent “super regions”
for the upper-level smoothing distribution $\phi$,
$\mathcal{H}_{1},\ldots,\mathcal{H}_{I}$ further subdivide $\mathcal{R}_{2}$
into smaller regions of constant likelihood ratio, which is required for the
lower-level discrete smoothing distribution $\mu$. Notably, thm2 implies that
we can compute the LP of the original discrete certificate for the adjusted
probability of $(p_{\bm{X},y}-\Delta)/(1-\Delta)$, and then multiply the
resulting optimal value by the constant $1-\Delta$. In particular, the
feasibility of the Linear Program and of computing $\bm{v},\tilde{\bm{v}}$ can
be delegated to the certificate for the underlying $\mu$. Analogously for the
upper bound for discrete lower-level $\mu$ we need: [Discrete Neyman-Pearson
upper bound] Partition the input space $\mathcal{X}^{D}$ into disjoint regions
$\mathcal{R}_{1},\ldots,\mathcal{R}_{I}$ of constant likelihood ratio
$\mu_{\bm{x}}(\bm{w})/\mu_{\tilde{\bm{x}}}(\bm{w})=c_{i}\\!\in\\!\mathbb{R}\cup\\{\infty\\}$
for all $\bm{w}\in\mathcal{R}_{i}$. Define
$\tilde{\bm{r}}_{i}\triangleq\Pr_{\bm{w}\sim\mu_{\tilde{\bm{x}}}}(\bm{w}\in\mathcal{R}_{i})$
and
$\bm{r}_{i}\triangleq\Pr_{\bm{w}\sim\mu_{\bm{x}}}(\bm{w}\in\mathcal{R}_{i})$.
Then eq:maximize is equivalent to the linear program
$\overline{p_{\tilde{\bm{x}},y}}=\max_{\bm{h}\in[0,1]^{I}}\bm{h}^{T}\tilde{\bm{r}}\text{s.t.}\bm{h}^{T}\bm{r}=p_{\bm{x},y}$,
where $\bm{h}$ represents the worst-case classifier. Then we can show:
[Hierarchical upper bound for discrete lower-level smoothing distributions]
Given discrete $\mathcal{X}$ and a smoothing distribution $\mu_{\bm{x}}$ for
discrete data, assume that the corresponding discrete certificate partitions
the input space $\mathcal{X}^{rD}$ into disjoint regions
$\mathcal{H}_{1},\ldots,\mathcal{H}_{I}$ of constant likelihood ratio:
$\mu_{\bm{x}}(\bm{w})/\mu_{\tilde{\bm{x}}}(\bm{w})=c_{i}$ for all
$\bm{w}\in\mathcal{H}_{i}$. Then we can compute an upper bound
$\overline{p_{\bm{\tilde{X}},y}}\geq p_{\tilde{\bm{X}},y}$ by solving the
following Linear Program (LP):

$\overline{p_{\bm{\tilde{X}},y}}=\max_{\bm{h}}\quad\bm{h}^{T}\tilde{\bm{\nu}}\cdot(1-\Delta)+\Delta\quad\text{s.t.}\quad\bm{h}^{T}\bm{\nu}=\frac{p_{\bm{X},y}}{1-\Delta},\quad
0\leq\bm{h}_{i}\leq 1$

where
$\tilde{\bm{\nu}}_{q}\triangleq\Pr_{\bm{W}\sim\mu_{\tilde{\bm{X}}_{\mathcal{C}}}}[\text{vec}(\bm{W})\in\mathcal{H}_{q}]$
and
$\bm{\nu}_{q}\triangleq\Pr_{\bm{W}\sim\mu_{\bm{X}_{\mathcal{C}}}}[\text{vec}(\bm{W})\in\mathcal{H}_{q}]$
for all $q\in\\{1,\ldots,I\\}$. This follows analogously to thm2, specifically
by applying the Neyman-Pearson-Lemma upper bound for discrete random variables
(lemma:discreteNPUpper) twice (first using the “super regions” and then the
“lower-level” regions).

##  6 Hierarchical randomized smoothing using ablation smoothing

There are different types of certificates based on randomized ablation. One
way is to use ablation smoothing to certify threat models where adversaries
control multiple pixels in an image levine2020robustness. liu2021pointguard
ablate nodes in a point cloud by deleting them. bojchevski2020efficient
experiment with ablation smoothing for individual attributes in a graph and
show that their method is superior. scholten2022randomized directly mask out
entire nodes in a graph, yielding certificates against adversaries that
arbitrarily perturb node attributes. Ablation smoothing differences. Note that
the ablation baseline levine2020robustness draws $k$ pixels to ablate from a
uniform distribution $\mathcal{U}(k,N)$ over all possible subsets with $k$
pixels of an image with $N$ total pixels. levine2020robustness derive that
$\Delta^{L}=1-\frac{\binom{N-r}{k}}{\binom{N}{k}}$. This is in contrast to our
approach, which selects pixels to add noise to with a probability of $p$ for
which we have $\Delta^{S}=1-p^{r}$. In other words, we choose the pixel
selection probability $p=1-\frac{k}{N}$ for every $k$, where $N$ is the number
of pixels in an image (i.e. we select $k$ pixels in expectation). This way,
increasing $k$ corresponds to decreasing $p$ and vice versa. In short, for
small $N$ we have $\Delta^{S}<\Delta^{L}$. The differences are negligible for
our discussion, see scholten2022randomized for a detailed comparison. Ablation
and additive noise hybrid. Notably, our method generalizes ablation and
additive noise smoothing into a common framework: The additional column
indicates which entity will be ablated, but instead of ablating entities we
only add noise to them. If we ablate the selected rows instead of adding
noise, we can restore certificates for ablation smoothing as a special case of
our framework: 22 We can first apply thm2 and obtain:

$\underline{p_{\bm{\tilde{X}},y}}=min_{\bm{h}}\quad\bm{h}^{T}\tilde{\bm{\nu}}\cdot(1-\Delta)\quad\text{s.t.}\quad\bm{h}^{T}\bm{\nu}=\frac{p_{\bm{X},y}-\Delta}{1-\Delta},\quad
0\leq\bm{h}_{i}\leq 1$

for $\tilde{\bm{v}}=1$ and ${\bm{v}}=1$ (since we always ablate selected rows
there is only a single lower-level region $\mathcal{H}$). Solving this LP
yields

$\underline{p_{\bm{\tilde{X}},y}}=\frac{p_{\bm{X},y}-\Delta}{1-\Delta}(1-\Delta)=p_{\bm{X},y}-\Delta$

Note one can analogously compute the upper bound
$\overline{p_{\bm{\tilde{X}},y}}$ by applying the corresponding
thm:discreteHS2:

$\overline{p_{\bm{\tilde{X}},y}}=\frac{p_{\bm{X},y}}{1-\Delta}(1-\Delta)+\Delta=p_{\bm{X},y}+\Delta$

Altogether in this special case we obtain for the multi-class certificate
$\underline{p_{\tilde{\bm{X}},y^{*}}}>\max_{y\neq
y^{*}}\overline{p_{\tilde{\bm{X}},y}}$:

$p_{\bm{X},y^{*}}-\Delta>p_{\bm{X},\tilde{y}}+\Delta$

where $y^{*}$ is the most likely class and $\tilde{y}$ the follow-up class.
That is in this special case with ablation smoothing we again obtain the
multi-class guarantees derived by [levine2020robustness]. Related work. In the
context of randomized ablation (a special case of our framework),
[jia2021intrinsic] also define three regions and build an ablation certificate
using the Neyman-Pearson Lemma. Their regions are, however, different to ours,
and their certificates are designed for poisoning threat models (whereas we
consider evasion settings). Moreover, [liu2021pointguard] define three regions
similar to ours, but their certificate is only for point cloud classifiers and
suffers from the same limitations of ablation certificates
levine2020robustness that we overcome: Our key idea is to further process the
region $\mathcal{R}_{2}$ and to utilize the budget $p_{\bm{X},y}-\Delta$ by
plugging it into the certificate for the lower-level smoothing distribution,
which yields stronger robustness-accuracy trade-offs under our threat model.
Discussion. In general, the difference between ablation certificates and
additive noise certificates is that ablation smoothing does not evaluate the
base classifier on the entire input space. Specifically, the supported sample
spaces around $\bm{X}$ and $\tilde{\bm{X}}$ are overlapping but not the
identical. In this context, additive noise certificates are technically
tighter for threat models where adversarial perturbations are bounded. Since
our certificates generalize ablation smoothing and additive noise smoothing we
inherit advantages and disadvantages of both methods: While the certificate
for the lower-level smoothing distribution may be tight, the upper-level
smoothing distribution prevents us from evaluating the base classifier on the
entire matrix space. We leave the development of tight and efficient
certificates for hierarchical randomized smoothing to future work. Gray-box
certificates. scholten2022randomized derive gray-box certificates by
exploiting the underlying message-passing principles of graph neural networks.
In detail, they derive a tighter gray-box $\Delta$ that accounts for the fact
that nodes may be disconnected from a target node when randomly deleting edges
in a graph (and disconnected nodes do not affect the prediction for the
target). Interestingly, we can directly plug their gray-box $\Delta$ into our
framework to obtain gray-box robustness certificates for hierarchical
randomized smoothing.

##  7 Extended discussion

Motivation for the indicator variable as additional column. We certify
robustness of classifiers operating on the higher-dimensional matrix space,
which allows our certificates to be efficient and highly flexible. In the
following we motivate this by considering the scenario of certifying
robustness on the original data space $\mathcal{X}^{N\times D}$. The main
problem is that it is technically challenging to obtain an efficient
certificate that depends only on the radius (and not on fixed $\bm{X}$ and
$\tilde{\bm{X}}$): Consider a discrete lower-level smoothing distribution,
fixed $\bm{X},\tilde{\bm{X}}$, smoothed matrix $\bm{W}$ and row
$i\in\mathcal{C}$ such that $\bm{x}_{i}\neq\tilde{\bm{x}}_{i}$ but
$\bm{x}_{i}=\bm{w}_{i}$. To certify robustness on the original space we have
to consider the likelihood ratio for the marginal distribution
$\Psi_{\bm{X}}(\bm{W})\triangleq\sum_{\boldsymbol{\tau}}\Psi_{\bm{X}}(\bm{W},\boldsymbol{\tau})$:

$\frac{\Psi_{\bm{X}}(\bm{W})}{\Psi_{\tilde{\bm{X}}}(\bm{W})}=\frac{\sum_{\boldsymbol{\tau}}\Psi_{\bm{X}}(\bm{W},\boldsymbol{\tau})}{\sum_{\boldsymbol{\tau}}\Psi_{\tilde{\bm{X}}}(\bm{W},\boldsymbol{\tau})}.$

Let $\mathcal{C}$ describe all perturbed and $\tilde{\mathcal{C}}$ all clean
rows. Similar due to independence we have:

$\frac{\Psi_{\bm{X}}(\bm{W})}{\Psi_{\tilde{\bm{X}}}(\bm{W})}\overset{(1)}{=}\frac{\prod_{i=1}^{n}\Psi_{\mathbf{x}_{i}}(\mathbf{w}_{i})}{\prod_{i=1}^{n}\Psi_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})}=\frac{\prod_{i\in\mathcal{C}}\Psi_{\mathbf{x}_{i}}(\mathbf{w}_{i})\prod_{i\in\tilde{\mathcal{C}}}\Psi_{\mathbf{x}_{i}}(\mathbf{w}_{i})}{\prod_{i\in\mathcal{C}}\Psi_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})\prod_{i\in\tilde{\mathcal{C}}}\Psi_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})}\overset{(2)}{=}\frac{\prod_{i\in\mathcal{C}}\Psi_{\mathbf{x}_{i}}(\mathbf{w}_{i})}{\prod_{i\in\mathcal{C}}\Psi_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})}$

where $(1)$ is due independence and $(2)$ due to
$\bm{x}_{i}=\tilde{\bm{x}}_{i}$ for clean rows $i\in\tilde{\mathcal{C}}$. Now
we can consider the perturbed row $i\in\mathcal{C}$ that we chose above and
consider the likelihood ratio:

$\frac{\Psi_{\mathbf{x}_{i}}(\mathbf{w}_{i})}{\Psi_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})}=\frac{\sum_{\boldsymbol{\tau}_{i}}\mu_{\mathbf{x}_{i}}(\mathbf{w}_{i}|\boldsymbol{\tau}_{i})\phi(\boldsymbol{\tau}_{i})}{\sum_{\boldsymbol{\tau}_{i}}\mu_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i}|\boldsymbol{\tau}_{i})\phi(\boldsymbol{\tau}_{i})}\overset{(3)}{=}\frac{\mu_{\mathbf{x}_{i}}(\mathbf{w}_{i})p+(1-p)}{p\mu_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})}$

where (3) is due to $\bm{w}_{i}=\bm{x}_{i}$ but
$\bm{w}_{i}\neq\tilde{\bm{x}}_{i}$. Due to the mixture, there are two ways
sampling $\bm{w}_{i}$ from $\Psi_{\bm{x}_{i}}$: First, we may not “select” row
$i$ (and then we do not smooth the row). Second, we may “select” row $i$ but
not add any noise on the selected row. Due to the mixture, we cannot simply
cancel out dimensions in which the rows $\bm{x}_{i}$ and $\tilde{\bm{x}}_{i}$
agree as in bojchevski2020efficient, leading to the problem that the
likelihood ratio and subsequently the point-wise certificate depends on
$\tilde{\bm{X}}$. Consequently, certifying robustness against an entire ball
of perturbations is technically challenging since the likelihood ratio will
change for every $\tilde{\bm{X}}$. We can resolve this problem by realizing
that for reasonable smoothing parameters, the probability
$\mu_{\tilde{\bm{x}}_{i}}(\bm{w}_{i})=\mu_{\tilde{\bm{x}}_{i}}(\bm{x}_{i})$ is
close to zero (probability for flipping exactly and only the dimensions where
$\tilde{\bm{x}}_{i}$ disagrees with $\bm{x}_{i}$ is close to zero for small
radii and high dimensional data). Specifically, the likelihood ratio for such
elements is approaching $\infty$ (e.g. for sparse smoothing with fixed radius
$r_{a},r_{d}$ we have $\mu_{\tilde{\mathbf{x}}_{i}}(\mathbf{w}_{i})\rightarrow
0$ for $D\rightarrow\infty$). Note that if we choose a $\bm{W}$ and row $i$
such that $\bm{x}_{i}\neq\tilde{\bm{x}}_{i}$ but
$\tilde{\bm{x}}_{i}=\bm{w}_{i}$, we can derive the same but the other way
around and find the likelihood ratio approaching $0$. In other words, we
cannot avoid regions of $0$, and $\infty$ likelihood ratio for such smoothing
distributions where we can flip “entire rows” from $\bm{x}_{i}$ to
$\tilde{\bm{x}}_{i}$. Certifying robustness on the higher-dimensional matrix
space instead is efficient, highly flexible and tight for higher dimensional
rows as for example in graphs where rows represent node features. Efficiency
of hierarchical randomized smoothing. We discuss (1) the efficiency of the
smoothing process, and (2) the efficiency of the certification. First,
regarding the efficiency of the smoothing process, note that hierarchical
smoothing only adds noise on a subset of all entities. Technically this can be
more efficient, although drawing from a bernoulli before adding noise may
require more time in practice than directly drawing (e.g. gaussian) noise for
all matrix elements simultaneously. Second, regarding certification
efficiency, our method only requires the additional computation of the term
$\Delta$ which takes constant time. Yet, closed-form solutions for the largest
certifiable radius may not be available, depending on the lower-level
smoothing distribution. For example we may not be able to derive a closed-form
radius for the multi-class certificate under hierarchical Gaussian smoothing
(see app:gaussian). Still, the certification process can be implemented
efficiently even without closed-form solutions, e.g. via bisection or just by
incrementally increasing $\epsilon$ to find the largest certifiable radius.
Non-uniform entity-selection probability $p$. In general we can also use a
non-uniform $p$: Let $p_{i}$ denote a node-specific (non-uniform) selection
probability for node $i$ that the adversary does not control. To get
robustness certificates we have to make the worst-case assumption, i.e. nodes
with lowest selection probability will be attacked first by the adversary. For
the certificate this means we have to find the largest possible $\Delta$ (see
also Theorem 1). Specifically, we have to choose $\Delta\triangleq
1-\prod_{i\in\mathcal{I}}p_{i}$ where $\mathcal{I}$ contains the indices of
the $r$ smallest probabilities $p_{i}$. Then we can proceed as in our paper to
compute certificates for the smoothing distribution with non-uniform $p$. If
we additionally know that certain nodes will be attacked before other nodes,
we do not have to consider the first nodes with smallest $p_{i}$ but rather
the first nodes with smallest $p_{i}$ under an attack order. This can be
useful for example if one cluster of nodes is attacked before another, then we
can increase the selection probability for nodes in the first cluster and
obtain stronger certificates. Node-specific smoothing distributions. If some
nodes are exposed to different levels of attacks than others, we can obtain
stronger certificates by making the lower-level smoothing distribution node-
specific. This would allow us to add less noise for nodes that are less
exposed to attacks, potentially leading to even better robustness-accuracy
trade-offs. Note, however, that this increases threat model complexity since
we will have to treat all nodes separately. Input-dependent randomized
smoothing. We derive certificates for a smoothing distribution that adds
random noise independently on the matrix elements. This is the most studied
class of smoothing distributions and allows us to integrate a whole suite of
existing distributions. Although there are first approaches adding non-
independent noise, the corresponding certificates are still in their infancy
sukenik2022intriguing. Data-dependent smoothing is a research direction
orthogonal to ours, and we consider the integration of such smoothing
distributions as future work. Finding optimal smoothing parameters. The
randomized smoothing framework introduces parameters of the smoothing
distribution as hyperparameters (see e.g. cohen2019certified). In general,
increasing the probabilities to add random noise increases the robustness
guarantees but also decreases the accuracy. We introduce an additional
parameter allowing to better control the robustness-accuracy trade-off.
Specifically, larger $p$ adds more noise and increases the robustness but also
decreases accuracy. In our exhaustive experiments we randomly explore the
entire space of possible parameter combinations. The reason for this is to
demonstrate Pareto-optimality, i.e. to find all ranges for which we can offer
both better robustness and accuracy. In practice one would have to conduct
significantly less experiments to find suitable parameters. Note that the task
of efficiently finding the optimal parameters is a research direction
orthogonal to deriving robustness certificates. Limitation of fixed threat
models. Recent adversarial attacks (see e.g. kollovieh2023assessing) go beyond
the notion of $\ell_{p}$-norm balls and assess adversarial robustness using
generative methods. Our approach is limited as we can only certify robustness
against bounded adversaries that perturb at most $r$ objects with a
perturbation strength that is bounded by $\epsilon$ under some
$\ell_{p}$-norm. Since current robustness certificates can only provide
provable guarantees under fixed threat models, certifying robustness against
adversaries that are not bounded by $\ell_{p}$-norms is out of the scope of
this paper.