# Algorithmic Stability of Heavy-Tailed Stochastic Gradient Descent on Least
Squares

Anant Raj
Coordinated Science Laboraotry
University of Illinois Urbana-Champaign.
Inria, Ecole Normale Supérieure
PSL Research University, Paris, France.
<EMAIL_ADDRESS>Melih Barsbey
Department of Computer Engineering
Boğaziçi University, Istanbul, Turkey.
<EMAIL_ADDRESS>
Mert Gürbüzbalaban
Department of Management
Science and Information Systems
Rutgers University, Piscataway, USA.
<EMAIL_ADDRESS>
Lingjiong Zhu
Department of Mathematics
Florida State University, FL, USA.
<EMAIL_ADDRESS>
Umut Şimşekli
Inria, CNRS, Ecole Normale Supérieure
PSL Research University, Paris, France.
<EMAIL_ADDRESS>

###### Abstract

Recent studies have shown that heavy tails can emerge in stochastic
optimization and that the heaviness of the tails has links to the
generalization error. While these studies have shed light on interesting
aspects of the generalization behavior in modern settings, they relied on
strong topological and statistical regularity assumptions, which are hard to
verify in practice. Furthermore, it has been empirically illustrated that the
relation between heavy tails and generalization might not always be monotonic
in practice, contrary to the conclusions of existing theory. In this study, we
establish novel links between the tail behavior and generalization properties
of stochastic gradient descent (SGD), through the lens of algorithmic
stability. We consider a quadratic optimization problem and use a heavy-tailed
stochastic differential equation as a proxy for modeling the heavy-tailed
behavior emerging in SGD. We then prove uniform stability bounds, which reveal
the following outcomes: (i) Without making any exotic assumptions, we show
that SGD will not be stable if the stability is measured with the squared-loss
$x\mapsto x^{2}$, whereas it in turn becomes stable if the stability is
instead measured with a surrogate loss $x\mapsto|x|^{p}$ with some $p<2$. (ii)
Depending on the variance of the data, there exists a _‘threshold of heavy-
tailedness’_ such that the generalization error decreases as the tails become
heavier, as long as the tails are lighter than this threshold. This suggests
that the relation between heavy tails and generalization is not globally
monotonic. (iii) We prove matching lower-bounds on uniform stability, implying
that our bounds are tight in terms of the heaviness of the tails. We support
our theory with synthetic and real neural network experiments.

##  1 Introduction

Over the last decade, understanding the generalization behavior in modern
machine learning settings has been one of the main challenges in statistical
learning theory. Here, the main goal has been deriving upper-bounds on the so-
called _generalization error_ , i.e., the gap between the true and the
empirical risks $|F(\theta)-\hat{F}(\theta,X)|$, which are respectively
defined as follows:

$\displaystyle F(\theta):=\mathbb{E}_{x\sim
P_{X}}[f(\theta,x)],\qquad\hat{F}(\theta,X):=\frac{1}{n}\sum_{i=1}^{n}f(\theta,x_{i}),$
(1)

where $\theta\in\mathbb{R}^{d}$ denotes the _parameter vector_ ,
$f:\mathbb{R}^{d}\times\mathcal{X}\mapsto\mathbb{R}_{+}$ is the _loss
function_ , $\mathcal{X}$ is the space of _data points_ , $P_{X}$ is the
unknown _data distribution_ , and finally $X=\\{x_{1},\dots,x_{n}\\}$ denotes
a (random) dataset with $n$ points, such that each $x_{i}$ is independently
and identically distributed (i.i.d.) from $P_{X}$.

The past few years have witnessed the development of a variety of mathematical
frameworks for analyzing the generalization error (see e.g., [31, 21] for
recent surveys). In the context of empirical risk minimization (ERM), i.e.,
solving $\min_{\theta\in\mathbb{R}^{d}}\hat{F}(\theta,X)$, one promising
direction has been to explicitly take into account the statistical properties
of the _optimization algorithm_ used during training, which is typically
chosen as stochastic gradient descent (SGD) that is based on the following
recursion:

$\displaystyle\theta_{k+1}=\theta_{k}-\eta\nabla\tilde{F}_{k+1}(\theta_{k},X),$
(2)

where $\eta$ is the step-size (or learning-rate), and
$\nabla\tilde{F}_{k}(\theta,X):=\frac{1}{b}\sum_{i\in\Omega_{k}}f(\theta,x_{i})$
is the stochastic gradient, with $\Omega_{k}\subset\\{1,\dots,n\\}$ being a
random subset drawn with or without replacement, and $b:=|\Omega_{k}|\ll n$
being the batch-size. In this line of research, [45] and [32] empirically
demonstrated that, perhaps surprisingly, a _heavy-tailed_ behavior can emerge
in SGD in different ways, and the heaviness of the tails correlates with the
generalization error, suggesting that heavier tails indicate better
generalization.

Theoretically investigating these empirical observations, [19] and [22]
explored the origins of the observed heavy-tailed behavior. They
simultaneously showed that, in online SGD111The framework of [22] can handle
stochastic optimization algorithms other than SGD as well. (i.e., when the
data is streaming), due the multiplicative nature of the gradient noise, i.e.,
$\nabla\tilde{F}_{k}(\theta,X)-\nabla\hat{F}(\theta,X)$, the distribution of
the iterates $\theta_{k}$ can converge to a heavy-tailed distribution as
$k\to\infty$. Furthermore, [19] showed that, when the loss $f$ is a quadratic
and the data distribution is Gaussian, the tails become monotonically heavier
when $\eta$ gets larger or $b$ gets smaller.

Due to the fact that analyzing the heavy-tailed behavior arising from (2) can
be highly non-trivial, relatively simpler heavy-tailed mathematical models
have been used as a proxy for the original heavy-tailed SGD recursion in
stationarity; e.g., SGD with heavy-tailed noise, i.e.,

$\displaystyle\theta_{k+1}=\theta_{k}-\eta_{k+1}\left[\nabla\hat{F}(\theta_{k},X)+\xi_{k+1}\right],\quad\text{with}\quad\mathbb{E}[\|\xi_{k}\|^{2}]=+\infty,\quad\forall
k=1,2,\dots,$ (3)

where $\xi_{k}$ denotes the heavy-tailed noise and $\eta_{k}$ denotes a
sequence of decreasing step-sizes. It has been revealed that another
interesting situation emerges in this setting, this time in the behavior of
the optimization error. Notably, [53, Remark 1] pointed out that, when the
loss function $f$ is chosen as a simple quadratic, i.e.,
$f(\theta,x)=\|\theta\|^{2}$, we have that
$\mathbb{E}[\|\theta_{k}-\theta_{\star}\|^{2}]=\mathbb{E}[\|\theta_{k}\|^{2}]=\mathbb{E}[f(\theta_{k},x)]=+\infty$
for all $k$, where $\theta_{\star}=0$ is the global minimum of $f$. While this
result might appear daunting as it might seemingly suggest that “SGD diverges”
under heavy-tailed perturbations, [50] refined this result and showed that, if
there exists $p\in[1,2)$ such that $\mathbb{E}[\|\xi_{k}\|^{p}]<\infty$, then
$\mathbb{E}[\|\theta_{k}-\theta_{\star}\|^{p}]$ converges to zero, for a class
of strongly convex losses $f$. This result is particularly remarkable, since
it shows that, even when the iterates may diverge under the ‘true’ loss
function $f$ (which SGD is originally trying to minimize), i.e.,
$\mathbb{E}[f(\theta_{k},x)]=+\infty$, they might still converge to the
_minimum of the original loss_ $\theta_{\star}$ when a surrogate loss function
$\tilde{f}$ is used for measuring the optimization error, which in this
example is $\tilde{f}(\theta,x)=\|\theta\|^{p}$ with $p<2$.

In an initial attempt for formalizing the relation between the tail behavior
and generalization, [46] also modeled the original heavy-tailed recursion (2)
by using a proxy and considered the following stochastic differential equation
(SDE) as a model (which can be seen as a continuous-time version of (3)):

$\displaystyle\mathrm{d}\theta_{t}=-\nabla\hat{F}(\theta_{t},X)\mathrm{d}t+\Sigma(\theta_{t})\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(4)

where $\Sigma:\mathbb{R}^{d}\mapsto\mathbb{R}^{d\times d}$ is a matrix-valued
function and $\mathrm{L}^{\alpha}_{t}$ denotes a heavy-tailed $\alpha$-stable
Lévy process, which is a random process parameterized by $\alpha\in(0,2]$,
such that a smaller $\alpha$ indicates heavier tails (we will make the
definition of $\mathrm{L}^{\alpha}_{t}$ precise in the next section). They
showed that, under several assumptions on the SDE (4), the worst-case
generalization error over the trajectory, i.e.,
$\sup_{t\in[0,1]}|\hat{F}(\theta_{t},X)-F(\theta_{t})|$ scales with the
intrinsic dimension of the trajectory $(\theta_{t})_{t\in[0,1]}$, which is
then upper-bounded as a particular function of the tail-exponent around a
local minimum, indicating that heavier tails imply lower generalization error.
Their results were later extended to discrete-time recursions as well in [23].
More recently, [3] linked heavy-tails to generalization through a notion of
compressibility in the over-parameterized regimes. Yet, these bounds require
several topological and statistical regularity assumptions that are hard to
verify in realistic settings, and the experiments in [3] illustrated that the
relation between the tail-exponent and the generalization error is not always
monotonic; hence, a generalization bound that requires less assumptions while
being more pertinent to the practical observations is still missing.

In this study, we aim at establishing novel links between tail behavior and
generalization and address the aforementioned shortcomings. We consider the
problem through the lens of _algorithmic stability_ [7, 20], and explore the
effects of heavy tails on the stability of SGD. Similar to recent work [1, 19,
27], in order to have a more explicit control over the problem, we limit our
scope to quadratic optimization, and consider the following SDE as a proxy for
heavy-tailed SGD:

$\mathrm{d}\theta_{t}=-\frac{1}{n}(X^{\top}X)\theta_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(5)

where $\Sigma\in\mathbb{R}^{d\times d}$ is a matrix that scales the noise, and
by a slight abuse of notation we represent the dataset as a matrix
$X\in\mathbb{R}^{n\times d}$, such that $i$-th row of $X$ is equal to $x_{i}$.
This SDE naturally arises from the ERM problem with the loss function being
$f(\theta,x)=(\theta^{\top}x)^{2}$. As the learning algorithm, we consider the
case where we assume that we have a sample from the stationary distribution of
(5) and analyze the stability of this sample. Our contributions are as
follows:

  * •

As opposed to classical SDEs driven by a Brownian motion (rather than
$\alpha$-stable Lévy processes $\mathrm{L}^{\alpha}_{t}$ as we consider here),
the stationary distribution of (5) does not admit a simple analytical closed-
form expression. As a remedy, we perform the stability analysis in the Fourier
domain, and introduce new proof techniques.

  * •

We prove upper-bounds on the stability of (5), which suggest that the
algorithm will not be stable, when $\alpha<2$ and the stability is measured
with respect to the quadratic loss $(\theta^{\top}x)^{2}$. We further show
that, when the stability is instead measured with respect to a surrogate loss
function $|\theta^{\top}x|^{p}$ with $p<\alpha<2$, the algorithm in turn
becomes stable, where the level of stability depends on $\alpha$, among
several other quantities. This result reveals a similar phenomenon to that of
[53] and [50] as discussed above. Furthermore, our results do not require any
non-trivial assumptions, compared to the existing heavy-tailed generalization
bounds [46, 23, 3].

  * •

Our theory further discloses an interesting property: depending on the
variance of the data distribution $P_{X}$, there exits an $\alpha_{0}>1$, such
that the algorithm becomes more stable as $\alpha\in[\alpha_{0},2]$ get
smaller, i.e., the tails get heavier up to a certain point determined by
$\alpha_{0}$. This result implies that the stability of the algorithm, hence
the generalization error will be monotonic with respect to the tail-exponent
$\alpha$ only when $\alpha$ is large enough. This outcome sheds more light on
the experimental results presented in [3], where the relation between the
tail-exponent and the generalization error is only partially monotonic.

  * •

We prove matching lower-bounds on the stability of (5), implying that our
stability bounds are tight in the tail-exponent $\alpha$.

We support our theory on both synthetic data and real experiments conducted on
standard benchmark datasets by using fully-connected and convolutional neural
networks. All the proofs and the implementation details are provided in the
Appendix.

##  2 Notation and Background

Notation. Consider a real-valued function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ defined on $\mathbb{R}^{d}$. The
Fourier transform of $f(\theta)$ for $\theta\in\mathbb{R}^{d}$ is denoted by
$\mathcal{F}f(u)$ and is defined as,
$\mathcal{F}f(u):=\int_{\mathbb{R}^{d}}f(\theta)e^{-\mathsf{i}u^{\top}\theta}d\theta$.
Similarly, the inverse Fourier transform of a function $\hat{f}(u)$ that is
from $\mathbb{R}^{d}$ to $\mathbb{R}$ is denoted by
$\mathcal{F}^{-1}\hat{f}(\theta)$ and is defined by,
$\mathcal{F}^{-1}\hat{f}(\theta):=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\hat{f}(u)e^{\mathsf{i}u^{\top}\theta}du$.
In both of these definitions, $\mathsf{i}:=\sqrt{-1}$.

$\alpha$-stable distributions. The $\alpha$-stable distribution appears as the
limiting distribution in the generalized central limit theorems for a sum of
i.i.d. random variables with infinite variance [30]. A scalar random variable
$X$ is called symmetric $\alpha$-stable, denoted by
$X\sim\mathcal{S}\alpha\mathcal{S}(\sigma)$, if its characteristic function
takes the form:
$\mathbb{E}\left[e^{\mathsf{i}uX}\right]=\exp\left(-\sigma^{\alpha}|u|^{\alpha}\right)$,
for any $u\in\mathbb{R}$, where $\sigma>0$ is known as the scale parameter
that measures the spread of $X$ around $0$ and $\alpha\in(0,2]$ which is known
as the tail-index that determines the tail thickness of the distribution and
the tail becomes heavier as $\alpha$ gets smaller. In general, the probability
density function of a symmetric $\alpha$-stable distribution,
$\alpha\in(0,2]$, does not yield closed-form expression except for a few
special cases. When $\alpha=1$ and $\alpha=2$, $\mathcal{S}\alpha\mathcal{S}$
reduces to the Cauchy and the Gaussian distributions, respectively. When
$0<\alpha<2$, the moments are finite only up to the order $\alpha$ in the
sense that $\mathbb{E}[|X|^{p}]<\infty$ if and only if $p<\alpha$, which
implies infinite variance. Moreover, $\alpha$-stable distribution can be
extended to the high-dimensional case for random vectors. One natural
extension is the rotationally symmetric $\alpha$-stable distribution. $X$
follows a $d$-dimensional rotationally symmetric $\alpha$-stable distribution
if it admits the characteristic function $\mathbb{E}\left[e^{\mathsf{i}\langle
u,X\rangle}\right]=e^{-\sigma^{\alpha}\|u\|_{2}^{\alpha}}$ for any
$u\in\mathbb{R}^{d}$. We refer to [41] for the details of $\alpha$-stable
distributions.

Lévy processes. Lévy processes are stochastic processes with independent and
stationary increments. Their successive displacements can be viewed as the
continuous-time analogue of random walks. Lévy processes include the Poisson
process, Brownian motion, the Cauchy process, and more generally stable
processes; see e.g. [6, 41, 2]. Lévy processes in general admit jumps and have
heavy tails which are appealing in many applications; see e.g. [12]. In this
paper, we will consider the rotationally symmetric $\alpha$-stable Lévy
process $\mathrm{L}^{\alpha}_{t}$ in $\mathbb{R}^{d}$ that is defined as
follows.

  1. (i)

$\mathrm{L}^{\alpha}_{0}=0$ almost surely;

  2. (ii)

For any $t_{0}<t_{1}<\cdots<t_{N}$, the increments
$\mathrm{L}^{\alpha}_{t_{n}}-\mathrm{L}^{\alpha}_{t_{n-1}}$ are independent;

  3. (iii)

The difference $\mathrm{L}^{\alpha}_{t}-\mathrm{L}^{\alpha}_{s}$ and
$\mathrm{L}^{\alpha}_{t-s}$ have the same distribution, with the
characteristic function $\exp(-(t-s)^{\alpha}\|u\|_{2}^{\alpha})$ for $t>s$;

  4. (iv)

$\mathrm{L}^{\alpha}_{t}$ has stochastically continuous sample paths, i.e. for
any $\delta>0$ and $s\geq 0$,
$\mathbb{P}(\|\mathrm{L}^{\alpha}_{t}-\mathrm{L}^{\alpha}_{s}\|>\delta)\rightarrow
0$ as $t\rightarrow s$.

When $\alpha=2$, $\mathrm{L}^{\alpha}_{t}=\sqrt{2}\mathrm{B}_{t}$, where
$\mathrm{B}_{t}$ is the standard $d$-dimensional Brownian motion.

Ornstein-Uhlenbeck processes. Ornstein-Uhlenbeck (OU) process [48] is a
$d$-dimensional Markov and Gaussian process that satisfies the SDE:

$\mathrm{d}X_{t}=-AX_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{B}_{t},$ (6)

where $A$ and $\Sigma$ are a $d\times d$ matrices and $\mathrm{B}_{t}$ is a
standard $d$-dimensional Brownian motion. The OU process is a special case of
the Langevin equation in physics [39], and has wide applications including for
example modeling the change in organismal phenotypes in evolutionary biology
[33], and the short-rate in the interest rate modeling in finance [49]. More
generally, we can consider an OU process driven by a Lévy process, for
example, replacing $\mathrm{B}_{t}$ in (6) by a rotationally symmetric
$\alpha$-stable Lévy process $\mathrm{L}^{\alpha}_{t}$ so that

$\mathrm{d}X_{t}=-AX_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t}.$
(7)

Under some mild conditions on $A$ and $\Sigma$, the OU process $X_{t}$ in (7)
admits a unique stationary distribution that can be fully characterized; see
e.g. [43, 34].

Algorithmic stability and generalization.  In this paper, we study the
generalization of the continuous-time heavy-tailed SGD by using the tools of
algorithmic stability. Several notions of stability have been defined in the
literature of statistical learning theory [7, 14]. We will use the notion of
algorithmic stability of the randomized algorithm $\mathcal{A}$ defined in
[20]. We denote the set $\mathcal{X}_{n}$ as the set of all possible size $n$
datapoints subsampled uniformly at random from $P_{X}$.

######  Definition 1 ([20], Definition 2.1).

For a loss function $f:\mathbb{R}^{d}\times\mathcal{X}\rightarrow\mathbb{R}$,
an algorithm $\mathcal{A}$ is $\varepsilon$-uniformly stable if

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}):=\sup_{X\cong\hat{X}}\sup_{z\in\mathcal{X}}~{}\mathbb{E}\left[f(\mathcal{A}(X),z)-f(\mathcal{A}(\hat{X}),z)\right]\leq\varepsilon,$
(8)

where the first supremum is taken over data $X,\hat{X}\in\mathcal{X}_{n}$ that
differ by one element, denoted by $X\cong\hat{X}$.

Since its introduction in statistical learning theory in [7], stability based
arguments have been useful in deriving generalization bound for several
learning algorithms [5, 13, 51, 35] and have also been extended to get
generalization bound for randomized algorithm like SGD and SGLD [20, 40, 26,
37, 11, 4, 9, 29, 15]. Here below, we provide a result from [20] which relates
algorithmic stability with the generalization performance of a randomized
algorithm.

######  Theorem 2 ([20], Theorem 2.2).

Suppose that $\mathcal{A}$ is an $\varepsilon$-uniformly stable algorithm,
then the expected generalization error is bounded by

$\displaystyle\left|\mathbb{E}_{\mathcal{A},X}~{}[\hat{F}(\mathcal{A}(X),X)-F(\mathcal{A}(X))]\right|\leq\varepsilon.$
(9)

In several of recent works [16, 8, 24], high probability bounds have been
obtained using algorithmic stability bounds.

##  3 Algorithmic Stability of Heavy-Tailed SGD on Least Squares Regression

In this section, we will investigate the effects of heavy-tails on algorithmic
stability. We consider the setting of least square regression with
$f(\theta,(x,y))=(\theta^{\top}x-y)^{2}/2$. We assume that we only have the
access to the data generation distribution $P_{X}$ via the generated training
samples and our goal is to learn a parameter vector $\theta\in\mathbb{R}^{d}$
which minimize the corresponding population risk. We denote the training data
by the matrix
$X=\left[x_{1}^{\top},x_{2}^{\top},\dots,x_{i}^{\top},\dots,x_{n}^{\top}\right]\in\mathbb{R}^{n\times
d}$ and $y=[y_{1},y_{2},\dots,y_{i},\dots,y_{n}]\in\mathbb{R}^{n}$, where $n$
is the number of data points, $d$ is the dimension of the problem, and
$x_{i}\in\mathbb{R}^{d}$, $y_{i}\in\mathbb{R}$ for all $i$. Training data
points are i.i.d. from the distribution $P_{X}$. We consider the ERM problem
as defined in (1):
$\min_{\theta\in\mathbb{R}^{d}}\frac{1}{2n}\sum_{i=1}^{n}(\theta^{\top}x_{i}-y_{i})^{2}$.

In the context of algorithmic stability, we assume that we have two training
datasets $(X,y)$ and $(\hat{X},\hat{y})$ that differ in only one data point.
Without loss of generality, we have
$\hat{X}=\left[x_{1}^{\top},x_{2}^{\top},\dots,\tilde{x}_{i}^{\top},\dots,x_{n}^{\top}\right]\in\mathbb{R}^{n\times
d}$ and $\hat{y}=[y_{1},y_{2},\dots,\tilde{y}_{i},\dots,y_{n}]$. For our ERM
problem, we consider the continuous-time heavy-tailed stochastic gradient
descent, which is represented by the following two SDEs that are driven by a
rotationally symmetric $\alpha$-stable Lévy process $\mathrm{L}^{\alpha}_{t}$
in $\mathbb{R}^{d}$,

$\displaystyle\mathrm{d}\theta_{t}$
$\displaystyle=-\frac{1}{n}\left(X^{\top}X\theta_{t}-X^{\top}y\right)\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(10) $\displaystyle\mathrm{d}\hat{\theta}_{t}$
$\displaystyle=-\frac{1}{n}\left(\hat{X}^{\top}\hat{X}\hat{\theta}_{t}-\hat{X}^{\top}\hat{y}\right)\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(11)

where $\Sigma\in\mathbb{R}^{d\times d}$ is a real-valued matrix.

Under mild conditions, the SDEs (10) and (11) have unique strong solutions,
which are Markov processes and they admit unique invariant distributions [43].
Thanks to the linearity of the drifts of these SDEs, the stationary
distribution is achieved very quickly, with an exponential rate [52]. Hence,
to ease our analysis, we will assume that we have two samples from the
stationary distributions of (10) and (11), say $\theta$ and $\hat{\theta}$. In
other words, we set our learning algorithm such that it gives a random sample
from the stationary distribution of the SDE determined by the dataset, i.e.,
$\mathcal{A}_{\text{cont}}((X,y))=\theta$, and
$\mathcal{A}_{\text{cont}}((\hat{X},\hat{y}))=\hat{\theta}$, where
$\mathcal{A}_{\text{cont}}$ denotes the _continuous-time_ heavy-tailed SGD
algorithm. In the rest of this section, we will derive stability bounds for
this learning algorithm.

###  3.1 Warm-up: the need for the surrogate loss

To motivate our analysis technique, let us first consider the following simple
setting, where we set $d=1$, so that we have
$f(\theta,(x,y))=(x\theta-y)^{2}$. In this specific case, when $\alpha>1$, we
can compute the stationary distributions of (10) and (11) in an explicit form.
With a slight abuse of notation, the distribution of $\theta_{t}$ converges to
a symmetric stable law: $(\delta/s)+\mathcal{S}\alpha\mathcal{S}((\alpha
s)^{-1/\alpha})$, where $s=(1/n)\sum_{i=1}^{n}x_{i}^{2}$ and
$\delta=(1/n)\sum_{i=1}^{n}x_{i}y_{i}$ with $\delta/s$ being the mean (and the
mode) of the stationary distribution, which coincides with the ordinary least-
squares solution. Similarly, the distribution of $\hat{\theta}_{t}$ converges
to
$(\hat{\delta}/\hat{s})+\mathcal{S}\alpha\mathcal{S}((\alpha\hat{s})^{-1/\alpha})$,
where $\hat{\delta}$ and $\hat{s}$ are defined analogously.

As a first observation, assume that we have a sample from the stationary
distribution of $\theta_{t}$, such that
$\theta\sim(\delta/s)+\mathcal{S}\alpha\mathcal{S}((\alpha s)^{-1/\alpha})$.
Considering this scheme as the algorithm, i.e.,
$\mathcal{A}_{\text{cont}}((X,y))=\theta$, a simple calculation shows that

$\displaystyle\mathbb{E}_{\mathcal{A}_{\text{cont}}(X,y)}\left[f\left(\mathcal{A}_{\text{cont}}((X,y)),(X,y)\right)\right]=\mathbb{E}_{\theta,(X,y)}\left[\frac{1}{n}\sum_{i=1}^{n}(x_{i}\theta-
y_{i})^{2}\right]=+\infty,$

since the variance of $\mathcal{S}\alpha\mathcal{S}((\alpha s)^{-1/\alpha})$
is infinite whenever $\alpha<2$. Therefore, it is clear that we cannot expect
any algorithmic stability in this scheme, as long as the stability is measured
with respect to the squared loss. However, as we will show in the sequel, in
turns out that if we instead measure the stability with respect to a surrogate
loss function, which in this case would be $|x\theta-y|^{p}$ for some
$p\in[1,\alpha)$, the algorithm becomes stable, even though it is based on a
distribution that concentrates near the optimum for squared loss.

On the other hand, we notice that the means of the stationary distributions,
i.e., $\delta/s$ and $\hat{\delta}/\hat{s}$ do not interact with the tail
exponent $\alpha$. Since, our main goal is to investigate the interplay
between the tail behavior and algorithmic stability, we will ignore this term
and assume that $y_{i}=0$ almost surely for all $i$. This way, we fall back to
the SDE given in (5).

In the light of these two observations, for the general case where $d\geq 1$,
we will use the following surrogate loss function to measure stability:

$\displaystyle f(x):=f(\theta,x):=|\theta^{\top}x|^{p},\quad\text{for
some}\quad p\in[1,2],$ (12)

which generalizes the original loss function. Note that, from now on we will
drop the notation $\tilde{f}$ for denoting surrogate losses for simplicity and
use a single notation for the loss function.

###  3.2 Algorithmic stability analysis in the Fourier domain

For $d\geq 2$, unfortunately we cannot identify the stationary distributions
of (10) and (11) in an explicit form. However, by using the theory of the
characterization of the stationary distribution for an Ornstein-Uhlenbeck
process driven by a Lévy process in the literature (see [43, 34] and the
background review in the Appendix), in the next lemma, we show that we can
characterize the stationary distribution of the Ornstein-Uhlenbeck process
driven by a rotationally symmetric $\alpha$-stable Lévy process in a semi-
explicit way:

$\mathrm{d}\theta_{t}=-A\theta_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(13)

where $A$ and $\Sigma$ are $d\times d$ real matrices.

######  Lemma 3.

Assume that $A$ is a real symmetric matrix with all the eigenvalues being
positive. Then (13) admits a unique stationary distribution $\pi$ whose
characteristic function is given by

$\int_{\mathbb{R}^{d}}e^{\mathsf{i}u^{\top}x}\pi(\mathrm{d}x)=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-sA}u\right\|_{2}^{\alpha}\mathrm{d}s\right).$
(14)

While Lemma 3 provides us information about the stationary distributions of
the SDEs (10) and (11), it considers the Fourier transforms of these
distributions, which makes this setting not amenable to conventional
algorithmic stability analysis tools.

As a remedy, we perform the stability analysis directly in the Fourier domain
and use the Fourier inversion theorem to compute stability bounds for
continuous-time heavy-tailed SGD. Our main approach is based on the following
observation. Let $g:\mathbb{R}^{d}\mapsto\mathbb{R}$ be a function, and $P$,
$Q$ be random variables in $\mathbb{R}^{d}$ with respective characteristic
functions $\psi_{P}$ and $\psi_{Q}$. If the Fourier inversion theorem holds on
$g$, then $g$ is the inverse Fourier transform of $\mathcal{F}g(\cdot)$.
Hence,

$\displaystyle\mathbb{E}\left[g(P)-g(Q)\right]$
$\displaystyle=\frac{1}{(2\pi)^{d}}\mathbb{E}\left[\int_{\mathbb{R}^{d}}\left(e^{\mathsf{i}u^{\top}P}-e^{\mathsf{i}u^{\top}Q}\right)\mathcal{F}g(u)~{}\mathrm{d}u\right]=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\mathbb{E}\left[e^{\mathsf{i}u^{\top}P}-e^{\mathsf{i}u^{\top}Q}\right]\mathcal{F}g(u)~{}\mathrm{d}u$
$\displaystyle=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left(\psi_{P}(u)-\psi_{Q}(u)\right)\mathcal{F}g(u)~{}\mathrm{d}u\leq\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}|\psi_{P}(u)-\psi_{Q}(u)||\mathcal{F}g(u)|~{}\mathrm{d}u.$
(15)

Hence, (15) enables us to utilize the result given in Lemma 3 and hence gives
us a way to perform stability analysis (as given in Definition 1) in the
Fourier domain.

Algorithmic stability via characteristic function.  By setting $y=\hat{y}=0$
and invoking Lemma 3, the characteristic functions of stationary distributions
corresponding to the SDEs (10) and (11) are respectively given as follows:

$\displaystyle\psi_{\theta}(u)$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}({X}^{\top}{X})}u\right\|_{2}^{\alpha}\mathrm{d}s\right),$
(16) $\displaystyle\psi_{\hat{\theta}}(u)$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}(\hat{X}^{\top}\hat{X})}u\right\|_{2}^{\alpha}\mathrm{d}s\right).$
(17)

From the Definition 1 and from (15) and (12), we have

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\mathbb{E}\left[\left|\theta^{\top}x\right|^{p}-\left|\hat{\theta}^{\top}x\right|^{p}\right]$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\cdot\left|\mathcal{F}\left[|x^{\top}\cdot|^{p}\right](u)\right|~{}\mathrm{d}u.$
(18)

In the remainder of this section, we will consider $\Sigma=I$ for convenience
with $I$ being the identity matrix. However, we provide bounds showing the
effect of $\Sigma$ in the Appendix.

One-dimensional case ($d=1$).  We first discuss the case where $d=1$ and
report it as a separate result since its proof is simpler and more
instructive. Following (15), as a first step, we prove a lemma, which relates
the characteristic functions of the stationary distributions by upper-bounding
$|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|$. For the sake of brevity, we
present this result in the Appendix (Lemma 8). By using this intermediate
result, we next prove upper- and lower-bounds on the stability of the
continuous time heavy-tailed SGD algorithm and discuss its behavior with
respect to $\alpha$ and $p$.

######  Theorem 4.

Consider the one dimensional loss function $f(x)=|\theta x|^{p}$. For any
$x\sim P_{X}$, if we have $|x|>R$ with probability $\delta_{1}$ and for any
$X$ sampled uniformly at random from the set $\mathcal{X}_{n}$, if we have
$\|X\|_{2}^{2}\leq\sigma^{2}n$ with probability $\delta_{2}$. Then,

  1. (i)

For $\alpha\in[1,2)$, the algorithm is not stable when $p\in[\alpha,2]$ i.e.
$\varepsilon_{\text{stab}(\mathcal{A}_{\text{cont}})}$ diverges. When
$\alpha=p=2$ then
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})\leq\frac{R^{4}}{\pi\sigma^{4}n}$
with probability at least $1-\delta_{1}-2\delta_{2}$.

  2. (ii)

For $p\in[1,\alpha)$, we have the following upper bound for the algorithmic
stability,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{2R^{p+2}}{\pi\sigma^{2}n}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{1}{\alpha}\left(\frac{1}{\alpha\sigma^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)=:c(\alpha),$

which holds with probability at least $1-\delta_{1}-2\delta_{2}$. Furthermore,
for some $\alpha_{0}>1$, if we have

$\displaystyle\sigma^{2}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$
(19)

where $\phi$ is the digamma function, then the map $\alpha\mapsto c(\alpha)$
is increasing for $\alpha\in[\alpha_{0},2)$.

  3. (iii)

The stability bound is tight in $\alpha$.

Informally, this result illuminates the following facts: (i) When subject to
heavy tails, i.e., $\alpha<2$, the algorithm is stable only when a surrogate
loss is used with $p<\alpha$. (ii) For $1\leq p<\alpha<2$, the stability level
$\varepsilon_{\text{stab}}$ is upper-bounded by a function of $\alpha$, $p$,
and the variance of the data distribution $\sigma^{2}$. Furthermore (and
perhaps more surprisingly), for a given _heavy-tailedness threshold_
$\alpha_{0}\in(1,2)$, if the data variance is sufficiently large as in (19),
the stability bound becomes monotonically increasing for
$\alpha\in[\alpha_{0},2)$, which indicates that as the algorithm becomes more
stable it gets heavier-tailed. However, this relation holds as long as the
heaviness of the tails does not exceed the threshold $\alpha_{0}$. (iii) We
further show that, there exists a data distribution $P_{X}$ such that
$\varepsilon_{\text{stab}}$ is lower-bounded by a function, which also depends
on $\alpha$, $p$, and $\sigma^{2}$. In the proved lower-bound, the terms
depending on $\alpha$ have the same order as of the ones given in the upper-
bound of Theorem 4. Hence our stability bound is tight in $\alpha$. Combined
with point (ii), this result suggests that the generalization error might not
be globally monotonic with respect to the heaviness of the tails under our
modeling strategy.

Multi-dimensional case ($d\geq 2$).  Now we will focus our attention to the
case of $d$ dimensions. We follow the same route as in Theorem 4, where we
first relate the characteristic functions of the stationary distributions. We
also present this result in the Appendix (Lemma 9). Based on Lemma 9, we next
provide stability bounds for the $d$-dimensional case.

######  Theorem 5.

Consider the $d$-dimensional loss function $f(x)=|\theta^{\top}x|^{p}$ such
that $\theta,x\in\mathbb{R}^{d}$. For any $x\sim P_{X}$ if $\|x\|_{2}\leq R$,
for any $X$ sampled uniformly at random from the set $\mathcal{X}_{n}$, if
$\frac{1}{n}\|X^{\top}Xu\|_{2}\geq\sigma_{\min}\|u\|_{2}$ for
$u\in\mathbb{R}^{d}$ and for any two $X\cong\hat{X}$ sampled from
$\mathcal{X}_{n}$ generating two stochastic process given by SDEs in equations
(10) and (11),
$\|x_{i}x_{i}^{\top}-\tilde{x}_{i}\tilde{x}_{i}^{\top}\|_{2}\leq 2\sigma$
holds with high probability. Then,

  1. (i)

For $\alpha\in(1,2)$, the algorithm is not stable when $p\in[\alpha,2]$ i.e.
$\varepsilon_{\text{stab}(\mathcal{A}_{\text{cont}})}$ diverges. When
$\alpha=p=2$ then with high probability
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})\leq\frac{2R^{2}}{\pi}\frac{\sigma}{n\sigma_{\min}}$.

  2. (ii)

For $p\in[1,\alpha)$, we have the following upper bound for the algorithmic
stability,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{8R^{p}}{\pi}\frac{\sigma}{n\alpha}\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\Gamma\left(1-\frac{p}{\alpha}\right)=:c(\alpha),$

which holds with high probability. Furthermore, for some $\alpha_{0}>1$, if we
have

$\sigma_{\min}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$

where $\phi$ is the digamma function, then the map $\alpha\mapsto c(\alpha)$
is increasing for $\alpha\in[\alpha_{0},2)$.

  3. (iii)

The stability bound is tight in $\alpha$.

The conclusions of Theorem 5 are almost identical to the ones of Theorem 4,
though its proof requires a more careful analysis, especially for the lower-
bound in (iii). The main differences here are that, we need the smallest
eigenvalue of the covariance matrix of $P_{X}$, i.e., $\sigma_{\min}$ to be
large enough, and we need a different condition on the second moment $\sigma$
of the distribution. Under these conditions, we obtain very similar stability
and monotonicity properties.

As a final remark, we note that our results do not require any non-trivial
topological or statistical assumptions in comparison with [46] and [3] that
suggested a globally monotonic relation for the generalization error and the
tail exponent $\alpha$. On the other hand, the rate $1/n$ in our bounds are in
line with the existing stability literature [20, 35].

Figure 1: Results of the synthetic data experiments with varying $a$,
$\alpha$, and $d$. Each experiment was repeated $200$ times with $n=1000$. The
lines correspond to the median, and the shaded areas are the interquartile
ranges.

##  4 Experiments

Synthetic data.  We first test the implications of the theoretical findings
presented above with synthetic data experiments. We assume that $y=0$ and
$P_{X}$ is a scaled uniform distribution $\mathcal{U}(-0.5a,0.5a)$, with $a$
determining the range of the distribution. We simulate the SDE presented in
(10) by using the Euler-Maruyama discretization, which yields the following
recursion:222Due to discretization error, the stationary distribution of (20)
differs from the one of the corresponding SDE. Yet, very recently the
difference between the two stationary distributions has shown to be of order
$\eta^{2/\alpha-1}$ (tight) in the Wasserstein-1 metric for $\alpha<2$ [10],
which gives a tolerable difference for small $\eta$.

$\displaystyle\theta_{k+1}$
$\displaystyle=\theta_{k}-\eta\frac{1}{n}\left(X^{\top}X\right)\theta_{k}-\eta^{1/\alpha}E_{k+1},$
(20)

where $\eta$ is the learning-rate and each $E_{k}\in\mathbb{R}^{d}$ is a
rotationally symmetric $\alpha$-stable random vector.

In the experiments, we systematically vary $a$ as well as the tail-index of
the additive noise, $\alpha$. For all experiments we set $\eta=0.1$ and ran
the algorithm for $3000$ iterations. We set $n=1000$ and varied $d$ to be
$100$ or $250$. The order of the loss function $f$ was selected to be $p=1$.
For each experimental setting, we repeated the experiment $200$ times, where
after sampling a population $N=100000$ observations, for each replication we
sampled $n=1000$ with replacement from within this population. The
generalization error was computed to be the difference between loss computed
on the replication sample of size $n$ and the population of size $N$. To
prevent numerical issues, the noise was scaled with a constant of $0.1$ in all
experiments, which corresponds to choosing $\Sigma=(1/10)I$.

The results are presented in Figure 1, and corroborate the trend predicted by
Theorem 5. As $a$ grows, the variance of the input increases, leading the map
$\alpha\mapsto c(\alpha)$ to become increasing for $\alpha\in[\alpha_{0},2)$
for some $\alpha_{0}$. Since $c(\alpha)$ is the upper bound for stability
$\varepsilon_{\text{stab}}$, this leads to the observed ‘V-shaped’ trend in
generalization error for higher $a$ values, where the inflection point
corresponds to $\alpha_{0}$ for a given experiment setting.

Experiments on image data.  In our second set of experiments, we consider a
real image classification task, where we use plain SGD (2) _without_ adding
explicit heavy-tailed noise and monitor the effect of the heavy-tails that are
_inherently_ introduced by SGD, as shown in [19, 22]. Here, we train two fully
connected neural networks (FCN) of different depths (4 vs. 6) as well as a
convolutional neural network (CNN) on the MNIST, CIFAR-10, and CIFAR-100
datasets [28, 25]. We train these models under different, constant learning
rates ($\eta$) and with batch sizes ($b$) of $50$ or $100$, producing models
trained under a wide range of $\eta/b$ values. The models are trained until
convergence, where the convergence criteria for MNIST and CIFAR-10 is a
training negative log-likelihood (NLL) of $<5\times 10^{-5}$ and a training
accuracy of $100\%$, and for CIFAR-100 these are a NLL of $<1\times 10^{-2}$
and a training accuracy of $>99\%$.

Figure 2: Test accuracy vs. mean estimated tail-index ($\hat{\alpha}$) for
each model. Color: training $\eta/b$ ratio.

For the estimation of the trained networks’ tail indices, we used the
multivariate estimator proposed in [36, Corollary 2.4], which is previously
used in various related neural network research [46, 19, 54, 3]. Since this
estimator assumes a stable distribution, after convergence we obtained $1000$
iterations of SGD and computed the average to be used in this estimation,
based on the generalized central limit theorem [19, Corollary 11], which
demonstrated that the ergodic averages of heavy-tailed SGD iterates converge
to a multivariate stable distribution. Before estimating the parameters, we
centered the parameters with median values. Each layer’s tail-index estimation
was conducted separately, which were in turn averaged to produce a single
tail-index for every model, as in [3]. See the Appendix for further details.

Previous literature demonstrated that (i) training neural networks with larger
$\eta/b$ values lead to heavy-tailed parameters [19] and (ii) networks with
heavier-tailed parameters are more likely to generalize [46, 3]. Here, Figure
2 demonstrates that networks with highest $\alpha$ (light-tails) consistently
perform worst in terms of generalization and the performance improves as the
$\alpha$ decreases until some _threshold_. This outcome is in line with the
predictions of our theoretical results, which suggest a ‘V-shaped’ behavior
for the relation between generalization and $\alpha$, as opposed to [46, 3].
As a final remark, here the values of $\alpha$ are larger compared with the
synthetic experiments; however, we shall emphasize that such values for
$\alpha$ still indicate strong heavy tails.

##  5 Conclusion

We established novel links between the tail behavior and generalization
properties of SGD building on the notion of algorithmic stability. We focused
on quadratic optimization and considered a heavy-tailed SDE previously
proposed as a proxy to SGD dynamics. We then proved uniform stability bounds
which uncover several phenomena about the effect of the heaviness of the tails
on the generalization. We also established lower bounds which show that our
stability bounds are tight in terms of the heaviness of the tails. We
supported our theory on a variety of experiments. Finally, due to its
theoretical nature, we believe that our study does not have a direct societal
consequence.

## Acknowledgments

We would like to thank Professor Atanas Stefanov from University of Alabama,
Birmingham for helpful discussions on Fourier transforms. Anant Raj is
supported by Marie-Curie Foundation via Marie-Curie Fellowship. Mert
Gürbüzbalaban’s research is partially supported by the Office of Naval
Research Award Number N00014-21-1-2244, National Science Foundation (NSF)
CCF1814888 and NSF DMS-2053485. Lingjiong Zhu is grateful to the partial
support from a Simons Foundation Collaboration Grant and the grant NSF
DMS-2053454 from the National Science Foundation. Umut Şimşekli’s research is
supported by the French government under management of Agence Nationale de la
Recherche as part of the “Investissements d’avenir” program, reference
ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

## References

  * Ali et al. [2020] Alnur Ali, Edgar Dobriban, and Ryan Tibshirani.  The implicit regularization of stochastic gradient flow for least squares.  In _International Conference on Machine Learning_ , pages 233–244. PMLR, 2020.
  * Applebaum [2009] David Applebaum.  _Lévy Processes and Stochastic Calculus_.  Cambridge University Press, second edition, 2009.
  * Barsbey et al. [2021] Melih Barsbey, Milad Sefidgaran, Murat A Erdogdu, Gaël Richard, and Umut Şimşekli.  Heavy Tails in SGD and Compressibility of Overparametrized Neural Networks.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 29364–29378. Curran Associates, Inc., 2021.
  * Bassily et al. [2020] Raef Bassily, Vitaly Feldman, Cristóbal Guzmán, and Kunal Talwar.  Stability of stochastic gradient descent on nonsmooth convex losses.  _Advances in Neural Information Processing Systems_ , 33:4381–4391, 2020.
  * Belkin et al. [2004] Mikhail Belkin, Irina Matveeva, and Partha Niyogi.  Regularization and semi-supervised learning on large graphs.  In _International Conference on Computational Learning Theory_ , pages 624–638. Springer, 2004.
  * Bertoin [1996] Jean Bertoin.  _Lévy Processes_.  Cambridge University Press, 1996.
  * Bousquet and Elisseeff [2002] Olivier Bousquet and André Elisseeff.  Stability and generalization.  _Journal of Machine Learning Research_ , 2(Mar):499–526, 2002.
  * Bousquet et al. [2020] Olivier Bousquet, Yegor Klochkov, and Nikita Zhivotovskiy.  Sharper bounds for uniformly stable algorithms.  In _Conference on Learning Theory_ , pages 610–626. PMLR, 2020.
  * Charles and Papailiopoulos [2018] Zachary Charles and Dimitris Papailiopoulos.  Stability and generalization of learning algorithms that converge to global optima.  In _International Conference on Machine Learning_ , pages 745–754. PMLR, 2018.
  * Chen et al. [2022] Peng Chen, Changsong Deng, Rene Schilling, and Lihu Xu.  Approximation of the invariant measure of stable SDEs by an Euler–Maruyama scheme.  _arXiv preprint arXiv:2205.01342_ , 2022.
  * Chen et al. [2018] Yuansi Chen, Chi Jin, and Bin Yu.  Stability and convergence trade-off of iterative optimization algorithms.  _arXiv preprint arXiv:1804.01619_ , 2018.
  * Cont and Tankov [2004] Rama Cont and Peter Tankov.  _Financial Modelling with Jump Processes_.  Chapman and Hall/CRC, 2004.
  * Cortes et al. [2012] Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh.  Algorithms for learning kernels based on centered alignment.  _Journal of Machine Learning Research_ , 13:795–828, 2012\.
  * Elisseeff et al. [2005] Andre Elisseeff, Theodoros Evgeniou, Massimiliano Pontil, and Leslie Pack Kaelbing.  Stability of randomized learning algorithms.  _Journal of Machine Learning Research_ , 6(3):55–79, 2005.
  * Farghly and Rebeschini [2021] Tyler Farghly and Patrick Rebeschini.  Time-independent generalization bounds for SGLD in non-convex settings.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Feldman and Vondrak [2019] Vitaly Feldman and Jan Vondrak.  High probability generalization bounds for uniformly stable algorithms with nearly optimal rate.  In _Conference on Learning Theory_ , pages 1270–1279. PMLR, 2019\.
  * Flatto [2019] Leopold Flatto.  The dixie cup problem and FKG inequality.  _High Frequency_ , 2(3-4):169–174, 2019.
  * Gelfand and Shilov [1969] Izrail Moiseevic Gelfand and Georgij Evgenevic Shilov.  _Generalized Functions. Vol. 1, Properties and Operations_.  Academic Press, 1969.
  * Gürbüzbalaban et al. [2021] Mert Gürbüzbalaban, Umut Şimşekli, and Lingjiong Zhu.  The heavy-tail phenomenon in SGD.  In _International Conference on Machine Learning_ , 2021.
  * Hardt et al. [2016] Moritz Hardt, Ben Recht, and Yoram Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In _International Conference on Machine Learning_ , pages 1225–1234. PMLR, 2016.
  * He and Tao [2020] Fengxiang He and Dacheng Tao.  Recent advances in deep learning theory.  _arXiv preprint arXiv:2012.10931_ , 2020.
  * Hodgkinson and Mahoney [2021] Liam Hodgkinson and Michael Mahoney.  Multiplicative noise and heavy tails in stochastic optimization.  In _International Conference on Machine Learning_ , pages 4262–4274. PMLR, 2021.
  * Hodgkinson et al. [2021] Liam Hodgkinson, Umut Şimşekli, Rajiv Khanna, and Michael W Mahoney.  Generalization properties of stochastic optimizers via trajectory analysis.  _arXiv preprint arXiv:2108.00781_ , 2021.
  * Klochkov and Zhivotovskiy [2021] Yegor Klochkov and Nikita Zhivotovskiy.  Stability and deviation optimal risk bounds with convergence rate $o(1/n)$.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Krizhevsky [2009] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, University of Toronto, 2009.
  * Kuzborskij and Lampert [2018] Ilja Kuzborskij and Christoph Lampert.  Data-dependent stability of stochastic gradient descent.  In _International Conference on Machine Learning_ , pages 2815–2824. PMLR, 2018.
  * Latorre et al. [2021] Fabian Latorre, Leello Tadesse Dadi, Paul Rolland, and Volkan Cevher.  The effect of the intrinsic dimension on the generalization of quadratic classifiers.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * LeCun et al. [2010] Yann LeCun, Corinna Cortes, and CJ Burges.  MNIST handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * Lei and Ying [2020] Yunwen Lei and Yiming Ying.  Fine-grained analysis of stability and generalization for stochastic gradient descent.  In _International Conference on Machine Learning_ , pages 5809–5819. PMLR, 2020.
  * Lévy [1937] Paul Lévy.  Théorie de l’addition des variables aléatoires.  _Gauthiers-Villars, Paris_ , 1937.
  * Liu and Theodorou [2019] Guan-Horng Liu and Evangelos A Theodorou.  Deep learning theory review: An optimal control and dynamical systems perspective.  _arXiv preprint arXiv:1908.10920_ , 2019.
  * Martin and Mahoney [2019] Charles H Martin and Michael W Mahoney.  Traditional and heavy tailed self regularization in neural network models.  In _International Conference on Machine Learning_ , 2019.
  * Martins [1994] Emilia P. Martins.  Estimating the rate of phenotypic evolution from comparative data.  _The American Naturalist_ , 144(2):193–209, 1994\.
  * Masuda [2004] Hiroki Masuda.  On multidimensional Ornstein-Uhlenbeck processes driven by a general Lévy process.  _Bernoulli_ , 10(1):97–120, 2004.
  * Maurer and Jaakkola [2005] Andreas Maurer and Tommi Jaakkola.  Algorithmic stability and meta-learning.  _Journal of Machine Learning Research_ , 6(6), 2005.
  * Mohammadi et al. [2015] Mohammad Mohammadi, Adel Mohammadpour, and Hiroaki Ogata.  On estimating the tail index and the spectral measure of multivariate $\alpha$-stable distributions.  _Metrika_ , 78(5):549–561, 2015.
  * Mou et al. [2018] Wenlong Mou, Liwei Wang, Xiyu Zhai, and Kai Zheng.  Generalization bounds of SGLD for non-convex learning: Two theoretical viewpoints.  In _Conference on Learning Theory_ , pages 605–638. PMLR, 2018.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Pavliotis [2014] Grigorios A Pavliotis.  _Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations_ , volume 60.  Springer, 2014.
  * Raginsky et al. [2017] Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky.  Non-convex learning via stochastic gradient Langevin dynamics: A nonasymptotic analysis.  In _Conference on Learning Theory_ , pages 1674–1703. PMLR, 2017\.
  * Samorodnitsky and Taqqu [1994] Gennady Samorodnitsky and Murad S. Taqqu.  _Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance_.  Chapman & Hall, New York, 1994.
  * Sathe and Upadhye [2020] Aastha M. Sathe and N. S. Upadhye.  Estimation of the parameters of multivariate stable distributions.  _Communications in Statistics - Simulation and Computation_ , 0(0):1–18, July 2020.  ISSN 0361-0918.  doi: 10.1080/03610918.2020.1784432.  URL https://doi.org/10.1080/03610918.2020.1784432.  Publisher: Taylor & Francis _eprint: https://doi.org/10.1080/03610918.2020.1784432.
  * Sato and Yamazato [1984] Ken-iti Sato and Makoto Yamazato.  Operator-self-decomposable distributions as limit distributions of processes of Ornstein-Uhlenbeck type.  _Stochastic Processes and their Applications_ , 17:73–100, 1984.
  * Simonyan and Zisserman [2015] Karen Simonyan and Andrew Zisserman.  Very Deep Convolutional Networks for Large-Scale Image Recognition.  _arXiv:1409.1556 [cs]_ , April 2015.  URL http://arxiv.org/abs/1409.1556.
  * Şimşekli et al. [2019] Umut Şimşekli, Mert Gürbüzbalaban, Thanh Huy Nguyen, Gaël Richard, and Levent Sagun.  On the heavy-tailed theory of stochastic gradient descent for deep neural networks.  _arXiv preprint arXiv:1912.00018_ , 2019.
  * Şimşekli et al. [2020] Umut Şimşekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu.  Hausdorff dimension, heavy tails, and generalization in neural networks.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 5138–5151. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/37693cfc748049e45d87b8c7d8b9aacd-Paper.pdf.
  * Tzagkarakis et al. [2018] George Tzagkarakis, John P Nolan, and Panagiotis Tsakalides.  Compressive sensing of temporally correlated sources using isotropic multivariate stable laws.  In _2018 26th European Signal Processing Conference (EUSIPCO)_ , pages 1710–1714. IEEE, 2018.
  * Uhlenbeck and Ornstein [1930] George Eugene Uhlenbeck and Leonard S. Ornstein.  On the theory of Brownian motion.  _Physical Review_ , 36(5):823–841, 1930.
  * Vasicek [1977] Oldrich Vasicek.  An equilibrium characterization of the term structure.  _Journal of Financial Economics_ , 5(2):177–188, 1977.
  * Wang et al. [2021] Hongjian Wang, Mert Gürbüzbalaban, Lingjiong Zhu, Umut Şimşekli, and Murat A Erdogdu.  Convergence rates of stochastic gradient descent under infinite noise variance.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Wu and Cheng [2021] Xinxing Wu and Qiang Cheng.  Algorithmic stability and generalization of an unsupervised feature selection algorithm.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Xie and Zhang [2020] Longjie Xie and Xicheng Zhang.  Ergodicity of stochastic differential equations with jumps and singular coefficients.  _Annales de l’Institut Henri Poincaré, Probabilités et Statistiques_ , 56(1):175–229, 2020.
  * Zhang et al. [2020] Jingzhao Zhang, Sai Praneeth Karimireddy, Andreas Veit, Seungyeon Kim, Sashank Reddi, Sanjiv Kumar, and Suvrit Sra.  Why are adaptive methods good for attention models?  In _Advances in Neural Information Processing Systems_ , volume 33, pages 15383–15393, 2020.
  * Zhou et al. [2020] Pan Zhou, Jiashi Feng, Chao Ma, Caiming Xiong, Steven Chu Hong Hoi, and Weinan E.  Towards theoretically understanding why SGD generalizes better than ADAM in deep learning.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 21285–21296. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/f3f27a324736617f20abbf2ffd806f6d-Paper.pdf.

Algorithmic Stability of Heavy-Tailed Stochastic Gradient Descent on Least
Squares
Appendix

The Appendix is organized as follows:

  * •

In Section A, we provide the background details about characterizing the
stationary distribution of a Lévy-driven OU process.

  * •

In Section B, we provide the proofs for the $1$-dimensional case.

  * •

In Section C, we prove the results for least square in $d$-dimension.

  * •

In Section D, we extend the $d$-dimensional result for general preconditioner
PSD $\Sigma$.

  * •

In Section E, we discuss useful results which we utilize in proving our
results for $d$-dimensional case.

  * •

In Section F, we provide further details about our experimental setup.

##  Appendix A Characterizing the Stationary Distribution of a Lévy-Driven OU
Process

In this section, we review the technical background of characterizing the
stationary distribution of an Ornstein-Uhlenbeck process driven by a general
Lévy process. Consider an Ornstein-Uhlenbeck process driven by a general Lévy
process

$\mathrm{d}\theta_{t}=-A\theta_{t}\mathrm{d}t+\mathrm{d}\mathrm{Z}_{t},$ (21)

where $\mathrm{Z}_{t}$ is a general Lévy process. One particular example is
$\mathrm{Z}_{t}=\Sigma\mathrm{L}^{\alpha}_{t}$ so that

$\mathrm{d}\theta_{t}=-A\theta_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t}.$
(22)

Under some regularity conditions on $A$ and the Lévy measure of $Z$,
$\theta_{t}$ in (21) admits a unique invariant distribution $\pi$, and the
class of all possible $\pi$’s forms the class of all $A$-self-decomposable
distributions; see [34] and the references therein.

Let $d\in\mathbb{N}$ and $\mathrm{Z}_{t}$ is a $d$-dimensional Lévy process
such that $Z_{0}=0$ a.s. and $\mathrm{Z}_{t}$ admits the generating triplet
$(b,C,\nu)$, that is, $b\in\mathbb{R}^{d}$, $C$ is a $d\times d$ symmetric
non-negative definite matrix and $\nu$ is a $\sigma$-finite measure on
$\mathbb{R}^{d}$ satisfying $\nu(\\{0\\})=0$ and
$\int_{\mathbb{R}^{d}}\min(1,\|z\|^{2})\nu(dz)<\infty$ for which
$\mathrm{Z}_{t}$ has the characteristic function

$\varphi_{t}(u):=\mathbb{E}\left[e^{\mathsf{i}\langle
u,\mathrm{Z}_{t}\rangle}\right]=\exp\left(t\left\\{iu^{\top}b-\frac{1}{2}u^{\top}Cu+\int_{\mathbb{R}^{d}}\left(e^{iu^{\top}z}-1-iu^{\top}z1_{\|z\|\leq
1}\right)\nu(dz)\right\\}\right)$ (23)

for any $u\in\mathbb{R}^{d}$ and $t>0$.

######  Lemma 6 (Theorems 4.1. and 4.2. in [43]).

Assume that $A$ is a $d\times d$ matrix such that the real parts of all its
eigenvalues are positive. Moreover, assume that
$\int_{\|z\|>1}\log\|z\|\nu(dz)<\infty$. Then, $\theta_{t}$ in (21) admits a
unique invariant distribution $\pi$ whose characteristic function is given by

$\int_{\mathbb{R}^{d}}e^{\mathsf{i}\langle
u,x\rangle}\pi(dx)=\exp\left(\int_{0}^{\infty}\log\varphi_{1}\left(e^{-sA^{\top}}u\right)ds\right),$
(24)

for any $u\in\mathbb{R}^{d}$. In particular, the generating triplet of the
limiting distribution is $(b_{\infty},C_{\infty},\nu_{\infty})$, where

$\displaystyle
b_{\infty}=A^{-1}b+\int_{\mathbb{R}^{d}}\int_{0}^{\infty}e^{-sA}z\left(1_{\|e^{-sA}z\|\leq
1}-1_{\|z\|\leq 1}\right)ds\nu(dz),$ (25) $\displaystyle
C_{\infty}=\int_{0}^{\infty}e^{-sA}Ce^{-sA^{\top}}ds,$ (26)
$\displaystyle\nu_{\infty}(E)=\int_{0}^{\infty}\nu\left(e^{sA}E\right)ds,\qquad\text{for
any $E\in\mathcal{B}(\mathbb{R}^{d})$}.$ (27)

By using Lemma 6, we can easily obtain the following result.

######  Lemma (Restatement of Lemma 3).

Assume that $A$ is a real symmetric matrix with all the eigenvalues being
positive. Then (22) admits a unique stationary distribution

$\int_{\mathbb{R}^{d}}e^{\mathsf{i}\langle
u,x\rangle}\pi(dx)=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-sA}u\right\|_{2}^{\alpha}ds\right),\qquad\text{for
any $u\in\mathbb{R}^{d}$}.$ (28)

###### Proof.

In our (22), $A$ is a real symmetric matrix with positive eigenvalues and
$\mathrm{Z}_{t}=\Sigma\mathrm{L}^{\alpha}_{t}$, where
$\mathrm{L}^{\alpha}_{t}$ is a rotationally symmetric $\alpha$-stable Lévy
process and moreover, $\alpha$-stable Lévy measure satisfies the condition
$\int_{\|z\|>1}\log\|z\|\nu(dz)<\infty$, so that the condition in Lemma 6 is
satisfied, and we conclude that (22) admits a unique stationary distribution
say $\pi$.

Moreover, in our case, the characteristic function of
$\mathrm{Z}_{t}=\Sigma\mathrm{L}^{\alpha}_{t}$ is given by

$\varphi_{1}(u)=\mathbb{E}\left[e^{\mathsf{i}\left\langle
u,\Sigma\mathrm{L}^{\alpha}_{1}\right\rangle}\right]=\mathbb{E}\left[e^{\mathsf{i}\left\langle\Sigma^{\top}u,\mathrm{L}^{\alpha}_{1}\right\rangle}\right]e^{-\|\Sigma^{\top}u\|_{2}^{\alpha}}.$
(29)

Therefore, by Lemma 6, the unique invariant distribution $\pi$ of (22) has the
following characteristic function:

$\int_{\mathbb{R}^{d}}e^{\mathsf{i}\langle
u,x\rangle}\pi(dx)=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-sA}u\right\|_{2}^{\alpha}ds\right).$
(30)

This completes the proof. ∎

######  Remark 7.

(i) We can also characterize the generating triplet for the limiting
distribution $\pi$ in the above lemma according to Lemma 6. In our case, $b=0$
in (25), $C=0$ in (26), and $\nu$ is the Lévy measure for
$\Sigma\mathrm{L}^{\alpha}_{t}$ in (27).

(ii) In general, it is not possible to further simplify the expression (28)
except for some special cases. For example, when $A=I$. For any $s\geq 0$,
$e^{-sI}=\sum_{k=0}^{\infty}\frac{(-sI)^{k}}{k!}=\sum_{k=0}^{\infty}\frac{(-s)^{k}}{k!}I=e^{-s}I$.
Therefore, when $A=I$, we can compute from (28) that

$\int_{\mathbb{R}^{d}}e^{\mathsf{i}\langle
u,x\rangle}\pi(dx)=\exp\left(-\int_{0}^{\infty}e^{-s\alpha}\left\|\Sigma^{\top}u\right\|_{2}^{\alpha}ds\right)=e^{-\frac{1}{\alpha}\left\|\Sigma^{\top}u\right\|_{2}^{\alpha}}.$
(31)

Hence, in this special case, the limiting distribution is
$\alpha^{\frac{-1}{\alpha}}\Sigma\mathrm{L}^{\alpha}_{1}$.

##  Appendix B Proofs for the 1-dimensional Case

In the next lemma, we first bound the difference between the characteristic
functions of the stationary distributions.

######  Lemma 8.

For two matrices $X\in\mathbb{R}^{n}$ and $\hat{X}\in\mathbb{R}^{n}$ as
defined earlier, the absolute value of difference between the characteristic
function for stationary distribution at any $u\in\mathbb{R}$ corresponding to
one-dimensional rotation invariant processes in equations (16) and (17)
($d=1$) is bounded as

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle\leq\frac{x_{i}^{2}-\tilde{x}_{i}^{2}}{\|\hat{X}\|^{2}}\left(|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\right)\right).$

###### Proof.

We can compute that

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle=\left|\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds\right)-\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|\hat{X}\|_{2}^{2}}u\right|^{\alpha}ds\right)\right|$
$\displaystyle=\left|\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds\right)\left(1-\exp\left(\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|{X}\|_{2}^{2}}u\right|^{\alpha}ds\right.\right.\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.\left.\left.-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|\hat{X}\|_{2}^{2}}u\right|^{\alpha}ds\right)\right)\right|$
(32)
$\displaystyle\leq\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds\right)\left|\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|\hat{X}\|_{2}^{2}}u\right|^{\alpha}ds\right|$
$\displaystyle=\exp\left(-|u|^{\alpha}\int_{0}^{\infty}e^{-s\frac{\alpha}{n}\|X\|_{2}^{2}}ds\right)|u|^{\alpha}\left|\int_{0}^{\infty}e^{-s\frac{\alpha}{n}\|X\|_{2}^{2}}ds-\int_{0}^{\infty}e^{-s\frac{\alpha}{n}\|\hat{X}\|_{2}^{2}}ds\right|$
$\displaystyle=|u|^{\alpha}\exp\left(-|u|^{\alpha}\int_{0}^{\infty}e^{-s\frac{\alpha}{n}\|X\|_{2}^{2}}ds\right)\left|\frac{n}{\alpha\|X\|_{2}^{2}}-\frac{n}{\alpha\|\hat{X}\|_{2}^{2}}\right|$
$\displaystyle=|u|^{\alpha}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|_{2}^{2}}\right)\left|\frac{n}{\alpha\|X\|_{2}^{2}}-\frac{n}{\alpha\|\hat{X}\|_{2}^{2}}\right|$
$\displaystyle=|u|^{\alpha}\frac{n|x_{i}^{2}-\tilde{x}_{i}^{2}|}{\alpha\|X\|_{2}^{2}\|\hat{X}\|_{2}^{2}}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|_{2}^{2}}\right)$
$\displaystyle=\frac{|x_{i}^{2}-\tilde{x}_{i}^{2}|}{\|\hat{X}\|^{2}}\left(|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\right)\right),$

which completes the proof. ∎

######  Theorem (Restatement of Theorem 4).

Consider the one dimensional loss function $f(x)=|\theta x|^{p}$. For any
$x\sim P_{X}$, if we have $|x|>R$ with probability $\delta_{1}$ and for any
$X$ sampled uniformly at random from the set $\mathcal{X}_{n}$, if we have
$\|X\|_{2}^{2}\leq\sigma^{2}n$ with probability $\delta_{2}$. Then,

  1. (i)

For $\alpha\in[1,2)$, the algorithm is not stable when $p\in[\alpha,2]$ i.e.
$\varepsilon_{\text{stab}(\mathcal{A}_{\text{cont}})}$ diverges. When
$\alpha=p=2$ then
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})\leq\frac{R^{4}}{\pi\sigma^{4}n}$
with probability at least $1-\delta_{1}-2\delta_{2}$.

  2. (ii)

For $p\in[1,\alpha)$, we have the following upper bound for the algorithmic
stability,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{2R^{p+2}}{\pi\sigma^{2}n}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{1}{\alpha}\left(\frac{1}{\alpha\sigma^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)=:c(\alpha),$

which holds with probability at least $1-\delta_{1}-2\delta_{2}$. Furthermore,
for some $\alpha_{0}>1$, if we have

$\displaystyle\sigma^{2}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$
(33)

where $\phi$ is the digamma function, then the map $\alpha\mapsto c(\alpha)$
is increasing for $\alpha\in[\alpha_{0},2)$.

  3. (iii)

The stability bound is tight in $\alpha$.

###### Proof.

Closed-form expression for the Fourier transform of function
$f(\theta)=|\theta x|^{p}$ has been given in [18]. However, we provide here
the result for the sake of completeness. Let us compute the Fourier transform
of the function $\theta\to|\theta x|^{p}$ for $p\in[1,2)$.

$\displaystyle\int_{-\infty}^{\infty}|\theta
x|^{p}e^{-\mathsf{i}u\theta}d\theta=$
$\displaystyle|x|^{p}\int_{-\infty}^{\infty}|\theta|^{p}e^{-\mathsf{i}u\theta}~{}d\theta$
$\displaystyle=$
$\displaystyle|x|^{p}\int_{0}^{\infty}(e^{\mathsf{i}u\theta}+e^{-\mathsf{i}u\theta})\theta^{p}~{}d\theta$
$\displaystyle=$
$\displaystyle|x|^{p}\left[\int_{0}^{\infty}e^{\mathsf{i}u\theta}\theta^{p}~{}d\theta+\int_{0}^{\infty}e^{-\mathsf{i}u\theta}\theta^{p}~{}d\theta\right]$
$\displaystyle=$ $\displaystyle
2|x|^{p}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|u|^{p+1}}.$
(34)

From [18],

$\displaystyle\int_{-\infty}^{\infty}|\theta
x|^{2}e^{-\mathsf{i}u\theta}d\theta=\frac{-2x^{2}}{u^{2}}\delta(u),$

where $\delta(u)$ is the Dirac-delta function.

First, we get the result for $p\in[1,2)$. We utilize the result from Lemma 8
and equation (34) to get,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\left[\frac{1}{(2\pi)}\int_{\mathbb{R}}|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\left|\int_{-\infty}^{\infty}|\theta
x|^{p}e^{-\mathsf{i}u\theta}~{}d\theta\right|~{}du\right]$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{|x_{i}^{2}-\tilde{x}_{i}^{2}|}{\|\hat{X}\|^{2}}\frac{n}{\alpha\|X\|^{2}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{\mathbb{R}}\left(|u|^{\alpha}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\right)\right)\frac{1}{|u|^{p+1}}~{}du$
$\displaystyle\leq\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2|x|^{p+2}}{\pi\|\hat{X}\|^{2}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{0}^{\infty}u^{\alpha-p-1}\frac{n}{\alpha\|X\|^{2}}\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|^{2}}\right)~{}du.$

In the above integral, by substituting $u^{\alpha}\frac{n}{\alpha\|X\|^{2}}$
with $t$ so that

$\displaystyle dt=u^{\alpha-1}\frac{n}{\|X\|^{2}}du,~{}\text{and
}\frac{1}{u^{p}}=\left(\frac{n}{\alpha\|X\|^{2}}\right)^{\frac{p}{\alpha}}t^{-p/\alpha},$
(35)

we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2|x|^{p+2}}{\pi\alpha\|\hat{X}\|^{2}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{n}{\alpha\|X\|^{2}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt.$
(36)

It is clear that, the above integral diverges for $p\geq\alpha$, hence the
algorithm is not stable for $p\in[1,2)$. Now, we check the case for $p=2$. For
$p=2$, we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{|x|^{4}}{\pi\|\hat{X}\|^{2}}\frac{n}{\|X\|^{2}}\int_{0}^{\infty}u^{\alpha-2}\exp\left(-u^{\alpha}\frac{n}{\|X\|_{2}^{2}}\right)\delta(u)~{}du.$

The above integral clearly diverges for $\alpha<2$. However, when $\alpha=2$,
then

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{|x|^{4}}{\pi\|\hat{X}\|^{2}}\frac{n}{\|X\|^{2}}.$

If we have $|x|>R$ for any $x\sim P_{X}$ with probability $\delta_{1}$ and
$\|X\|_{2}^{2}\leq\sigma^{2}n$ for any $X$ sampled uniformly from the set
$\mathcal{X}_{n}$ with probability $\delta_{2}$, then for $\alpha=2$ and
$p=2$,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{R^{4}}{\pi\sigma^{4}n},$ (37)

with probability at least $1-\delta_{1}-2\delta_{2}$. This proves the part (i)
of our result.

Next, let us prove the part (ii). From equation (36), for $p<\alpha$, we have

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2|x|^{p+2}}{\pi\alpha\|\hat{X}\|^{2}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{n}{\alpha\|X\|^{2}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2|x|^{p+2}}{\pi\alpha\|\hat{X}\|^{2}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{n}{\alpha\|X\|^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

If we have $|x|>R$ for any $x\sim P_{X}$ with probability $\delta_{1}$ and
$\|X\|_{2}^{2}\leq\sigma^{2}n$ for any $X$ sampled uniformly at random from
the set $\mathcal{X}_{n}$ with probability $\delta_{2}$, then with probability
at least $1-\delta_{1}-2\delta_{2}$, the following holds:

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{2R^{p+2}}{\pi\sigma^{2}n}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{1}{\alpha}\left(\frac{1}{\alpha\sigma^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)=c(\alpha).$

Now, consider the function,

$\displaystyle\Lambda(\alpha)=\frac{1}{\alpha}\left(\frac{1}{\alpha\sigma^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

We can compute that

$\displaystyle\partial_{\alpha}\log\Lambda(\alpha)=\frac{p}{\alpha^{2}}\left[\log\alpha+\log\sigma^{2}-1-\frac{\alpha}{p}+\phi\left(1-\frac{p}{\alpha}\right)\right],$

where $\phi$ is the digamma function. For any arbitrary $\alpha_{0}$, if we
choose

$\sigma^{2}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$

then $\partial_{\alpha}\log\Lambda(\alpha)>0$ for $\alpha\in[\alpha_{0},2)$.
Hence, for all $\alpha_{1},\alpha_{2}\in[\alpha_{0},2)$,
$\alpha_{1}<\alpha_{2}\Rightarrow\Lambda(\alpha_{1})\leq\Lambda(\alpha_{2})$.
This proves that $c(\alpha)$ is an increasing map in $\alpha$.

(iii). Now we show that the bound on the stability is tight for some
appropriately chosen $P_{X}$. Note that we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\left[\frac{1}{(2\pi)}\int_{\mathbb{R}}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|\left|\int_{-\infty}^{\infty}|\theta
x|^{p}e^{-\mathsf{i}u\theta}~{}d\theta\right|~{}du\right]$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\underbrace{\left[\frac{1}{\pi}\int_{0}^{\infty}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|\left|\int_{-\infty}^{\infty}|\theta
x|^{p}e^{-\mathsf{i}u\theta}~{}d\theta\right|~{}du\right]}_{:=\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})}.$

Hence, let us consider $u\geq 0$. From equation (32), we have

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle=\Bigg{|}\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds\right)$
$\displaystyle\qquad\qquad\cdot\left(1-\exp\left(\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|{X}\|_{2}^{2}}u\right|^{\alpha}ds-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|\hat{X}\|_{2}^{2}}u\right|^{\alpha}ds\right)\right)\Bigg{|}$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left|e^{-s\frac{1}{n}\|X\|_{2}^{2}}u\right|^{\alpha}ds\right)\left(1-\exp\left(-|u|^{\alpha}\left[\frac{n}{\alpha\|\hat{X}\|_{2}^{2}}-\frac{n}{\alpha\|{X}\|_{2}^{2}}\right]\right)\right)$
$\displaystyle=\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|_{2}}\right)\left(1-\exp\left(-|u|^{\alpha}\left[\frac{n\overbrace{(x_{i}^{2}-\tilde{x}_{i}^{2})}^{:=\delta}}{\alpha\|X\|_{2}^{2}\|\hat{X}\|_{2}^{2}}\right]\right)\right)$
$\displaystyle=\exp\left(-|u|^{\alpha}\frac{n}{\alpha\|X\|_{2}}\right)\left[\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k!}|u|^{k\alpha}\left(\frac{n\delta}{\alpha\|X\|_{2}^{2}\|\hat{X}\|_{2}^{2}}\right)^{k}\right].$

Hence,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\left[\frac{1}{\pi}\int_{0}^{\infty}|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\left|\int_{-\infty}^{\infty}|\theta
x|^{p}e^{-\mathsf{i}u\theta}~{}d\theta\right|~{}du\right]$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\exp\left(-u^{\alpha}\frac{n}{\alpha\|X\|_{2}^{2}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdot\left[\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k!}u^{k\alpha-p-1}\left(\frac{n\delta}{\alpha\|X\|_{2}^{2}\|\hat{X}\|_{2}^{2}}\right)^{k}\right]du.$
(38)

To simplify the above term, we do need to compute the integral:

$\displaystyle\int_{0}^{\infty}\exp\left(-\frac{u^{\alpha}n}{\alpha\|X\|_{2}^{2}}\right)\frac{u^{k\alpha}}{u^{p+1}}~{}du.$

Let us use the substitution:

$\displaystyle\frac{u^{\alpha}n}{\alpha\|X\|_{2}^{2}}=t$

such that

$\displaystyle dt=\frac{u^{\alpha-1}n}{\|X\|_{2}^{2}}du\text{~{}~{}~{}~{}and
~{}}u=\left(\frac{\alpha\|X\|_{2}^{2}}{n}\right)^{1/\alpha}t^{1/\alpha}.$

Hence,

$\displaystyle\int_{0}^{\infty}\exp\left(-\frac{u^{\alpha}n}{\alpha\|X\|_{2}^{2}}\right)\frac{u^{k\alpha}}{u^{p+1}}~{}du$
$\displaystyle=\int_{0}^{\infty}e^{-t}u^{\alpha(k-1)-p}\frac{\|X\|_{2}^{2}}{n}~{}dt$
$\displaystyle=\frac{\|X\|_{2}^{2}}{n}\left(\frac{\alpha\|X\|_{2}^{2}}{n}\right)^{k-1-\frac{p}{\alpha}}\int_{0}^{\infty}e^{-t}t^{k-1-\frac{p}{\alpha}}dt$
$\displaystyle=\frac{\|X\|_{2}^{2}}{n}\left(\frac{\alpha\|X\|_{2}^{2}}{n}\right)^{k-1-\frac{p}{\alpha}}\Gamma\left(k-\frac{p}{\alpha}\right).$

This implies,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\exp\left(-\frac{|u|^{\alpha}n}{\alpha\|X\|_{2}^{2}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\sum_{k=1}^{\infty}\left[\frac{(-1)^{k+1}}{k!}u^{k\alpha-p-1}\left(\frac{n\delta}{\alpha\|X\|_{2}^{2}\|\hat{X}\|_{2}^{2}}\right)^{k}\right]~{}du$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\sum_{k=1}^{\infty}\left[\frac{(-1)^{k+1}}{k!}\frac{\delta^{k}}{\alpha\|\hat{X}\|_{2}^{2k}}\left(\frac{n}{\alpha\|X\|_{2}^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(k-\frac{p}{\alpha}\right)\right]$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\sum_{k=1}^{\infty}(-1)^{k+1}\gamma_{k},$

where
$\gamma_{k}=\frac{1}{k!}\frac{\delta^{k}}{\alpha\|\hat{X}\|_{2}^{2k}}\left(\frac{n}{\alpha\|X\|_{2}^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(k-\frac{p}{\alpha}\right)$.
Let us now compute $\frac{\gamma_{k+1}}{\gamma_{k}}$. We have

$\displaystyle\frac{\gamma_{k+1}}{\gamma_{k}}=\frac{\delta}{\|\hat{X}\|_{2}^{2}}\frac{\Gamma\left(k+1-\frac{p}{\alpha}\right)}{(k+1)\Gamma\left(k-\frac{p}{\alpha}\right)}=\frac{\delta}{\|\hat{X}\|_{2}^{2}}\frac{k-\frac{p}{\alpha}}{k+1}=\frac{\delta}{\|\hat{X}\|_{2}^{2}}\left(1-\frac{1+\frac{p}{\alpha}}{k+1}\right).$

Hence, we have the following,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\left[\gamma_{1}\sum_{k=1}^{\infty}\left(-\frac{\delta}{\|\hat{X}\|_{2}^{2}}\right)^{k-1}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\gamma_{1}\left(1+\frac{p}{\alpha}\right)\sum_{k=1}^{\infty}\left(-\frac{\delta}{\|\hat{X}\|_{2}^{2}}\right)^{k-1}\prod_{j=1}^{k-1}\left(\frac{1}{j+1}\right)\right]$
$\displaystyle\geq\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\gamma_{1}\sum_{k=1}^{\infty}\left(-\frac{\delta}{\|\hat{X}\|_{2}^{2}}\right)^{k-1}$
$\displaystyle=\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{\gamma_{1}}{1+\frac{\delta}{\|\hat{X}\|_{2}^{2}}}$
$\displaystyle\geq\frac{2|x|^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{\delta}{\alpha\|\hat{X}\|_{2}^{2}}\left(\frac{n}{\alpha\|X\|_{2}^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)$
(39)

Now, let us assume that $P_{X}$ is a distribution with discrete support in
range $\sigma^{2}$ to $R$ with $C$ number of support points equally spaced.
Hence, with probability $(1-1/C)$, $\delta\geq c$ for some positive constant
$c$. Hence, with with high probability,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})\geq\frac{2\sigma^{2p}c}{R^{2}\pi\alpha
n}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{\alpha}\left(\frac{1}{\alpha
R^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

This completes the proof. ∎

##  Appendix C Proofs for Least-Square in $d$-dimension

We start by proving the following lemma, relating the characteristic functions
of the two distributions.

######  Lemma 9.

For two matrices $X\in\mathbb{R}^{n\times d}$ and
$\hat{X}\in\mathbb{R}^{n\times d}$ as defined earlier, the absolute value of
difference between the characteristic functions of the stationary
distributions at any $u\in\mathbb{R}^{d}$ corresponding to $d$ dimensional
rotation invariant processes in equations (16) and (17) is bounded as

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle\leq\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right),$

where $\sigma_{\min}$ is the smaller of the smallest of singular values of the
matrices $\frac{1}{n}X^{\top}X$ and $\frac{1}{n}\hat{X}^{\top}\hat{X}$, and
$x_{i}x_{i}^{\top}-\tilde{x}_{i}{\tilde{x}_{i}}^{\top}=\sigma_{1}v_{1}v_{1}^{\top}+\sigma_{2}v_{2}v_{2}^{\top}$
where $v_{1}$ and $v_{2}$ are orthogonal vectors.

###### Proof.

We can compute that

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle=\left|\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)-\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right)\right|$
$\displaystyle\leq\underbrace{\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)}_{:=B}\underbrace{\left|\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right|}_{:=C}.$

We first consider the term $C$ in the above equation. From Lemma 12, we have
for two positive numbers $a$ and $b$, and for some $1\leq\alpha\leq 2$, we
have

$\displaystyle|a^{\alpha}-b^{\alpha}|\leq|a-b|(a^{\alpha-1}+b^{\alpha-1}).$

Now,

$\displaystyle C$
$\displaystyle=\left|\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right|$
$\displaystyle=\left|\int_{0}^{\infty}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}-\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}\right)ds\right|$
$\displaystyle\leq\int_{0}^{\infty}\left|\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}-\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}\right|\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle\leq\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u-e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle=\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}\left(I-e^{s\frac{1}{n}X^{\top}X-s\frac{1}{n}\hat{X}^{\top}\hat{X}}\right)u\right\|_{2}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds.$

Now, we have from the definitions,

$\displaystyle
X^{\top}X-\hat{X}^{\top}\hat{X}=x_{i}x_{i}^{\top}-\tilde{x}_{i}\tilde{x}_{i}^{\top}.$

Hence,

$\displaystyle C$
$\displaystyle\leq\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}\left(I-e^{s\frac{1}{n}\left(x_{i}x_{i}^{\top}-\tilde{x}_{i}\tilde{x}_{i}^{\top}\right)}\right)u\right\|_{2}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds.$

We recall that the 2-norm $\|\cdot\|_{2}$ for a matrix
$D\in\mathbb{R}^{d\times d}$ is defined as follows:

$\displaystyle\|D\|_{2}=\sup_{u,\|u\|_{2}=1}\|Du\|_{2}$

for $u\in\mathbb{R}^{d}$. We notice that

$\min\left(\frac{1}{n}\left\|\left(X^{\top}X\right)u\right\|_{2},\frac{1}{n}\left\|\left(\hat{X}^{\top}\hat{X}\right)u\right\|_{2}\right)\geq\sigma_{\min}\|u\|_{2}.$

Hence,

$\displaystyle C$
$\displaystyle\leq\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}\left\|I-e^{-s\frac{1}{n}\left(\tilde{x}_{i}\tilde{x}_{i}^{\top}-{x}_{i}x_{i}^{\top}\right)}\right\|_{2}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle=\int_{0}^{\infty}\underbrace{\left\|I-e^{s\frac{1}{n}\left(x_{i}x_{i}^{\top}-\tilde{x_{i}}\tilde{x_{i}}^{\top}\right)}\right\|_{2}}_{:=D}\underbrace{\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}+\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)}_{:=E}ds.$
(40)

Let us consider the term $D$ first.
$\tilde{x}_{i}{\tilde{a}_{i}}^{\top}-x_{i}x_{i}^{\top}$ is a rank 2 matrix.
Consider the two non-zero eigenvalues of this matrix are $\sigma_{1}$ and
$\sigma_{2}$. Hence,
$\tilde{x}_{i}{\tilde{a}_{i}}^{\top}-x_{i}x_{i}^{\top}=\sigma_{1}v_{1}v_{1}^{\top}+\sigma_{2}v_{2}v_{2}^{\top}$
where $v_{1}$ and $v_{2}$ are the eigenvectors. Then,

$\displaystyle
I-e^{-s\frac{1}{n}\left(\tilde{x}_{i}{\tilde{a}_{i}}^{\top}-x_{i}x_{i}^{\top}\right)}=\left(1-e^{\frac{-s\sigma_{1}}{n}}\right)v_{1}v_{1}^{\top}+\left(1-e^{\frac{-s\sigma_{2}}{n}}\right)v_{2}v_{2}^{\top}.$

Hence,

$\displaystyle
D=\left\|I-e^{-s\frac{1}{n}\left(\tilde{x}_{i}{\tilde{a}_{i}}^{\top}-x_{i}x_{i}^{\top}\right)}\right\|_{2}$
$\displaystyle\leq\left\|\left(1-e^{\frac{-s\sigma_{1}}{n}}\right)v_{1}v_{1}^{\top}\right\|_{2}+\left\|\left(1-e^{\frac{-s\sigma_{2}}{n}}\right)v_{2}v_{2}^{\top}\right\|_{2}$
$\displaystyle\leq\frac{s\sigma_{1}}{n}+\frac{s\sigma_{2}}{n},$

where $v_{1}$ and $v_{2}$ are orthogonal vectors with
$\|v_{1}\|_{2}=\|v_{2}\|_{2}=1$ and $v_{1}^{\top}v_{2}=0$. By definition, we
have

$\displaystyle\frac{1}{n}\left\|X^{\top}Xu\right\|_{2}\geq\sigma_{\min}\|u\|_{2},\quad\text{and}\quad\frac{1}{n}\left\|\hat{X}^{\top}\hat{X}u\right\|_{2}\geq\sigma_{\min}\|u\|_{2}.$

This gives,

$\displaystyle\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}\leq
e^{-s\sigma_{\min}},\quad\text{and}\quad\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}\leq
e^{-s\sigma_{\min}},$

which implies that

$\displaystyle E\leq 2\|u\|_{2}^{\alpha}e^{-s\alpha\sigma_{\min}}.$ (41)

Therefore,

$\displaystyle C$
$\displaystyle\leq\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n}\int_{0}^{\infty}se^{-s\alpha\sigma_{\min}}ds=\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n\alpha\sigma_{\min}}.$
(42)

Hence, we have

$\displaystyle|\psi_{1}(u)-\psi_{2}(u)|$
$\displaystyle\leq\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)$
$\displaystyle\leq\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\|u\|_{2}^{\alpha}\int_{0}^{\infty}e^{-s\alpha\sigma_{\min}}ds\right),$
(43)

where the last inequality is due to the definition of $\sigma_{\min}$. Hence
we conclude that

$\displaystyle|\psi_{1}(u)-\psi_{2}(u)|$
$\displaystyle\leq\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right),$
(44)

which completes the proof. ∎

######  Theorem (Restatement of Theorem 5).

Consider the $d$-dimensional loss function $f(x)=|\theta^{\top}x|^{p}$ such
that $\theta,x\in\mathbb{R}^{d}$. For any $x\sim P_{X}$ if $\|x\|_{2}\leq R$,
for any $X$ sampled uniformly at random from the set $\mathcal{X}_{n}$, if
$\frac{1}{n}\|X^{\top}Xu\|_{2}\geq\sigma_{\min}\|u\|_{2}$ for
$u\in\mathbb{R}^{d}$ and for any two $X\cong\hat{X}$ sampled from
$\mathcal{X}_{n}$ generating two stochastic process given by SDEs in equations
(10) and (11),
$\|x_{i}x_{i}^{\top}-\tilde{x}_{i}\tilde{x}_{i}^{\top}\|_{2}\leq 2\sigma$
holds with high probability. Then,

  1. (i)

For $\alpha\in(1,2)$, the algorithm is not stable when $p\in[\alpha,2]$ i.e.
$\varepsilon_{\text{stab}(\mathcal{A}_{\text{cont}})}$ diverges. When
$\alpha=p=2$ then with high probability
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})\leq\frac{2R^{2}}{\pi}\frac{\sigma}{n\sigma_{\min}}$.

  2. (ii)

For $p\in[1,\alpha)$, we have the following upper bound for the algorithmic
stability,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{8R^{p}}{\pi}\frac{\sigma}{n\alpha}\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\Gamma\left(1-\frac{p}{\alpha}\right)=c(\alpha),$

which holds with high probability. Furthermore, for some $\alpha_{0}>1$, if we
have

$\sigma_{\min}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$

where $\phi$ is the digamma function, then the map $\alpha\rightarrow
c(\alpha)$ is increasing for $\alpha\in[\alpha_{0},2)$.

  3. iii

The stability bound is tight in $\alpha$.

###### Proof.

We have $d$-dimensional loss function for an $x\in\mathbb{R}^{d}$ sampled
uniformly at random from $P_{X}$, $f(\theta)=|\theta^{\top}x|^{p}$ Let us
denote the Fourier transform of $f$, $\mathcal{F}f(u)$ as $h(u)$. For an
orthogonal matrix $A$ such that $Ae_{1}=\frac{x}{\|x\|_{2}}$, we have from the
results in Lemma 13,

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|u_{1}|^{p+1}}\quad\text{for~{}}p\in[1,2),$
(45)

and

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{1},u_{2},\cdots,u_{d})\frac{2}{u_{1}^{2}}\quad\text{for~{}}p=2,$
(46)

where $\delta$ is the Dirac-delta function. Let us first consider the case
when $p\in[1,2)$. From equation (18),

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right||h(u)|~{}du$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du.$

In the above equation, let us apply the change of variable $u=Av$ and use
result from Lemma 13 (equations (45)) and we get the following,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|Av\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|Av\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}(Av)^{\top}\theta}~{}d\theta\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|h(Av)\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left[\left(\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\right)\right.~{}$
$\displaystyle\qquad\qquad\qquad\left.\cdot\left(\left|2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(v_{2},\cdots,v_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|v_{1}|^{p+1}}\right|\right)\right]~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{(\sigma_{1}+\sigma_{2})}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{-\infty}^{\infty}|v_{1}|^{\alpha}\exp\left(\frac{-|v_{1}|^{\alpha}}{\alpha\sigma_{\min}}\right)\frac{1}{|v_{1}|^{p+1}}~{}dv_{1}$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha\sigma_{\min}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{0}^{\infty}v_{1}^{\alpha-p-1}\exp\left(\frac{-v_{1}^{\alpha}}{\alpha\sigma_{\min}}\right)~{}dv_{1}.$

In the above integral, by substituting
$\frac{v^{\alpha}}{\alpha\sigma_{\min}}$ with $t$ so that

$\displaystyle dt=v^{\alpha-1}\frac{1}{\sigma_{\min}}dv,~{}\text{and
}\frac{1}{v^{p}}=\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}t^{-p/\alpha},$
(47)

we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}$
$\displaystyle\qquad\qquad\qquad\cdot\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt.$
(48)

It is clear that, the above integral diverge for $p\geq\alpha$, hence the
algorithm is not stable for $p\in[1,2)$. Now, we check the case for $p=2$. For
$p=2$, we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{2}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du.$

In the above equation, we make change of variable $u=Av$ and use result from
Lemma 13 (equations (46)), we get the following,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|Av\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\exp\left(-\frac{\|Av\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{2}e^{\mathsf{i}(Av)^{\top}\theta}~{}d\theta\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|h(Av)\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2}{\pi}\int_{\mathbb{R}^{d}}\frac{(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\|x\|_{2}^{2}\delta(v_{1},v_{2},\cdots,v_{d})\frac{2}{v_{1}^{2}}~{}dv.$

The above integral clearly diverges for $\alpha<2$. However, when $\alpha=2$,
then

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{\|x\|_{2}^{2}}{\pi}\frac{(\sigma_{1}+\sigma_{2})}{n\sigma_{\min}}.$

Now, if $\sigma$ is the upper bound on $\sigma_{1}$ and $\sigma_{2}$ for all
$X\cong\hat{X}\in\mathcal{X}_{n}$ and $\|x\|_{2}\leq R$ for $x\sim P_{X}$ with
high probability, then,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{2R^{2}}{\pi}\frac{\sigma}{n\sigma_{\min}},$

holds with high probability. This proves the part (i) of our claim.

Next, we will prove part (ii) when $p<\alpha$. We have from equation (48),

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

Now, if $\sigma$ is the upper bound on $\sigma_{1}$ and $\sigma_{2}$ for all
$X\cong\hat{X}\in\mathcal{X}_{n}$ and $\|x\|_{2}\leq R$ for $x\sim P_{X}$ with
high probability then,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{8R^{p}}{\pi}\frac{\sigma}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)$

holds with high probability. Now, consider the function,

$\displaystyle\Lambda(\alpha)=\frac{1}{\alpha}\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

We can compute that

$\displaystyle\partial_{\alpha}\log\Lambda(\alpha)=\frac{p}{\alpha^{2}}\left[\log\alpha+\log\sigma_{\min}-1-\frac{\alpha}{p}+\phi\left(1-\frac{p}{\alpha}\right)\right],$

where $\phi$ is the digamma function. For any arbitrary $\alpha_{0}$, if we
choose

$\sigma_{\min}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)\right),$

then $\partial_{\alpha}\log\Lambda(\alpha)>0$ for $\alpha\in[\alpha_{0},2)$.
Hence, for all $\alpha_{1},\alpha_{2}\in[\alpha_{0},2)$,
$\alpha_{1}<\alpha_{2}\Rightarrow\Lambda(\alpha_{1})\leq\Lambda(\alpha_{2})$.
This proves that $c(\alpha)$ is an increasing map in $\alpha$.

This completes the proof till part (ii). Now, we will prove tightness result
in $\alpha$. Let us have the following construction. Consider a one
dimensional distribution $P_{X}$ supported in a ring such that the density
function $\int_{A}p(x)dx\leq\eta$ such that
$A=\\{x:|x|\geq\sqrt{\sigma_{\min}d\log d}\text{~{}or~{}}|x|\leq R\\}$. The
empirical covariance matrix $X^{\top}X$ is a diagonal matrix. Hence, from the
results in [17]. with high probability $1-\delta$, we have

$\displaystyle\frac{1}{n}\|X^{\top}Xu\|\geq\sigma_{\min}\|u\|_{2}^{2}.$

Exact expression for $\delta$ is given in [17]. Similarly, for the dataset
$\hat{X}$, the similar condition holds,

$\displaystyle\frac{1}{n}\|\hat{X}^{\top}\hat{X}u\|\geq\sigma_{\min}\|u\|_{2}^{2}$

with high probability $1-\delta$. We have,

$\displaystyle\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|$
$\displaystyle=\left|\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)-\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right)\right|$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)\left|1-\exp\left(-\int_{0}^{\infty}\left|\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}-\left\|e^{-s\frac{1}{n}{X}^{\top}{X}}u\right\|_{2}^{\alpha}\right|ds\right)\right|.$
(49)

From equation (18),

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right||h(u)|~{}du$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)\right|\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du.$

In the above equation, let us apply the change of variable $u=Av$ where $A$ is
the orthogonal matrix defined earlier and we get the following,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(Av)-\psi_{\hat{\theta}}(Av)\right|\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}(Av)^{\top}\theta}~{}d\theta\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(Av)-\psi_{\hat{\theta}}(Av)\right|\left|h(Av)\right|~{}dv$

Since, $A$ is orthogonal matrix, we can see that

$\displaystyle\left|\psi_{\theta}(Av)-\psi_{\hat{\theta}}(Av)\right|=\left|\psi_{\theta}(v)-\psi_{\hat{\theta}}(v)\right|.$

Hence,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(v)-\psi_{\hat{\theta}}(v)\right|\left|h(Av)\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(v)-\psi_{\hat{\theta}}(v)\right|\Bigg{(}\Bigg{|}2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(v_{2},\cdots,v_{d})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|v_{1}|^{p+1}}\Bigg{|}\Bigg{)}dv.$

Let us denote,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left|\psi_{\theta}(v)-\psi_{\hat{\theta}}(v)\right|\Bigg{(}\Bigg{|}2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(v_{2},\cdots,v_{d})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|v_{1}|^{p+1}}\Bigg{|}\Bigg{)}dv.$

Now, we use the property of Dirac-delta function. From our construction,
$X^{\top}X$ and $\hat{X}^{\top}\hat{X}$ are diagonal matrices. Let us denote
$X^{\top}X=\text{diag}(a_{1},a_{2},\cdots,a_{d})$. Similarly, we denote
$\hat{X}^{\top}\hat{X}=\text{diag}(\hat{a}_{1},\hat{a}_{2},\cdots,\hat{a}_{d})$.
Hence, we have

$\displaystyle\exp\left(-\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}v\right\|_{2}^{\alpha}ds\right)$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left(\sum_{i=1}^{d}e^{-2(s/n)a_{i}}v_{i}^{2}\right)^{\frac{\alpha}{2}}ds\right).$

From the construction, the matrix $X^{\top}X$ and $\hat{X}^{\top}X$ are both
diagonal and differ at two diagonal elements with probability $(1-1/d)$. They
differ at one diagonal element with probability $1/d$. Let’s assume that
$x_{i}$ has non-zero element at dimension 1 and $\tilde{x}_{i}$ has non-zero
element either at dimension 1 or at 2 (without loss of generality). Hence,
with high probability,

$\displaystyle\left|1-\exp\left(-\int_{0}^{\infty}\left|\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}-\left\|e^{-s\frac{1}{n}{X}^{\top}{X}}u\right\|_{2}^{\alpha}\right|ds\right)\right|$
$\displaystyle=1-\exp\left(-\int_{0}^{\infty}\left|\left(\sum_{i=1}^{d}e^{-2(s/n)a_{i}}v_{i}^{2}\right)^{\alpha/2}-\left(\sum_{i=1}^{d}e^{-2(s/n)\hat{a}_{i}}v_{i}^{2}\right)^{\alpha/2}\right|ds\right)$

Combining everything together and using the property of Dirac-delta function
we get,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\Bigg{[}\exp\left(-\int_{0}^{\infty}e^{-(s\alpha/n)a_{1}}|v_{1}|^{\alpha}ds\right)$
$\displaystyle\qquad\qquad\qquad\cdot\left[1-\exp\left(-\int_{0}^{\infty}\left(e^{-(s\alpha/n)a_{1}}-e^{-(s\alpha/n)\tilde{a}_{1}}\right)|v_{1}|^{\alpha}ds\right)\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdot\Bigg{(}\Bigg{|}2\|x\|_{2}^{p}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|v_{1}|^{p+1}}\Bigg{|}\Bigg{)}\Bigg{]}dv_{1}$
$\displaystyle=\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\Bigg{[}\exp\left(-\int_{0}^{\infty}e^{-(s\alpha/n)a_{1}}v_{1}^{\alpha}ds\right)$
$\displaystyle\qquad\qquad\qquad\cdot\frac{1}{v_{1}^{p+1}}\left[1-\exp\left(-\int_{0}^{\infty}\left(e^{-(s\alpha/n)a_{1}}-e^{-(s\alpha/n)\tilde{a}_{1}}\right)v_{1}^{\alpha}ds\right)\right]\Bigg{]}dv_{1}$
$\displaystyle=\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\Bigg{[}\exp\left(-v_{1}^{\alpha}\frac{n}{\alpha
a_{1}}\right)$
$\displaystyle\qquad\qquad\cdot\frac{1}{v_{1}^{p+1}}\left[1-\exp\left(-v_{1}^{\alpha}\Big{[}\frac{n}{\alpha
a_{1}}-\frac{n}{\alpha\hat{a}_{1}}\Big{]}\right)\right]\Bigg{]}dv_{1}$
$\displaystyle=\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\Bigg{[}\frac{1}{v_{1}^{p+1}}\exp\left(-v_{1}^{\alpha}\frac{n}{\alpha
a_{1}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\cdot\Big{[}1-\exp\left(\frac{n(\|x_{i}\|^{2}-\|\tilde{x}_{i}\|^{2})}{\alpha
a_{1}\hat{a}_{1}}\right)\Big{]}\Bigg{]}dv_{1}.$

Let us denote $\delta=\|x_{i}\|^{2}-\|\tilde{x}_{i}\|^{2}$. Hence,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdot\int_{0}^{\infty}\Bigg{[}\frac{1}{v_{1}^{p+1}}\exp\left(-v_{1}^{\alpha}\frac{n}{\alpha
a_{1}}\right)\left[1-\exp\left(\frac{n\delta}{\alpha
a_{1}\hat{a}_{1}}\right)\right]\Bigg{]}dv_{1}$
$\displaystyle=\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\int_{0}^{\infty}\exp\left(-v_{1}^{\alpha}\frac{n}{\alpha
a_{1}}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdot\left[\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k!}v_{1}^{k\alpha-p-1}\left(\frac{n\delta}{\alpha
a_{1}\hat{a}_{1}}\right)^{k}\right]dv_{1}.$

The above equation is just reduction to the computation of one dimensional
case which we did in equation (38). We apply similar argument that we did
apply in computing the lower bound in equation (38). Hence, with high
probability, we get (equation (39))

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})\geq\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{\delta}{\alpha\hat{a}_{1}^{2}}\left(\frac{n}{\alpha
a_{1}^{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)$

Here, we also assume that $P_{X}$ is a distribution with discrete support in
range $\sigma_{\min}$ to $R$ with $C$ number of support points equally spaced.
Hence, with probability $(1-1/C)$, $\delta\geq c$ for some positive constant
$c$.

By construction and the result from [17], we know that $nCd\log d\geq
a_{1}\geq n\sigma_{\min}$ for some positive constant $C$ with high
probability. This also holds for $\hat{a}_{1}$. Hence, for some positive
constant $C_{1}$ and $C_{2}$ ($C_{1}$ and $C_{2}$ has dependence on the
dimension) ,

$\displaystyle\Phi_{x,X,\hat{X}}(\mathcal{A}_{\text{cont}})\geq\frac{C_{1}}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\left(\frac{1}{\alpha
C_{2}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

This completes the proof. ∎

##  Appendix D Case for General P.S.D $\Sigma$ (Preconditioning)

In this section, we would discuss the effect of general positive semidefinite
matrix $\Sigma$. As in equations (50) and (51), we consider two SDEs
corresponding to a rotationally symmetric $\alpha$-stable Lévy process
$\mathrm{L}^{\alpha}_{t}$ in $\mathbb{R}^{d}$,

$\displaystyle\mathrm{d}\theta_{t}$
$\displaystyle=-\frac{1}{n}(X^{\top}X)\theta_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(50) $\displaystyle\mathrm{d}\hat{\theta}_{t}$
$\displaystyle=-\frac{1}{n}(\hat{X}^{\top}\hat{X})\hat{\theta}_{t}\mathrm{d}t+\Sigma\mathrm{d}\mathrm{L}^{\alpha}_{t},$
(51)

where $\Sigma\in\mathbb{R}^{d\times d}$ is a real valued P.S.D matrix. The
corresponding characteristic functions are given by as in equations (52) and
(53) (see Lemma 3),

$\displaystyle\psi_{\theta}(u)$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}({X}^{\top}{X})}u\right\|_{2}^{\alpha}ds\right),$
(52) $\displaystyle\psi_{\hat{\theta}}(u)$
$\displaystyle=\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}(\hat{X}^{\top}\hat{X})}u\right\|_{2}^{\alpha}ds\right).$
(53)

We assume that the largest and smallest eigenvalues of the matrix $\Sigma$ is
$\lambda_{\max}$ and $\lambda_{\min}$.

######  Lemma 10.

For two matrices $X\in\mathbb{R}^{n\times d}$ and
$\hat{X}\in\mathbb{R}^{n\times d}$ as defined earlier, the absolute value of
difference between the characteristic functions of the stationary
distributions at any $u\in\mathbb{R}^{d}$ corresponding to $d$ dimensional
rotation invariant processes in equations (52) and (53) is bounded as

$\displaystyle|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\leq\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right),$

where $\sigma_{\min}$ is the smaller of the smallest of singular values of the
matrices $\frac{1}{n}X^{\top}X$ and $\frac{1}{n}\hat{X}^{\top}\hat{X}$, and
$x_{i}x_{i}^{\top}-\tilde{x}_{i}{\tilde{x}_{i}}^{\top}=\sigma_{1}v_{1}v_{1}^{\top}+\sigma_{2}v_{2}v_{2}^{\top}$
where $v_{1}$ and $v_{2}$ are orthogonal vectors.

###### Proof.

We can compute that

$\displaystyle|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|$
$\displaystyle=\left|\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)-\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right)\right|$
$\displaystyle\leq\underbrace{\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)}_{:=B}\underbrace{\left|\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right|}_{:=C}.$

We first consider the term $C$ in the above equation. From Lemma 12, we have
for two positive numbers $a$ and $b$, and for some $1\leq\alpha\leq 2$, we
have

$\displaystyle|a^{\alpha}-b^{\alpha}|\leq|a-b|(a^{\alpha-1}+b^{\alpha-1}).$

Now,

$\displaystyle C$
$\displaystyle=\left|\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}ds\right|$
$\displaystyle=\left|\int_{0}^{\infty}\left(\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}-\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha}\right)ds\right|$
$\displaystyle\leq\int_{0}^{\infty}\left|\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}-\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}\right|\left(\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle\leq\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u-\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}\left(\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle=\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}\left(I-e^{s\frac{1}{n}X^{\top}X-s\frac{1}{n}\hat{X}^{\top}\hat{X}}\right)u\right\|_{2}\left(\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|\Sigma^{\top}e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds$
$\displaystyle\leq\lambda_{\max}^{\alpha}\underbrace{\int_{0}^{\infty}\left\|e^{-s\frac{1}{n}X^{\top}X}\left(I-e^{s\frac{1}{n}X^{\top}X-s\frac{1}{n}\hat{X}^{\top}\hat{X}}\right)u\right\|_{2}\left(\left\|e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha-1}+\left\|e^{-s\frac{1}{n}\hat{X}^{\top}\hat{X}}u\right\|_{2}^{\alpha-1}\right)ds}_{\text{This
term has been analyzed as an upper bound on term C in
Lemma~{}\ref{lem:dd_char_func}(Equation~{}\eqref{eq:upper_bnd_C}).}}.$

Using the result directly from equation (42), we have,

$\displaystyle
C\leq\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n\alpha\sigma_{\min}}.$
(54)

Let us consider the term B. Using the similar arguments as in Lemma 9
(equation (43)), we have,

$\displaystyle\exp\left(-\int_{0}^{\infty}\left\|\Sigma^{\top}e^{-s\frac{1}{n}X^{\top}X}u\right\|_{2}^{\alpha}ds\right)$
$\displaystyle\leq\exp\left(-\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}\int_{0}^{\infty}e^{-s\alpha\sigma_{\min}}ds\right)$
$\displaystyle=\exp\left(-\frac{\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right).$
(55)

Hence, we have the final result,

$\displaystyle|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\leq\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right),$
(56)

which completes the proof. ∎

######  Theorem 11.

Consider the $d$-dimensional loss function $f(x)=|\theta^{\top}x|^{p}$ such
that $\theta,x\in\mathbb{R}^{d}$. For any $x\sim P_{X}$ if $\|x\|_{2}\leq R$,
for any $X$ sampled uniformly at random from the set $\mathcal{X}_{n}$, if
$\frac{1}{n}\|X^{\top}Xu\|_{2}\geq\sigma_{\min}\|u\|_{2}$ for
$u\in\mathbb{R}^{d}$ and for any two $X\cong\hat{X}$ sampled from
$\mathcal{X}_{n}$ generating two stochastic process given by SDEs in equations
(50) and (51),
$\|x_{i}x_{i}^{\top}-\tilde{x}_{i}\tilde{x}_{i}^{\top}\|_{2}\leq 2\sigma$
holds with high probability. Then,

  1. (i)

For $\alpha\in(1,2)$, the algorithm is not stable when $p\in[\alpha,2]$ i.e.
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$ diverges. When
$\alpha=p=2$ then with high probability
$\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})\leq\frac{2R^{2}}{\pi}\frac{\lambda_{\max}^{2}\sigma}{n\sigma_{\min}}$.

  2. (ii)

For $p\in[1,\alpha)$, we have the following upper bound for the algorithmic
stability,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{8R^{p}}{\pi}\lambda_{\min}^{p}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\frac{\sigma}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)$
$\displaystyle=c(\alpha),$

which holds with high probability. Furthermore, for some $\alpha_{0}>1$, if we
have

$\sigma_{\min}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)-\alpha_{0}^{2}\log\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)\right),$

where $\phi$ is the digamma function, then the map $\alpha\rightarrow
c(\alpha)$ is increasing for $\alpha\in[\alpha_{0},2)$.

###### Proof.

We have $d$-dimensional loss function for an $x\in\mathbb{R}^{d}$ sampled
uniformly at random from $P_{X}$, $f(\theta)=|\theta^{\top}x|^{p}$ Let us
denote the Fourier transform of $f$, $\mathcal{F}f(u)$ as $h(u)$. For an
orthogonal matrix $A$ such that $Ae_{1}=\frac{x}{\|x\|_{2}}$, we have from the
results in Lemma 13,

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|u_{1}|^{p+1}}\quad\text{for~{}}p\in[1,2),$
(57)

and

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{1},u_{2},\cdots,u_{d})\frac{2}{u_{1}^{2}}\quad\text{for~{}}p=2,$
(58)

where $\delta$ is the Dirac-delta function. Let us first consider the case
when $p\in[1,2)$. From equation (18) and Lemma 10,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)||h(u)|~{}du$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}|\psi_{\theta}(u)-\psi_{\hat{\theta}}(u)|\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{\alpha}^{\alpha}}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\exp\left(-\frac{\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du.$

In the above equation, we make change of variable $u=Av$ and use the result
from Lemma 13 (equation (57)) to get the following,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|Av\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\exp\left(-\frac{\lambda_{\min}^{\alpha}\|Av\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{p}e^{\mathsf{i}(Av)^{\top}\theta}~{}d\theta\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|h(Av)\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\left[\left(\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\right)\right.~{}$
$\displaystyle\qquad\qquad\left.\cdot\left(\left|2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(v_{2},\cdots,v_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|v_{1}|^{p+1}}\right|\right)\right]~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2\|x\|_{2}^{p}}{\pi}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\frac{\lambda_{\max}^{\alpha}(\sigma_{1}+\sigma_{2})}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{-\infty}^{\infty}|v_{1}|^{\alpha}\exp\left(\frac{-\lambda_{\min}^{\alpha}|v_{1}|^{\alpha}}{\alpha\sigma_{\min}}\right)\frac{1}{|v_{1}|^{p+1}}~{}dv_{1}$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\frac{\lambda_{\max}^{\alpha}(\sigma_{1}+\sigma_{2})}{n\alpha\sigma_{\min}}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{0}^{\infty}v_{1}^{\alpha-p-1}\exp\left(\frac{-\lambda_{\min}^{\alpha}v_{1}^{\alpha}}{\alpha\sigma_{\min}}\right)~{}dv_{1}.$

In the above integral, by substituting
$\frac{\lambda_{\min}^{\alpha}v^{\alpha}}{\alpha\sigma_{\min}}$ with $t$ so
that

$\displaystyle
dt=\lambda_{\min}^{\alpha}v^{\alpha-1}\frac{1}{\sigma_{\min}}dv,~{}\text{and
}\frac{1}{v^{p}}=\lambda_{\min}^{p}\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}t^{-p/\alpha},$
(59)

we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\lambda_{\min}^{p}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}$
$\displaystyle\qquad\qquad\qquad\cdot\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt.$
(60)

It is clear that, the above integral diverge for $p\geq\alpha$, hence the
algorithm is not stable for $p\in[1,2)$. Now, we check the case for $p=2$. For
$p=2$, we have,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|u\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\cdot\exp\left(-\frac{\lambda_{\min}^{\alpha}\|u\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{2}e^{\mathsf{i}u^{\top}\theta}~{}d\theta\right|~{}du.$

In the above equation, we make change of variable $u=Av$ aand use the result
from Lemma 13 (equation (58)) to get the following,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\Bigg{\\{}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|Av\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\exp\left(-\frac{\lambda_{\min}^{\alpha}\|Av\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|\int_{\mathbb{R}^{d}}|\theta^{\top}x|^{2}e^{\mathsf{i}(Av)^{\top}\theta}~{}d\theta\right|~{}dv\Bigg{\\}}$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{1}{(2\pi)^{d}}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{2(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\left|h(Av)\right|~{}dv$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{2}{\pi}\int_{\mathbb{R}^{d}}\lambda_{\max}^{\alpha}\frac{(\sigma_{1}+\sigma_{2})\|v\|_{2}^{\alpha}}{n\alpha\sigma_{\min}}\exp\left(-\frac{\lambda_{\min}^{\alpha}\|v\|_{2}^{\alpha}}{\alpha\sigma_{\min}}\right)\|x\|_{2}^{2}\delta(v_{1},v_{2},\cdots,v_{d})\frac{2}{v_{1}^{2}}~{}dv.$

In the last equation, we used the result from Lemma 13. The above integral
clearly diverges for $\alpha<2$. However, when $\alpha=2$, then

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{\|x\|_{2}^{2}}{\pi}\frac{\lambda_{\max}^{2}(\sigma_{1}+\sigma_{2})}{n\sigma_{\min}}.$

Now, if $\sigma$ is the upper bound on $\sigma_{1}$ and $\sigma_{2}$ for all
$X\cong\hat{X}\in\mathcal{X}_{n}$ and $\|x\|_{2}\leq R$ for $x\sim P_{X}$ with
high probability then,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle\leq\frac{2R^{2}}{\pi}\frac{\lambda_{\max}^{2}\sigma}{n\sigma_{\min}},$

holds with high probability. This proves part (i) of our claim.

Next, we will prove part (ii) when $p<\alpha$. We have from equation (60),

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\Bigg{\\{}\frac{4\|x\|_{2}^{p}}{\pi}\lambda_{\min}^{p}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}\Gamma(p+1)$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\cdot\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\int_{0}^{\infty}t^{-p/\alpha}e^{-t}~{}dt\Bigg{\\}}$
$\displaystyle=\sup_{X\cong\hat{X}}\sup_{x\in\mathcal{X}}\frac{4\|x\|_{2}^{p}}{\pi}\lambda_{\min}^{p}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\frac{(\sigma_{1}+\sigma_{2})}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

Now, if $\sigma$ is the upper bound on $\sigma_{1}$ and $\sigma_{2}$ for all
$X\cong\hat{X}\in\mathcal{X}_{n}$ and $\|x\|_{2}\leq R$ for $x\sim P_{X}$ with
high probability then,

$\displaystyle\varepsilon_{\text{stab}}(\mathcal{A}_{\text{cont}})$
$\displaystyle=\frac{8R^{p}}{\pi}\lambda_{\min}^{p}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\frac{\sigma}{n\alpha}\Gamma(p+1)\cos\left(\frac{(p-1)\pi}{2}\right)\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right)$

holds with high probability. Now, consider the function,

$\displaystyle\Lambda(\alpha)=\frac{1}{\alpha}\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)^{\alpha}\left(\frac{1}{\alpha\sigma_{\min}}\right)^{\frac{p}{\alpha}}\Gamma\left(1-\frac{p}{\alpha}\right).$

We can compute that

$\displaystyle\partial_{\alpha}\log\Lambda(\alpha)=\log\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)+\frac{p}{\alpha^{2}}\left[\log\alpha+\log\sigma_{\min}-1-\frac{\alpha}{p}+\phi\left(1-\frac{p}{\alpha}\right)\right],$

where $\phi$ is the digamma function. For any arbitrary $\alpha_{0}$, if we
choose

$\sigma_{\min}\geq\exp\left(1+\frac{2}{p}-\log\alpha_{0}-\phi\left(1-\frac{p}{\alpha_{0}}\right)-\alpha_{0}^{2}\log\left(\frac{\lambda_{\max}}{\lambda_{\min}}\right)\right),$

then $\partial_{\alpha}\log\Lambda(\alpha)>0$ for $\alpha\in[\alpha_{0},2)$.
Hence, for all $\alpha_{1},\alpha_{2}\in[\alpha_{0},2)$,
$\alpha_{1}<\alpha_{2}$ it follows that
$\Lambda(\alpha_{1})\leq\Lambda(\alpha_{2})$. This proves that $c(\alpha)$ is
an increasing map in $\alpha$. This completes the proof. ∎

##  Appendix E Useful Results

Here below, we provide a trivial result which we have utilized in the proof of
main result.

######  Lemma 12.

For any two positive numbers $a$ and $b$, and for some $0<\alpha\leq 2$, we
have

$\displaystyle|a^{\alpha}-b^{\alpha}|\leq|a-b|(a^{\alpha-1}+b^{\alpha-1}).$
(61)

###### Proof.

When $a=b>0$, the result is obviously true. Without loss of generality, let us
assume that $a>b>0$ and by considering the RHS of (61), we get

$\displaystyle|a-b|(a^{\alpha-1}+b^{\alpha-1})$
$\displaystyle=(a-b)(a^{\alpha-1}+b^{\alpha-1})$
$\displaystyle=a^{\alpha}+ab^{\alpha-1}-a^{\alpha-1}b-b^{\alpha}$
$\displaystyle=|a^{\alpha}-b^{\alpha}|+ab^{\alpha-1}-a^{\alpha-1}b.$

Since, we have assumed that $a>b>0$ and $\alpha>0$, hence
$ab^{\alpha-1}-a^{\alpha-1}b>0$ always which essentially means,

$\displaystyle|a^{\alpha}-b^{\alpha}|\leq|a-b|(a^{\alpha-1}+b^{\alpha-1}).$

Same argument can be given while assuming $b>a>0$. This completes the proof. ∎

######  Lemma 13 (Fourier transform of $|\theta^{\top}x|^{p}$).

Consider the function $f(\theta)=|\theta^{\top}x|^{p}$ for $p\in[1,2]$ and
$h(u)$ denotes the Fourier transform of $f(\theta)$ where
$u=[u_{1},\cdots,u_{d}]$ is a vector in $d$-dimension. Given an unitary matrix
$A\in\mathbb{R}^{d\times d}$ such that $A^{\top}A=AA^{\top}=I$ where $I$ is an
identity matrix in $\mathbb{R}^{d\times d}$ and $Ae_{1}=\frac{x}{\|x\|_{2}}$
where $e_{i}$ is vector in $\mathbb{R}^{d}$ with all entries set to $0$ except
$i$th entry which is set to 1, we have

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|u_{1}|^{p+1}}\quad\text{for~{}}p\in[1,2),$

and

$\displaystyle
h(Au)=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{1},u_{2},\cdots,u_{d})\frac{2}{u_{1}^{2}}\quad\text{for~{}}p=2,$

where $\delta$ is the Dirac-delta function.

###### Proof.

We recall that the Fourier transform is given by

$\displaystyle\mathcal{F}f(u)=\int_{\mathbb{R}^{d}}f(\theta)e^{-\mathsf{i}u^{\top}\theta}d\theta.$

Let

$\displaystyle h(u):=\mathcal{F}[|\langle
x,\cdot\rangle|^{p}]=\|x\|_{2}^{p}\mathcal{F}\left[\left|\left\langle\frac{x}{\|x\|_{2}},\cdot\right\rangle\right|^{p}\right].$

We consider now an unitary matrix $A\in\mathbb{R}^{d\times d}$ such that
$A^{\top}A=AA^{\top}=I$ where $I$ is an identity matrix in
$\mathbb{R}^{d\times d}$ and $Ae_{1}=\frac{x}{\|x\|_{2}}$ where $e_{i}$ is
vector in $\mathbb{R}^{d}$ with all entries set to $0$ except $i$th entry
which is set to 1. Now let us compute $h(Au)$.

$\displaystyle h(Au)$
$\displaystyle=\|x\|_{2}^{p}\int_{\mathbb{R}^{d}}\left|\left\langle\frac{x}{\|x\|_{2}},\theta\right\rangle\right|^{p}e^{-\mathsf{i}(Au)^{\top}\theta}~{}d\theta$
$\displaystyle=\|x\|_{2}^{p}\int_{\mathbb{R}^{d}}\left|\left\langle
Ae_{1},\theta\right\rangle\right|^{p}e^{-\mathsf{i}(Au)^{\top}\theta}~{}d\theta.$

In the above integral we substitute, $\beta=A^{\top}\theta$. Hence, when
$p\in[1,2)$, we have

$\displaystyle h(Au)$
$\displaystyle=\|x\|_{2}^{p}\int_{\mathbb{R}^{d}}|\langle
e_{1},\beta\rangle|^{p}e^{-\mathsf{i}u^{\top}\beta}~{}d\beta$
$\displaystyle=\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\int_{-\infty}^{\infty}|\beta_{1}|^{p}~{}e^{-\mathsf{i}u_{1}\beta_{1}}~{}d\beta_{1}$
$\displaystyle=\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\int_{0}^{\infty}\left(e^{-\mathsf{i}u_{1}\beta_{1}}+e^{\mathsf{i}u_{1}\beta_{1}}\right)\beta_{1}^{p}~{}d\beta_{1}$
$\displaystyle=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\Gamma(p+1)\cos\left(\frac{(p+1)\pi}{2}\right)\frac{1}{|u_{1}|^{p+1}}.$

When $p=2$, we have

$\displaystyle h(Au)$
$\displaystyle=\|x\|_{2}^{p}\int_{\mathbb{R}^{d}}|\langle
e_{1},\beta\rangle|^{2}e^{-\mathsf{i}u^{\top}\beta}~{}d\beta$
$\displaystyle=\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{2},\cdots,u_{d})\int_{-\infty}^{\infty}|\beta_{1}|^{2}~{}e^{-\mathsf{i}u_{1}\beta_{1}}~{}d\beta_{1}$
$\displaystyle=2\|x\|_{2}^{p}(2\pi)^{d-1}\delta(u_{1},u_{2},\cdots,u_{d})\frac{2}{u_{1}^{2}}.$
(62)

This completes the proof. ∎

##  Appendix F Further Details on Experiment Settings and Resources

This section contains further details regarding the experiments presented in
the main paper. As the synthetic data experiment setting was fully described
in the text, most of the information below will pertain to the real data
experiments, building up on the information already presented in the main
paper. See the accompanying code regarding the implementation of the
experiments described.

###  F.1 Datasets

The real data experiments involved a supervised learning setting, where images
are classified into a number of predefined class labels. Each model
architecture with given hyperparameters were trained on MNIST [28], CIFAR10,
and CIFAR100 [25] data sets333MNIST and CIFAR10/100 data sets have been shared
under Creative Commons Attribution-Share Alike 3.0 license and MIT License
respectively.. The MNIST data set includes $28\times 28$ black and white
handwritten digits, with digits ranging from $0$ to $9$. The data set in its
original form includes $60000$ training and $10000$ test samples. CIFAR10 and
CIFAR100 are also image classification dataset comprising $32\times 32$ color
images of objects or animals, making up $10$ and $100$ classes respectively.
There are $50000$ training and $10000$ test images in either of these data
sets, and the instances are divided among classes equally. We used the
standard train-test splits in all data sets.

###  F.2 Models

We used three different architectures in our experiments: a fully connected
network with 4 hidden layers (FCN4), another fully connected network with 6
hidden layers (FCN6), and a convolutional neural network (CNN). In both FCN
architectures, all hidden layer widths were $2048$. All architectures featured
ReLU activation functions. Batch normalization, dropout, residual layers, or
any explicit regularization term in the loss function were not used in any
part of the experiments. The architecture we chose for our CNN model closely
follows that of VGG11 model [44], with the significant difference that only a
single linear layer with a softmax output follows the convolutional layers
presented below:

$64,M,128,M,256,256,M,512,512,M,512,512,M.$

Here, integers describe the number of filters for 2-dimensional convolutional
layers - for which the kernel sizes are $3\times 3$. $M$ stands for $2\times
2$ max-pooling operations with a stride value of $2$. This architecture was
slightly modified for the MNIST experiments by removing the first max-pooling
layer due to the smaller dimensions of the MNIST images. The Table 1 describes
the number of different parameters used for each model-dataset combination.

| FCN4 | FCN6 | CNN
---|---|---|---
MNIST | 14,209,024 | 22,597,632 | 9,221,696