# SGD with Clipping is Secretly Estimating the Median Gradient

Fabian Schaipp  Guillaume Garrigos  Umut Şimşekli  Robert M. Gower

###### Abstract

There are several applications of stochastic optimization where one can
benefit from a robust estimate of the gradient. For example, domains such as
distributed learning with corrupted nodes, the presence of large outliers in
the training data, learning under privacy constraints, or even heavy-tailed
noise due to the dynamics of the algorithm itself. Here we study SGD with
robust gradient estimators based on estimating the median. We first consider
computing the median gradient across _samples_ , and show that the resulting
method can converge even under heavy-tailed, state-dependent noise. We then
derive iterative methods based on the stochastic proximal point method for
computing the geometric median and generalizations thereof. Finally we propose
an algorithm estimating the median gradient across _iterations_ , and find
that several well known methods – in particular different forms of clipping –
are particular cases of this framework.

Machine Learning, ICML

figuret

##  1 Introduction

Many problems in machine learning can be represented by

$\min_{\bm{w}\in\mathbb{R}^{d}}\ell(\bm{w}),$

an optimization problem where $\bm{w}$ is a vector of parameters and
$\ell:\mathbb{R}^{d}\to\mathbb{R}$ is a loss function. To tackle this problem,
gradient-based optimization algorithms have been the de facto choice in many
application domains (Nemirovsky & Yudin, 1983; Bottou, 2010; Bottou et al.,
2018), where we are only given access to a _noisy oracle_ of the gradient,
i.e.,

$\bm{g}_{t}=\nabla\ell(\bm{w}_{t})+\bm{\xi}_{t}(\bm{w}_{t}).$

Here $t$ denotes the iterations and
$\bm{\xi}_{t}(\bm{w}_{t})\in\mathbb{R}^{d}$ is noise vector that is sampled
from a distribution that depends on the parameter $\bm{w}_{t}$. With this
noise oracle, stochastic gradient descent (SGD) with learning rate $\eta_{t}$
is given by

$\displaystyle\bm{w}_{t+1}=\bm{w}_{t}-\eta_{t}\bm{g}_{t}.$ (SGD)

While SGD has proven useful in numerous applications, its performance heavily
relies on the ‘quality’ of the noisy oracles $\bm{g}_{t}$, and can tragically
degrade when the noise $\bm{\xi}_{t}$ has significant outliers, exhibits heavy
tails (i.e., $\|\bm{\xi}_{t}\|$ can be very large with non-negligible
probability), or has been corrupted in an adversarial way. We provide two
examples.

  1. (i)

Heavy-tailed gradients. Recent studies have provided theoretical and empirical
evidence that heavy tails can naturally arise in stochastic optimization. From
a theoretical perspective (Gurbuzbalaban et al., 2021; Hodgkinson & Mahoney,
2021; Pavasovic et al., 2023) showed that the parameter-dependent nature of
$\bm{\xi}_{t}$ can result in heavy tails, even in very basic settings such as
linear regression with Gaussian data. On the other hand, it has been
empirically observed that the gradients for training transformer architectures
on language tasks are more heavy-tailed compared to, for example,
convolutional models for image data (Zhang et al., 2020b; Kunstner et al.,
2023). While the presence of heavy-tails in gradients might be beneficial in
certain settings (Simsekli et al., 2020), from an optimization perspective, it
mostly introduces non-trivial challenges, which might even make the algorithm
diverge unless additional care is taken (Zhang et al., 2020b; Gorbunov et al.,
2020).

  2. (ii)

Corrupted nodes. A well-known problem in distributed learning is when some
computation nodes are malicious and can communicate adversarial updates, which
results in inaccurate gradient oracles that can significantly misguide the
optimization procedure. Similar to the heavy-tailed setting, depending on how
much malicious are the nodes, SGD can be impractical unless certain
modifications are made (Mhamdi et al., 2018). We refer to (Karimireddy et al.,
2021) and references therein for an overview of more robust aggregation
techniques.

In order to make the optimization algorithms more robust to such noisy
oracles, many techniques have been developed for different domains under
various conditions. In the context of heavy-tailed gradients, clipping is
often employed (Zhang et al., 2020a; Puchkin et al., 2023; Koloskova et al.,
2023). In distributed learning robust aggregation rules are needed to protect
against corrupted or malicious nodes. (Data et al., 2018; Karimireddy et al.,
2021; Khirirat et al., 2023). This line of research is also closely related to
error feedback in federated learning (Seide et al., 2014; Karimireddy et al.,
2019; Richtárik et al., 2021) and learning under differential privacy
constraints (Khirirat et al., 2023).

#### Robustness of median.

While attracting increasing attention in stochastic optimization thanks to
modern machine learning applications, taming the impact of strong noise has
already been considered in various other domains. Indeed, estimation under
heavy-tails and corruptions has been a long-standing topic in robust
statistics (Huber, 1981). Being a vast research field that has produced
numerous algorithms designed for different tasks, one of the recurring themes
in robust statistics is the use of the _sample median_ and its variations as
opposed to the _sample mean_ (Minsker, 2015), whenever an aggregation of
random variables is needed. For instance, in the presence of heavy tails, the
sample median has been shown to be a significantly more robust estimator of
the true mean, whereas the sample mean can be vulnerable to strong noise
(Lugosi & Mendelson, 2019).

#### Sample median SGD.

Inspired by tools from robust statistics, some notions of median have been
utilized in stochastic optimization as well (Acharya et al., 2022; Yin et al.,
2018). One such approach is to compute a sample median over a sample of
gradients. That is, let

$\bm{g}_{t}^{(i)}=\nabla\ell(\bm{w}_{t})+\bm{\xi}^{(i)}_{t}(\bm{w}_{t})\quad\mbox{for}\quad
i=1,\dots,n,$

where $\\{\bm{\xi}^{(i)}_{t}(\bm{w}_{t})\in\mathbb{R}^{d}\\}_{i=1}^{n}$ is a
set of i.i.d. noise vectors.

Given these sampled gradients, stochastic optimization algorithms typically
aggregate them by computing their sample mean, or average
$\bm{g}_{t}=\frac{1}{n}\sum_{i=1}^{n}\bm{g}_{t}^{(i)}$, and then use this
$\bm{g}_{t}$ in (SGD). Alternatively, to improve robustness, one can aggregate
these gradients $\\{\bm{g}_{t}^{(i)}\\}_{i=1}^{n}$ using a multidimensional
variant of the median. While this approach has paved the way for powerful
optimization algorithms in terms of robustness, it often introduces a
significant computational burden, since a multivariate median needs to be
computed at every iteration (Weiszfeld & Plastria, 2009; Vardi & Zhang, 2000).

Motivation. Even though the clipping-based approaches developed in
optimization theory and the median-based approaches inspired from robust
statistics aim at solving similar tasks, they where designed in different
mathematical frameworks under different conditions. In this study, our main
goal is to bridge the gap between these two seemingly different branches and
to bring a unified theoretical perspective that can shed further light on both
directions by leveraging tools from both optimization and robust statistics.

Contributions. We first revisit median gradient-based optimization and
consider a generic notion of median which include several multivariate medians
that have been extensively used in practice. By exploiting the statistical
properties of the sample median estimator, we show that SGD with
(approximately computed) median gradients indeed yields a very powerful
algorithm, which converges with $1/\sqrt{T}$ rate ($T$ being the number of
iterations) even when the gradient oracles are heavy-tailed with diverging
second-order moments. We further illustrate that, in addition to having
infinite variance, the noise vectors $\bm{\xi}_{t}(\bm{w})$ can even have a
multiplicative dependence on $\bm{w}$, a scenario which is highly challenging
and cannot be directly covered by existing theoretical results, see e.g.,
(Wang et al., 2021; Zhang et al., 2020b; Puchkin et al., 2023).

Encouraged by this result, we then focus on the practical computation of the
medians. To this end, we solve the optimization-based formulation of the
median using Stochastic Proximal Point (SPP) (Asi & Duchi, 2019; Davis &
Drusvyatskiy, 2019; Toulis et al., 2020). As an alternative, we also consider
iterative solvers based on the stochastic subgradient method.

We then depart from previous work by approximately computing the median
gradient _across_ iterations. That is, the median over
$\\{\bm{g}_{t}\\}_{t=0}^{T}.$ We use this online median gradient within the
SGD framework. When using the Stochastic Proximal Point method for computing
the median across iterations, we recover several known clipping-based
optimization algorithms as special cases. In particular, we recover different
forms of clipping (Gorbunov et al., 2020; Zhang et al., 2020a), such as Clip21
for distributed learning (Khirirat et al., 2023), and signed gradient methods
(Riedmiller & Braun, 1993; Bernstein et al., 2018). We also show that (non-
clipped) heavy-ball momentum (Polyak, 1964) is in fact an online mean
estimator.

Our proposed framework allows to unify the different motivations for
previously developed techniques, and illustrates that clipping-based
optimization algorithms essentially runs a “hidden” estimation of the median
across iterations. In particular, the observation that _centered clipping_
(Karimireddy et al., 2021) and Clip21 are an iterative estimator of a
geometric median appears to be new, and might have consequences for its
application in distributed learning.

We finally illustrate our theory on experiments where we compare the
effectiveness of sample median, sample mean, and several online median
estimators. Our results underline that for heavy-tailed noise, using the
sample median is highly effective in contrast to the sample mean which is
unstable and often does not converge. Our experiments also show that our
online median estimates are robust, while being less expensive to compute. We
further compare different clipping techniques for training transformer
architectures on language modeling tasks, and show that they improve upon the
performance of SGD with momentum.

Notation. Throughout the paper, $\|\cdot\|_{p}$ denotes the $\ell_{p}$-norm,
given by
$\|\bm{z}\|_{p}:=\left(\sum_{i=1}^{d}|\bm{z}|_{i}^{p}\right)^{\frac{1}{p}}$
when $p\in[1,\infty)$, and
$\|\bm{z}\|_{\infty}=\max\limits_{i=1,\dots,d}|\bm{z}_{i}|$ when $p=\infty$.
When $p=2$ we recover the usual Euclidean $\ell_{2}$-norm $\|\cdot\|_{2}$,
which we will often simply denote by $\|\cdot\|$. For a matrix $\bm{A}$, we
denote its Frobenius norm by
$\|\bm{A}\|_{F}:=\sqrt{\sum_{i,j}\bm{A}_{ij}^{2}}$.

##  2 Preliminaries

###  2.1 Heavy-tailed Random Variables

A random variable $X$ is said to be _heavy-tailed_ if the tails of its
distribution decay slower than the tails of an exponential distribution. One
important example, is the family of symmetric $\alpha$-stable distributions
(Nolan, 2020, Def. 1.4), denoted $\mathcal{P}_{\alpha,\sigma}$: it is
parameterized by a stability parameter $\alpha\in(0,2]$ and a scale parameter
$\sigma>0$. For $\alpha=1$, it is the Cauchy distribution, while for
$\alpha=2$ it is a Gaussian. The parameter $\sigma$ controls the spread of the
distribution. When $\sigma=1$, we call the distribution _standard_ and denote
it by $\mathcal{P}_{\alpha}:=\mathcal{P}_{\alpha,1}$. The parameter $\alpha$
on the other hand controls the heaviness of the tails: for $\alpha<2$ and
$X\sim\mathcal{P}_{\alpha}$ has infinite variance, whereas for $\alpha\leq 1$
even the mean of $X$ is infinite. However, for any $\alpha$, the median is
equal to $0$.

###  2.2 Mean Estimation and (Sample) Median

Estimating the mean of a distribution, given a finite number of samples, is
one of the core problems studied in statistics. When the distribution in
question has heavy-tails, or has significant outliers, the sample mean turns
out to be a poor estimator. We refer to (Lugosi & Mendelson, 2019) for a
survey of alternative techniques. In this section, we will introduce basic
definitions for median.

In this section we denote by $z$ a real-value random value, and with
$\bm{z}\in\mathbb{R}^{d}$ a random vector. We always denote the expectation
with respect to the distribution of $z$ (or $\bm{z}$ respectively) by
$\mathbb{E}$.

#### Multivariate median.

The one-dimensional median is defined by

$\displaystyle\texttt{median}{}(z)\in\underset{m\in\mathbb{R}}{\operatorname{argmin}}\;\mathbb{E}\left[|m-z|\right].$
(1)

In contrast to the mean, the median is always defined (though it may not be
unique) (Cramér, 2016).

Since we are interested in approximating the median gradient, we need a notion
of median that generalizes to random vectors $\bm{z}\in\mathbb{R}^{d}$. For
this, we can define the $\ell_{p}$-median:

$\texttt{median}_{\ell_{p}}(\bm{z})\in\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathbb{E}\left[\|\bm{m}-\bm{z}\|_{p}\right].$
(2)

This recovers several existing notions of a multivariate median: for $p=1$ we
obtain a _componentwise median_ (or $\ell_{1}$-median) and for $p=2$ we obtain
the _geometric median_ (or $\ell_{2}$-median) (Fréchet, 1948; Weiszfeld &
Plastria, 2009; Minsker, 2015; Cohen et al., 2016). Clearly, for $d=1$ all
notions conincide with the definition in (1).

#### Sample median.

Given $n$ realizations $\bm{z}^{(1)},\ldots,\bm{z}^{(n)}$ of $\bm{z}$, the
sample median is defined as

$\texttt{median}_{\ell_{p},i\in[n]}(\bm{z}^{(i)})\in\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\frac{1}{n}\sum_{i=1}^{n}\|\bm{m}-\bm{z}^{(i)}\|_{p}.$
(3)

We remark that (3) is a special case of (2), by setting the distribution of
$\bm{z}$ as a sum of Dirac measures
$\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{z}_{i}}$.

It is generally accepted that the sample median is a more robust estimator
than the sample mean (Huber, 1981; Hampel et al., 2005). This is nicely
illustrated by the notion of a breakdown point (Donoho & Huber, 1983; Lopuhaa
& Rousseeuw, 1991), which is the smallest fraction of a sample that, if
arbitrarily corrupted, can arbitrarily change the value of the estimator. The
breakdown point of the sample median is roughly $\tfrac{1}{2}$, meaning that
at least half of the samples need to be outliers for the median to diverge.
The sample mean however can be arbitrarily corrupted by a single outlier, thus
its breakdown point is one over $\tfrac{1}{n}$.

If the samples $\bm{z}^{(i)}$ are obtained by taking block-wise means over
observations, (3) is called the median-of-means estimator (Nemirovsky & Yudin,
1983; Lugosi & Mendelson, 2019), which is known to be a robust mean estimator
(Lugosi & Mendelson, 2019).

#### Relation to M-estimation.

More generally, for a function $\mathcal{D}:\mathbb{R}^{d}\to\mathbb{R}$ we
will later consider the problem

$\displaystyle\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathbb{E}\left[\mathcal{D}(\bm{m}-\bm{z})\right].$
(4)

We assume that $\mathcal{D}$ is a closed, proper, convex function and that (4)
admits a solution (see Appendix B for a discussion). For
$\mathcal{D}=\tfrac{1}{2}\|\cdot\|_{2}^{2}$, the solution to (4) is clearly
the mean $\mathbb{E}\left[\bm{z}\right]$, and for $\mathcal{D}=\|\cdot\|_{2}$,
the solution is the geometric median as defined in (2). Problem (4) can also
be seen as a special case of _M-estimators_ (Huber, 1981, Sec. 3.2).

The choice of $\mathcal{D}$ determines properties of the associated estimator:
if the function $\mathcal{D}$ heavily penalizes large values, like the
quadratic $\mathcal{D}=\tfrac{1}{2}\|\cdot\|_{2}^{2}$, then the solution to
(4) will be sensitive to outliers. In contrast $\mathcal{D}=\|\cdot\|_{2}$ and
$\mathcal{D}=\|\cdot\|_{1}$ grow only linearly, and thus (4) will be less
sensitive to outliers. We return to problem (4) later in Section 4.

For SGD, its update direction is a sampled gradient $\bm{g}_{t}$ which aims to
approximate the full gradient $\nabla\ell(\bm{w}_{t})$. This section motivates
that in the presence of large outliers it might be necessary to use more
robust gradient estimators as update direction. A natural alternative is to
use the sample median together with SGD (Alistarh et al., 2018; Acharya et
al., 2022), as we do next.

##  3 Warm-up: Sample Median Gradient

In this section, we showcase that SGD with an (approximate) sample median is
robust to heavy-tailed noise, whereas SGD using the sample mean can fail.
Consider the following setup: at every iteration $t\geq 0$, we have access to
$n$ noisy oracle gradients
$\big{(}\bm{g}_{t}^{(1)},\ldots,\bm{g}_{t}^{(n)}\big{)}$, and we apply

$\displaystyle\bm{w}_{t+1}=\bm{w}_{t}-\eta(\bm{m}_{t}+\bm{e}_{t}),$ (5)

where $\bm{m}_{t}=\texttt{median}_{\ell_{p},i\in[n]}(\bm{g}_{t}^{(i)})$, and
$\bm{e}_{t}\in\mathbb{R}^{d}$ represents an error we may commit in computing
the sample median.

######  Assumption 3.1.

There exists $\sigma_{1},\sigma_{2},\delta_{1},\delta_{2}\geq 0$ for every
$\bm{w}\in\mathbb{R}^{d}$, such that the sample median
$\bm{m}=\texttt{median}_{\ell_{p},i\in[n]}(\bm{g}^{(i)})$ verifies

$\displaystyle\begin{split}&\|\mathbb{E}\left[\bm{m}\right]-\nabla\ell(\bm{w})\|^{2}\leq\delta_{1}^{2}+\delta_{2}\|\nabla\ell(\bm{w})\|^{2},\\\
&\mathbb{E}\left[\|\bm{m}-\mathbb{E}\left[\bm{m}\right]\|^{2}\right]\leq\sigma_{1}^{2}+\sigma_{2}\|\nabla\ell(\bm{w})\|^{2}.\end{split}$
(6)

######  Assumption 3.2.

There exists $\varepsilon>0$ such that, for every $t\geq 0$, the error in
computing the sample median has finite variance:
$\mathbb{E}\left[\|\bm{e}_{t}\|^{2}\ |\
\bm{g}_{t}^{(1)},\dots,\bm{g}_{t}^{(n)}\right]\leq\varepsilon^{2}$.

Based on these two assumptions, we establish the following complexity results
for (5).

######  Proposition 3.3.

Let $\ell:\mathbb{R}^{d}\to\mathbb{R}$ be $L$-smooth and let Assumption 3.1
hold with $\delta_{2}\leq\frac{1}{8}$. Take
$\eta\leq\frac{1}{8L(1+2\sigma_{2}+2\delta_{2})}$ and consider $\bm{w}_{t}$ a
sequence generated by (5) where $\bm{e}_{t}$ verifies Assumption 3.2. Then,
for every $T\geq 1$ it holds

$\min_{0\leq t\leq
T-1}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]=\mathcal{O}\left(\tfrac{1}{\eta
T}+\eta+\delta_{1}^{2}+\varepsilon^{2}\right).$

######  Proposition 3.4.

Let $\ell:\mathbb{R}^{d}\to\mathbb{R}$ be convex and $G$-Lipschitz, with
${\rm{argmin}}~{}\ell\neq\emptyset$. Let Assumptions 3.1 and Assumption 3.2
hold. Without loss of generality assume that $\delta_{2}=\sigma_{2}=0$. Then,
for $\eta=\frac{1}{(\delta+\varepsilon)T+\sqrt{T}}$ and every $T\geq 1,$ the
iterates $\bm{w}_{t}$ from (5) satisfy

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]=\mathcal{O}\left(\frac{1}{\sqrt{T}}+\delta_{1}+\varepsilon\right),$

where $\bar{\bm{w}}_{T}:=\tfrac{1}{T}\sum_{t=0}^{T-1}\theta^{t}\bm{w}_{t}$
with $\theta=\tfrac{T}{T+2}$.

The proofs are provided in Appendix A. Let us now discuss our main Assumption
3.1. Assume for simplicity that $\delta_{2}=\sigma_{2}=0$, that
$\bm{e}_{t}=0$, and consider that
$\bm{g}_{t}^{(i)}=\nabla\ell(\bm{w}_{t})+\bm{\xi}_{t}^{(i)}$, where
$\bm{\xi}_{t}^{(i)}$ is a noise term. Then, Assumption 3.1 becomes111See Lemma
A.2 in the appendix for a proof.

$\displaystyle\mathbb{E}\left[\|\texttt{median}_{\ell_{p},i\in[n]}(\bm{\xi}_{t}^{(i)})\|^{2}\right]\leq\nu^{2}<+\infty.$
(7)

In other words, we need the second moment of the sample median of
$\bm{\xi}_{t}^{(i)}$ to be uniformly bounded. We want to stress that our
choice of using the median is crucial when the $\bm{\xi}_{t}^{(i)}$ come from
a heavy-tailed distribution, such as a (univariate) standard $\alpha$-stable
distribution, with $\alpha\in(1,2)$. Indeed, in this case we know that the
variance of the sample median is _finite_ for sufficiently large sample size
$n$ (see (Bickel, 1967, Theorem 2.2) and (Nolan, 2020, Theorem 1.2))222We know
more: The sample median of a parent distribution with median $\xi$ is
asymptotically normal with location $\xi$ (Chu & Hotelling, 1955). For Cauchy
distributions ($\alpha=1$), the sample median has a finite variance for $n\geq
5$ (Chu & Hotelling, 1955, Remark 1), while first and second moment of the
sample mean are undefined.. On the contrary, if we were taking $\bm{m}_{t}$ as
a more classical sample _mean_ we would be in trouble, because the variance of
the sample mean would be _infinite_ 333This is due to the fact that, by
definition of $\alpha$-stability, the sample mean of i.i.d sample follows the
same $\alpha$-stable distribution (up to affine translation). everywhere. See
Fig. 1 for an illustration of this phenomenon.

This simple discussion illustrates that Assumption 3.1 can be satisfied under
mild conditions on the noise. We exploit this in the next corollary, which
shows that the $\ell_{1}$-sample median can even accommodate state-dependent
noise. Let us define first our model for such noise.

######  Assumption 3.5.

For every $\bm{w}\in\mathbb{R}^{d}$, we have that
$\bm{g}=\nabla\ell(\bm{w})+\bm{\Sigma}(\bm{w})\bm{\zeta}$, where
$\bm{\Sigma}(\bm{w})\in\mathbb{R}^{d\times d}$ is a matrix and $\bm{\zeta}$ is
a random vector whose components are drawn i.i.d. from a standard symmetric
$\alpha$-stable distribution with $\alpha~{}>~{}1$. There exists constants
$c_{1},c_{2}>0$ such that

$\|\bm{\Sigma}(\bm{w})\|_{F}^{2}\leq
c_{1}+c_{2}\|\nabla\ell(\bm{w})\|^{2}\quad\mbox{for all
}\bm{w}\in\mathbb{R}^{d}.$

As opposed to the following assumption that has been commonly made in the
recent literature (see e.g., (Zhang et al., 2020b; Puchkin et al., 2023;
Sadiev et al., 2023))
$\mathbb{E}\left[\|\bm{g}-\nabla\ell(\bm{w})\|^{p}\right]\leq\sigma^{2}$ for
some $p<2$ where $\sigma^{2}>0$ cannot depend on $\bm{w}$, Assumption 3.5 can
accommodate much stronger noise. To illustrate this, consider $d=1$ and
$\ell(\bm{w})=\bm{w}^{2}/2$. Assumption 3.5 allows the noise model
$\bm{\Sigma}(\bm{w})\bm{\zeta}$ with
$\bm{\Sigma}(\bm{w})=\sqrt{c_{1}+c_{2}\bm{w}^{2}}$, which results in strong
noise, since $\bm{\zeta}$ has infinite variance if $\alpha<2$, and this is
amplified by $\bm{\Sigma}(\bm{w})$ when $\bm{w}$ is large. Such strong noise
makes vanilla SGD diverge, as we will illustrate in our experiments.

######  Corollary 3.6.

Let $n\geq 3$ be odd and $p=1$. Consider the iterates from (5) with
$\bm{e}_{t}=0$, and suppose that
$\bm{g}_{t}^{(i)}=\nabla\ell(\bm{w}_{t})+\bm{\Sigma}(\bm{w}_{t})\bm{\zeta}_{t}^{(i)}$
verifies Assumption 3.5.

  1. (i)

If $\ell$ is smooth and $\eta=1/\sqrt{T}$ then we have that $\min_{0\leq t\leq
T-1}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]=\mathcal{O}(1/\sqrt{T})$.

  2. (ii)

If $\ell$ is convex and Lipschitz and $\eta=1/\sqrt{T}$ then we have that
$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]=\mathcal{O}(1/\sqrt{T})$,
where $\bar{\bm{w}}_{T}$ is defined as in Proposition 3.4.

This result illustrates how powerful the sample median can be when used in
gradient-based optimization. However, the method (5) has clear drawbacks: (i)
it requires to draw multiple samples per iteration (multiple calls of the
oracle), and (ii) we need to compute the sample median up to some tolerance.
For instance, the geometric median can not be computed in closed form and
there exists a vast body of literature on how to approximate it (Cohen et al.,
2016).

We tackle (ii) in the next section by developing an efficient and inexpensive
iterative method for approximately solving the general problem (4) which
includes the sample median as a special case. We tackle issue (i) later in
Section 5.

Figure 1: First and second moment of sample median/mean for standard
($\alpha=1.1$)-stable distribution. We approximate $\mathbb{E}$ by running
$1000$ trials. Even for small sample size, the sample median has small
variance (cf. (7)), in contrast to the sample mean.

##  4 Iterative Methods for Median Estimation

We now focus on the following gradient estimation problem: let $\mathcal{G}$
denote the (fixed) distribution over noisy gradients $\bm{g}$, and consider
the problem of approximating

$\displaystyle\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathbb{E}_{\bm{g}\sim\mathcal{G}}\left[\mathcal{D}(\bm{m}-\bm{g})\right].$
(8)

In the previous section, we used the sample median. It follows that solving
(8) with $\mathcal{D}=\|\cdot\|_{p}$, and with $\mathcal{G}$ being the
empirical distribution over a sample $(\bm{g}^{(1)},\ldots,\bm{g}^{(n)})$, the
solution to (8) is the sample median (3).

For solving (8), we consider only stochastic iterative methods, that can be
applied for any distribution $\mathcal{G}$. The most commonly used such method
is SGD, which has been used for computing the median (Cardot et al., 2011).
But here we advocate for a different method: The stochastic proximal point
method (SPP) (Asi & Duchi, 2019; Davis & Drusvyatskiy, 2019). In cases where
SPP has a closed form update, it is often preferred to SGD since is easier to
tune (Asi & Duchi, 2019; Milzarek et al., 2022). Fortunately, as we show next,
the SPP method enjoys closed-form updates when
$\mathcal{D}\in\\{\|\cdot\|_{1},\|\cdot\|_{2},\tfrac{1}{2}\|\cdot\|^{2}\\}$.

The update of SPP with stepsize $\tau>0$ applied to (8) can be written as (see
details in Appendix B.3)

$\displaystyle\bm{m}_{t+1}\;=\;\bm{g}_{t}+\mathrm{prox}_{\tau\mathcal{D}}(\bm{m}_{t}-\bm{g}_{t}),$
(9)

where we denote by $\mathrm{prox}_{\mathcal{D}}$ the proximal operator of
$\mathcal{D}$,

$\mathrm{prox}_{\mathcal{D}}(\bm{m}):=\underset{\bm{y}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathcal{D}(\bm{y})+\frac{1}{2}\|\bm{y}-\bm{m}\|^{2}.$

Because we are interested in the $\ell_{p}$-median, we first characterize the
SPP update (9) for when $\mathcal{D}$ is the $\ell_{p}$-norm, and show it
simply require to project the sampled gradient onto a certain ball of radius
$\tau$ centered at $\bm{m}_{t}$.

######  Proposition 4.1.

Let $\mathcal{D}=\|\cdot\|_{p}$ be the $\ell_{p}$-norm for $p\in[1,\infty].$
Let $q\in[1,\infty]$ such that $\frac{1}{p}+\frac{1}{q}=1$, and let
$\mathbb{B}_{q}(\bm{m},\tau)\;:=\;\left\\{\bm{y}\;:\;\|\bm{m}-\bm{y}\|_{q}\leq\tau\right\\}.$
The SPP update (9) is given by

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{m}_{t}\;+\;\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},\tau)}\big{(}\bm{g}_{t}-\bm{m}_{t}\big{)}$
(10)
$\displaystyle=\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m}_{t},\tau)}\big{(}\bm{g}_{t}\big{)}.$
(11)

Figure 2: One SPP update (11) for $\mathcal{D}=\|\cdot\|_{1}$ (blue) or
$\mathcal{D}=\|\cdot\|_{2}$ (green) amounts to project onto a ball of the dual
norm, centered at the current iterate $\bm{m}_{t}$ and of radius $\tau>0$. The
thick lines represent the projections of $\bm{g}_{t}$ for all possible values
of $\tau$.

Using Proposition 4.1 we can now develop closed form updates of the SPP method
for when $\mathcal{D}$ is the $\ell_{1}$\- or the $\ell_{2}$-norm, with both
methods being related to clipping.

#### Vectorwise clipping.

If we choose $\mathcal{D}=\|\cdot\|_{2}$, SPP performs projections onto the
$\ell^{2}$-ball which is also known as clipping (Khirirat et al., 2023;
Koloskova et al., 2023).

######  Corollary 4.2.

For $\mathcal{D}=\|\cdot\|_{2}$ update (10) is given by

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{m}_{t}+\mathrm{clip}_{\tau,2}(\bm{g}_{t}-\bm{m}_{t}),$ (12)

where
$\mathrm{clip}_{\tau,2}(\bm{v}):=\frac{\tau}{\max\\{\tau,\|\bm{v}\|_{2}\\}}\bm{v}$.

This corollary shows that the SPP method for computing the $\ell_{2}$-median
can be written as $\bm{m}_{t+1}=\beta_{t}\bm{m}_{t}+(1-\beta_{t})\bm{g}_{t}$,
with
$\beta_{t}:=1-\tfrac{\tau}{\max\\{\tau,\|\bm{g}_{t}-\bm{m}_{t}\|_{2}\\}}$. In
this sense, (12) is momentum where the momentum coefficient depends on how
_far_ $\bm{g}_{t}$ is from the previous estimate $\bm{m}_{t}$. This allows the
iterates to be more robust to outliers, as samples which are too far from
$\bm{m}_{t}$ will be down weighted (see also Fig. 2).

#### Componentwise clipping.

If we choose $\mathcal{D}=\|\cdot\|_{1}$, SPP performs projections onto the
$\ell^{\infty}$-ball, also known as componentwise clipping. Here as well,
clipping protects against large individual entries by shrinking them
independently.

######  Corollary 4.3.

For $\mathcal{D}=\|\cdot\|_{1}$ update (10) is given by

$\displaystyle\bm{m}_{t+1}=\bm{m}_{t}+\mathrm{clip}_{\tau,1}(\bm{g}_{t}-\bm{m}_{t}),$
(13)

where
$\mathrm{clip}_{\tau,1}(\bm{v}):=\left(\min\\{\max\\{\bm{v}_{i},-\tau\\},\tau\\}\right)_{i=1}^{d}$.

####  SPP and momentum.

It is easy to see that the solution to problem (8) with
$\mathcal{D}=\tfrac{1}{2}\|\cdot\|^{2}$ is the mean
$\mathbb{E}\left[\bm{g}\right]$. Computing (9) for this choice of
$\mathcal{D}$, we obtain

$\displaystyle\bm{m}_{t+1}=(1-\tfrac{\tau}{1+\tau})\bm{m}_{t}+\tfrac{\tau}{1+\tau}\bm{g}_{t}.$
(14)

Interestingly, this is the (heavy-ball) momentum gradient estimator (Polyak,
1964; Sutskever et al., 2013), where the momentum coefficient is given by
$\beta=1/(1+\tau)\in[0,1).$

#### Rates for SPP

Now that we have presented how to compute the iterates (9) in practice, let us
state its convergence properties. Since all these methods are instantiates of
the SPP method, we can use standard results such as (Davis & Drusvyatskiy,
2019, Thm. 4.4).

######  Proposition 4.4.

Let $\mathcal{D}$ be a norm, let
$\varphi(\bm{m}):=\mathbb{E}_{\bm{g}\sim\mathcal{G}}\left[\mathcal{D}(\bm{m}-\bm{g})\right]$,
and let
$\bm{m}_{*}\in\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\varphi(\bm{m})$.
If $\bm{m}_{t}$ is generated by (9) then for
$\bar{\bm{m}}_{T}=\frac{1}{T}\sum_{t=1}^{T}\bm{m}_{t}$ we have

$\mathbb{E}\left[\varphi(\bar{\bm{m}}_{T})-\varphi(\bm{m}_{*})\right]=O\left(\frac{\|\bm{m}_{0}-\bm{m}_{*}\|^{2}}{\tau
T}+\tau\right).$

In Appendix B we give further examples for $\mathcal{D}$, such as the Huber
loss. We also compare SPP with SGD for solving (8), showing that sign-SGD is
related to $\ell_{1}$-median estimation.

##  5 Gradient Estimator in the Wild

Our final objective is to use robust gradient estimators within a SGD-type
method. Because at each iteration the weights $\bm{w}_{t}$ are updated, the
distribution of the gradients (denoted by $\mathcal{G}$ in the previous
section) also changes at each iteration of SGD. Our strategy is to interweave
updates in the weights $\bm{w}_{t}$ with SPP updates in the gradient
estimators $\bm{m}_{t}.$ That is, we sample _one_ stochastic gradient
$\bm{g}_{t}$ per iteration and set

$\displaystyle\begin{array}[]{llr}\bm{m}_{t+1}&=\quad\bm{g}_{t}+\mathrm{prox}_{\tau\mathcal{D}}(\bm{m}_{t}-\bm{g}_{t}),\\\
\bm{w}_{t+1}&=\quad\bm{w}_{t}-\eta_{t}\bm{m}_{t+1}.\\\ \end{array}$ (17)

Formally, $\bm{m}_{t}$ is now an estimate of the median of the gradients
_across_ iterations, which is in contrast to (5), where $\bm{m}_{t}$ was an
approximation of the sample median at fixed $\bm{w}_{t}$. With the results
from Section 4, we immediately obtain the following special cases (here
$p\in\\{1,2\\}$):

(i) if $\mathcal{D}=\frac{1}{2}\|\cdot\|^{2}$ :
$\begin{cases}\bm{m}_{t+1}=(1-\tfrac{\tau}{1+\tau})\bm{m}_{t}+\tfrac{\tau}{1+\tau}\bm{g}_{t}\\\
\bm{w}_{t+1}=\bm{w}_{t}-\eta_{t}\bm{m}_{t+1},\end{cases}$

(ii) if $\mathcal{D}=\|\cdot\|_{p}$ :
$\begin{cases}\bm{m}_{t+1}=\bm{m}_{t}+{\rm
clip}_{p,\tau}(\bm{g}_{t}-\bm{m}_{t})\\\
\bm{w}_{t+1}=\bm{w}_{t}-\eta_{t}\bm{m}_{t+1}.\end{cases}$

Equation (i) is SGD with momentum (14); notice how the SPP step size $\tau$
defines the momentum coefficient. Instead, (ii) corresponds to clipping, and
here the SPP step size determines the clipping threshold. In other words, SGD
with momentum is to online mean estimation what clipping is to online median
estimation.

A slightly different approach would be to restart $\bm{m}_{t}=\bm{0}$ in every
iteration (_cold start_). In this case, (i) reduces to SGD with no momentum,
while in (ii) the clipping operator is applied to $\bm{g}_{t}$, instead of
$\bm{g}_{t}-\bm{m}_{t}$, which is also a common technique in practice (Zhang
et al., 2020a). We discuss this relation below in more detail.

#### Connection to Error Feedback.

Methods such as (17) are also used in the distributed learning setting to
improve communication and protect against corrupted nodes (Karimireddy et al.,
2021). These methods are often called Error Feedback (EF) methods (Seide et
al., 2014; Richtárik et al., 2021) and generally involve an update of the form
$\bm{m}_{t+1}=\bm{m}_{t}+\mathcal{C}_{t}(\bm{g}_{t}-\bm{m}_{t})$ where
$\mathcal{C}_{t}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ is a compression
operator. Projections onto balls can be seen as a compression, since one can
encode the resulting vector in fewer bits (Safaryan et al., 2021). In
particular, centered clipping from (Karimireddy et al., 2021), and the Clip21
method (Khirirat et al., 2023) are using update (12) in the distributed
setting, and are thus secretly estimating the geometric median of the
gradients returned across all nodes.

####  SGD with standard clipping.

Even in the non-distributed setting, a standard technique to avoid training
instabilities is to directly clip the gradient $\bm{g}_{t}$ (Zhang et al.,
2020a). For momentum coefficient $\beta\in[0,1)$ and $c>0$, let

$\displaystyle\begin{split}\bm{m}_{t+1}&=\beta\bm{m}_{t}+(1-\beta)\min\\{1,\tfrac{c}{\|\bm{g}_{t}\|}\\}\bm{g}_{t},\\\
\bm{w}_{t+1}&=\bm{w}_{t}-\eta_{t}\bm{m}_{t+1}.\end{split}$ (18)

We will refer to (18) as clipped-SGD. The method has been studied extensively,
often for the case $\beta=0$, see (Zhang et al., 2020a; Gorbunov et al., 2020;
Koloskova et al., 2023; Puchkin et al., 2023) and references therein.

We present two motivations for using clipped-SGD. If $\beta=0$ then (18) is
the same as (12) with cold start (as explained above). When $\beta\neq 0$,
clipped-SGD can also be viewed as a SPP method for computing a robust mean
estimator given by

$\displaystyle\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathbb{E}_{\bm{g}\sim\mathcal{G}}\left[\tfrac{1}{2}\|\bm{m}-\min\\{1,\tfrac{c}{\|\bm{g}\|}\\}\bm{g}\|^{2}\right].$
(19)

The solution to (19) is
$\mathbb{E}_{\bm{g}\sim\mathcal{G}}\left[\min\\{1,\tfrac{c}{\|\bm{g}\|}\\}\bm{g}\right]$,
which is known as the Catoni-Giulini mean estimator (Catoni & Giulini, 2018;
Lugosi & Mendelson, 2019). One step of SPP (9) applied to (19) gives (18) with
$\beta=\tfrac{1}{1+\tau}$.

##  6 Experiments

We compare the iterative gradient estimation techniques for two experimental
setups, namely (i) estimating a vector/gradient from a fixed distribution, and
(ii) methods of the form (17), that is, iterative gradient estimation while
simultaneously learning the weights.

We study both settings for synthetic data where we can control the level of
heavy-tailedness, and for language modelling tasks with transformer models. In
the following, we use VClip to refer to vectorwise clipping (12) and CClip for
componentwise clipping (13).

We consider the same three language modeling tasks as studied in (Kunstner et
al., 2023): an encoder-only transformer for the PTB, a Transformer-XL model
for the WikiText-2, and fine-tuning a DistillBERT for question-answering on
the SQuAD dataset. We refer to Appendix C for details.

###  6.1 Gradient Estimation with Fixed Weights

#### Synthetic data.

In the first experiment, we verify the hypothesis that the choice of
$\mathcal{D}$ matters for problem (4) – in particular, we show that estimating
the median is more stable when the noise distribution is heavy-tailed.

We consider the problem of estimating a fixed vector $\bm{\hat{g}}$ from the
oracle $\bm{g}\sim\bm{\hat{g}}+\bm{\xi}$ with $\bm{\xi}\sim\mathcal{P}$, under
varying degree of heavy-tailedness of $\mathcal{P}$. For this purpose, we
choose the standard, unskewed $\alpha$-stable distribution
$\mathcal{P}_{\alpha}$.

_Setup._ We generate $\bm{\hat{g}}\in\mathbb{R}^{d}$ with $d=10$ where each
component is generated i.i.d standard Gaussian. In each iteration, a sample
$\bm{g}_{t}$ is generated as follows: each coordinate of $(\bm{g}_{t})_{i}$ is
(independently) sampled from $\mathcal{P}_{\alpha}$ with location
$(\hat{\bm{g}})_{i}$ and varying values for $0<\alpha\leq 2$. Importantly, the
median of $\mathcal{P}_{\alpha}$ is equal to the location (i.e.
$\bm{\hat{g}}$) for all values of $\alpha$. Recall that for $\alpha\leq 1$ the
mean of $\mathcal{P}_{\alpha}$ is not defined, and otherwise equal to the
median.

We run the SPP method (9) for
$\mathcal{D}\in\\{\tfrac{1}{2}\|\cdot\|^{2},\|\cdot\|_{2},\|\cdot\|_{1}\\}$,
with a fixed step size $\tau=0.01$ and track the relative $\ell_{2}$-error
$\frac{\|\bm{m}_{t}-\bm{\hat{g}}\|_{2}}{\|\bm{\hat{g}}\|_{2}}$. The
corresponding methods are called momentum, VClip, and CClip. For each method
and distribution, we run $1000$ iterations and $50$ different seeds.

_Discussion._ From Fig. 3 we find that as the noise becomes more heavy-tailed
(as $\alpha$ decreases), momentum fails to produce accurate estimates of
$\bm{\hat{g}}$. On the other hand, both VClip and CClip are robust to heavy
tails, as expected, as we showed that VClip and CClip are estimating the
median.

Figure 3: $\tau=0.01$ Left: Final error for varying values of $\alpha$ (from
left to right, distributions are more heavy-tailed). Shaded area marks minimal
and maximal value over the $50$ independent runs. Right: Convergence plot for
all methods for $\tfrac{1}{\alpha}\in[0.5,1.5]$ (higher value of
$\tfrac{1}{\alpha}$ corresponds to heavier tails).

#### Language modelling.

We also run momentum, VClip and CClip in order to estimate the full gradient
for the language modeling tasks, where the weights are _fixed at
initialization_. The setup is identical to the one described in (Kunstner et
al., 2023), in particular Figure 1 therein. We present the results in Section
C.1.

###  6.2 Training with Online Gradient Estimation

#### Least squares with heavy-tailed noise.

We consider the noise-perturbed least-squares problem
$\min_{\bm{w}\in\mathbb{R}^{d}}\tfrac{1}{2}\|\bm{w}\|^{2}$ and use a gradient
oracle given by $\bm{g}\sim\bm{w}+\bm{\xi}\in\mathbb{R}^{d}$, where
$\bm{\xi}\sim\mathcal{P}$ is the noise term. We consider two setups:

  1. (S1)

The $d$ components of $\bm{\xi}$ are i.i.d. distributed according to a
standard $\alpha$-stable distribution with location zero and $\alpha=1.1$.
Note that $\mathbb{E}_{\bm{\xi}\sim\mathcal{P}}\left[\bm{\xi}\right]=\bm{0}$.

  2. (S2)

Same as (S1), but the noise is state-dependent, given by
$\bm{\xi}_{t}=\sigma\cdot\bm{\xi}$ where $\sigma=\sqrt{1+\|\bm{w}_{t}\|^{2}}$.

  3. (S3)

$\bm{\xi}$ follows an elliptically contoured $\alpha$-stable distribution
(Nolan, 2013) with location zero and $\alpha=1.1$. In particular, the
components of $\bm{\xi}$ are not independent.

We run six different methods: on the one hand, the method depicted in Section
3, (5), where in each iteration we draw $n=5$ samples. We run three variations
of (5), namely using the $\ell_{1}$\- and $\ell_{2}$-sample median as well as
the sample mean for $\bm{m}_{t}$. The $\ell_{2}$-median is approximately
computed using the method proposed in (Vardi & Zhang, 2000).

(a) Setting (S1)

(b) Setting (S2)

(c) Setting (S3)

Figure 4: Least squares with heavy-tailed noise that is independent over
components (left), state-dependent (middle), and where components of noise are
dependent (right). All methods in grey use $n=5$ samples per iteration, all
others only one. Note that for (S2), the sample mean immediately diverges, and
hence does not appear in the plot. Shaded area depicts minimum and maximum
value over 50 repetitions.

On the other hand, we run our online median estimator methods (17), with
$\mathcal{D}=\|\cdot\|_{2}$ and $\mathcal{D}=\|\cdot\|_{1}$, that is VClip and
CClip. Importantly, these methods only receive one sample per iteration, and
they do not involve any median computations (in contrast to $\ell_{1}$\- and
$\ell_{2}$-sample median). We also run clipped-SGD (cf. (18)) with $\beta=0.9$
and $c=50$.444This choice is a tradeoff between slow convergence (large $c$)
and high instability (small $c$). We run all methods with a learning rate
$\eta_{t}=0.01$. For VClip and CClip we set $\tau=1$. We always run 50
repetitions and report averaged metrics.

_Discussion._ Fig. 4 shows that both $\ell_{1}$\- and $\ell_{2}$-sample median
attain the smallest objective, while the sample mean does not converge as
predicted by our theory (see discussion after Proposition 3.3). From the
online methods, CClip performs best in the scenario where the noise is
componentwise independent, while VClip is better in the other case. We find
that both VClip and CClip reach much smaller loss values than clipping the
stochastic gradient directly, and thus clipping only the increment is better
suited here. We also find that the $\ell_{1}$-sample median can deal with
state-dependent noise (S2) as predicted by our theory.

#### Language modeling.

We now consider (17), that is training while simultaneously estimating the
gradient, for the three language modelling tasks. We compare SGD-M, VClip and
CClip, and clipped-SGD. We also add Adam as a baseline, but focus on the
comparison among the SGD-type methods. For all methods, we tune the learning
rate on a $\log_{10}$-scaled grid (reported in Table 1), displaying the one
that attained smallest final training loss (averaged over $5$ seeds).

We choose a standard momentum/clipping parameters for all tasks (without
tuning): we set $\beta=0.9$ for SGD-M, $\tau=0.1$ for V/CCLip, and
$\beta=0.9,~{}c=1$ for clipped-SGD.

_Discussion._ Fig. 5 shows that VClip improves over SGD-M for PTB and
WikiText-2, and is on par for SQuAD. CClip performs always worse, which might
be due to the noise being not componentwise independent. However, VClip can
not close the gap to Adam, and also performs worse than clipped-SGD for the
WikiText-2 dataset.

(a) PTB

(b) WikiText-2

(c) SQuAD

Figure 5: Training loss for each method, with tuned learning rate. Shaded area
depicts minimum and maximum value over 5 seeds.

(a) PTB

(b) WikiText-2

(c) SQuAD

Figure 6: Validation score (measured as perplexity, or F1 score) for each
method, with tuned learning rate. Shaded area depicts minimum and maximum
value over 5 seeds.

##  7 Conclusion

We provide an analysis of SGD using the median of sampled gradients, and
showed convergence guarantees under potentially heavy-tailed and state-
dependent noise. In contrast, under the same assumption, SGD with the sample
mean does not converge. We further provide a new approach to compute medians
using the stochastic proximal point method. Our final algorithm proposes to
use online estimates of the median gradient, requiring only one sample per
iteration. We observe that this general framework intersects with many
existing works on robust stochastic optimization, where the connection with
median estimation was not known previously.

## References

  * Acharya et al. (2022) Acharya, A., Hashemi, A., Jain, P., Sanghavi, S., Dhillon, I. S., and Topcu, U.  Robust training in high dimensions via block coordinate geometric median descent.  In _Proceedings of The 25th International Conference on Artificial Intelligence and Statistics_ , volume 151 of _Proceedings of Machine Learning Research_ , pp. 11145–11168. PMLR, 28–30 Mar 2022.
  * Alistarh et al. (2018) Alistarh, D., Allen-Zhu, Z., and Li, J.  Byzantine stochastic gradient descent.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper_files/paper/2018/file/a07c2f3b3b907aaf8436a26c6d77f0a2-Paper.pdf.
  * Asi & Duchi (2019) Asi, H. and Duchi, J. C.  Stochastic (approximate) proximal point methods: convergence, optimality, and adaptivity.  _SIAM Journal on Optimization_ , 29(3):2257–2290, 2019.  ISSN 1052-6234.
  * Bauschke & Combettes (2017) Bauschke, H. H. and Combettes, P. L.  _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_.  Springer, 2nd edition edition, 2017.  ISBN 978-1-4419-9467-7.
  * Beck (2017) Beck, A.  _First-order methods in optimization_ , volume 25 of _MOS-SIAM Series on Optimization_.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Optimization Society, Philadelphia, PA, 2017.  ISBN 978-1-611974-98-0.  doi: 10.1137/1.9781611974997.ch1.
  * Bernstein et al. (2018) Bernstein, J., Wang, Y.-X., Azizzadenesheli, K., and Anandkumar, A.  signSGD: Compressed optimisation for non-convex problems.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 560–569. PMLR, 10–15 Jul 2018.  URL https://proceedings.mlr.press/v80/bernstein18a.html.
  * Bickel (1967) Bickel, P. J.  Some contributions to the theory of order statistics.  In _Proc. Fifth Berkeley Sympos. Math. Statist. and Probability (Berkeley, Calif., 1965/66), Vol. I: Statistics_ , pp. 575–591. Univ. California Press, Berkeley, CA, 1967.
  * Bottou (2010) Bottou, L.  Large-scale machine learning with stochastic gradient descent.  In _Proceedings of COMPSTAT’2010_ , pp. 177–186. Physica-Verlag/Springer, Heidelberg, 2010.
  * Bottou et al. (2018) Bottou, L., Curtis, F. E., and Nocedal, J.  Optimization methods for large-scale machine learning.  _SIAM Review_ , 60(2):223–311, 2018.  ISSN 0036-1445.  doi: 10.1137/16M1080173.
  * Burke & Ferris (1993) Burke, J. and Ferris, M.  Weak Sharp Minima in Mathematical Programming.  _SIAM Journal on Control and Optimization_ , 31(5):1340–1359, 1993.  ISSN 0363-0129.  doi: 10.1137/0331063.
  * Cardot et al. (2011) Cardot, H., Cénac, P., and Zitt, P.-A.  Efficient and fast estimation of the geometric median in hilbert spaces with an averaged stochastic gradient algorithm, 2011.
  * Catoni & Giulini (2018) Catoni, O. and Giulini, I.  Dimension-free PAC-bayesian bounds for the estimation of the mean of a random vector.  February 2018.  doi: 10.48550/ARXIV.1802.04308.
  * Chu & Hotelling (1955) Chu, J. T. and Hotelling, H.  The moments of the sample median.  _The Annals of Mathematical Statistics_ , 26(4):593–606, 1955.  ISSN 00034851.  URL http://www.jstor.org/stable/2236373.
  * Cohen et al. (2016) Cohen, M. B., Lee, Y. T., Miller, G., Pachocki, J., and Sidford, A.  Geometric median in nearly linear time.  In _Proceedings of the Forty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’16, pp. 9–21, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450341325.
  * Cornejo et al. (1997) Cornejo, O., Jourani, A., and Zălinescu, C.  Conditioning and Upper-Lipschitz Inverse Subdifferentials in Nonsmooth Optimization Problems.  _Journal of Optimization Theory and Applications_ , 95(1):127–148, 1997.  ISSN 0022-3239, 1573-2878.  doi: 10.1023/A:1022687412779.
  * Cramér (2016) Cramér, H.  _Mathematical Methods of Statistics (PMS-9), Volume 9_.  Princeton University Press, 2016.  ISBN 9781400883868.
  * Data et al. (2018) Data, D., Song, L., and Diggavi, S.  Data encoding for byzantine-resilient distributed gradient descent.  In _2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pp. 863–870, 2018.  doi: 10.1109/ALLERTON.2018.8636017.
  * Davis & Drusvyatskiy (2019) Davis, D. and Drusvyatskiy, D.  Stochastic model-based minimization of weakly convex functions.  _SIAM Journal on Optimization_ , 29(1):207–239, 2019.  ISSN 1052-6234.  doi: 10.1137/18M1178244.
  * Demidovich et al. (2023) Demidovich, Y., Malinovsky, G., Sokolov, I., and Richtárik, P.  A Guide Through the Zoo of Biased SGD, 2023.
  * Donoho & Huber (1983) Donoho, D. and Huber, P. J.  The notion of breakdown point.  In _A Festschrift for Erich L. Lehmann_ , Wadsworth Statist./Probab. Ser., pp. 157–184. Wadsworth, Belmont, CA, 1983.  ISBN 0-534-98044-9.
  * Ferris (1991) Ferris, M. C.  Finite termination of the proximal point algorithm.  _Mathematical Programming_ , 50(1-3):359–366, 1991.  ISSN 0025-5610, 1436-4646.  doi: 10.1007/BF01594944.
  * Fréchet (1948) Fréchet, M. R.  Les éléments aléatoires de nature quelconque dans un espace distancié.  1948\.  URL https://api.semanticscholar.org/CorpusID:124022420.
  * Gorbunov et al. (2020) Gorbunov, E., Danilova, M., and Gasnikov, A.  Stochastic optimization with heavy-tailed noise via accelerated gradient clipping.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems_ , volume 33, pp. 15042–15053. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper_files/paper/2020/file/abd1c782880cc59759f4112fda0b8f98-Paper.pdf.
  * Gurbuzbalaban et al. (2021) Gurbuzbalaban, M., Simsekli, U., and Zhu, L.  The heavy-tail phenomenon in sgd.  In Meila, M. and Zhang, T. (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 3964–3975. PMLR, 18–24 Jul 2021.
  * Hampel et al. (2005) Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J., and Stahel, W. A.  _Robust Statistics: The Approach Based on Influence Functions_.  Wiley, March 2005.  ISBN 9781118186435.  doi: 10.1002/9781118186435.
  * Hodgkinson & Mahoney (2021) Hodgkinson, L. and Mahoney, M.  Multiplicative noise and heavy tails in stochastic optimization.  In _International Conference on Machine Learning_ , pp. 4262–4274. PMLR, 2021.
  * Huber (1981) Huber, P. J.  _Robust statistics_.  Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, Inc., New York, 1981.  ISBN 0-471-41805-6.
  * Karimireddy et al. (2019) Karimireddy, S. P., Rebjock, Q., Stich, S., and Jaggi, M.  Error feedback fixes SignSGD and other gradient compression schemes.  In Chaudhuri, K. and Salakhutdinov, R. (eds.), _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 3252–3261. PMLR, 09–15 Jun 2019.  URL https://proceedings.mlr.press/v97/karimireddy19a.html.
  * Karimireddy et al. (2021) Karimireddy, S. P., He, L., and Jaggi, M.  Learning from history for byzantine robust optimization.  In Meila, M. and Zhang, T. (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 5311–5319. PMLR, 18–24 Jul 2021.  URL https://proceedings.mlr.press/v139/karimireddy21a.html.
  * Khirirat et al. (2023) Khirirat, S., Gorbunov, E., Horváth, S., Islamov, R., Karray, F., and Richtárik, P.  Clip21: Error feedback for gradient clipping, 2023.
  * Koloskova et al. (2023) Koloskova, A., Hendrikx, H., and Stich, S. U.  Revisiting gradient clipping: Stochastic bias and tight convergence guarantees.  In Krause, A., Brunskill, E., Cho, K., Engelhardt, B., Sabato, S., and Scarlett, J. (eds.), _Proceedings of the 40th International Conference on Machine Learning_ , volume 202 of _Proceedings of Machine Learning Research_ , pp. 17343–17363. PMLR, 23–29 Jul 2023.  URL https://proceedings.mlr.press/v202/koloskova23a.html.
  * Kunstner et al. (2023) Kunstner, F., Chen, J., Lavington, J. W., and Schmidt, M.  Noise is not the main factor behind the gap between sgd and adam on transformers, but sign descent might be.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=a65YK0cqH8g.
  * Lemaire (1992) Lemaire, B.  About the Convergence of the Proximal Method.  In Oettli, P. D. W. and Pallaschke, P. D. D. (eds.), _Advances in Optimization_ , number 382 in Lecture Notes in Economics and Mathematical Systems, pp. 39–51. Springer Berlin Heidelberg, 1992.  ISBN 978-3-540-55446-2 978-3-642-51682-5.  doi: 10.1007/978-3-642-51682-5˙4.
  * Lemaire (1998) Lemaire, B.  Well-posedness, conditioning and regularization of minimization, inclusion and fixed-point problems.  _Pliska Studia Mathematica Bulgarica_ , 12(1):71p–84p, 1998.
  * Lewis & Pang (1998) Lewis, A. S. and Pang, J.-S.  Error bounds for convex inequality systems.  _Generalized convexity, generalized monotonicity: recent results_ , pp. 75–110, 1998.
  * Lopuhaa & Rousseeuw (1991) Lopuhaa, H. P. and Rousseeuw, P. J.  Breakdown points of affine equivariant estimators of multivariate location and covariance matrices.  _The Annals of Statistics_ , 19(1):229–248, 1991\.  ISSN 00905364.  URL http://www.jstor.org/stable/2241852.
  * Lugosi & Mendelson (2019) Lugosi, G. and Mendelson, S.  Mean estimation and regression under heavy-tailed distributions: A survey.  _Foundations of Computational Mathematics_ , 19(5):1145–1190, August 2019.  ISSN 1615-3383.  doi: 10.1007/s10208-019-09427-x.
  * Maritz & Jarrett (1978) Maritz, J. and Jarrett, R.  A note on estimating the variance of the sample median.  _Journal of the American Statistical Association_ , 73(361):194–196, 1978.
  * Mhamdi et al. (2018) Mhamdi, E. M. E., Guerraoui, R., and Rouault, S.  The Hidden Vulnerability of Distributed Learning in Byzantium.  In _Proceedings of the 35th International Conference on Machine Learning_ , pp. 3518–3527. PMLR, 2018.
  * Milzarek et al. (2022) Milzarek, A., Schaipp, F., and Ulbrich, M.  A semismooth Newton stochastic proximal point algorithm with variance reduction.  April 2022.
  * Minsker (2015) Minsker, S.  Geometric median and robust estimation in banach spaces.  _Bernoulli_ , 21(4), November 2015.  ISSN 1350-7265.  doi: 10.3150/14-bej645.
  * Nemirovsky & Yudin (1983) Nemirovsky, A. S. and Yudin, D. B. a.  _Problem complexity and method efficiency in optimization_.  A Wiley-Interscience Publication. John Wiley & Sons, Inc., New York, 1983\.  ISBN 0-471-10345-4.  Translated from the Russian and with a preface by E. R. Dawson, Wiley-Interscience Series in Discrete Mathematics.
  * Nolan (2013) Nolan, J. P.  Multivariate elliptically contoured stable distributions: theory and estimation.  _Computational Statistics_ , 28(5):2067–2089, January 2013.  ISSN 1613-9658.  doi: 10.1007/s00180-013-0396-7.
  * Nolan (2020) Nolan, J. P.  _Univariate Stable Distributions_.  Springer International Publishing, 2020.  doi: 10.1007/978-3-030-52915-4.
  * Pavasovic et al. (2023) Pavasovic, K. L., Durmus, A., and Simsekli, U.  Approximate heavy tails in offline (multi-pass) stochastic gradient descent.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.  URL https://openreview.net/forum?id=UkPeUXML7s.
  * Polyak (1964) Polyak, B. T.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, 1964.  ISSN 0041-5553.
  * Puchkin et al. (2023) Puchkin, N., Gorbunov, E., Kutuzov, N., and Gasnikov, A.  Breaking the heavy-tailed noise barrier in stochastic optimization problems.  November 2023.  doi: 10.48550/ARXIV.2311.04161.
  * Richtárik et al. (2021) Richtárik, P., Sokolov, I., and Fatkhullin, I.  EF21: A new, simpler, theoretically better, and practically faster error feedback.  In Beygelzimer, A., Dauphin, Y., Liang, P., and Vaughan, J. W. (eds.), _Advances in Neural Information Processing Systems_ , 2021.
  * Riedmiller & Braun (1993) Riedmiller, M. and Braun, H.  A direct adaptive method for faster backpropagation learning: the rprop algorithm.  In _IEEE International Conference on Neural Networks_ , pp. 586–591 vol.1, 1993.  doi: 10.1109/ICNN.1993.298623.
  * Sadiev et al. (2023) Sadiev, A., Danilova, M., Gorbunov, E., Horváth, S., Gidel, G., Dvurechensky, P., Gasnikov, A., and Richtárik, P.  High-probability bounds for stochastic optimization and variational inequalities: the case of unbounded variance.  _arXiv preprint arXiv:2302.00999_ , 2023.
  * Safaryan et al. (2021) Safaryan, M., Shulgin, E., and Richtárik, P.  Uncertainty principle for communication compression in distributed and federated learning and the search for an optimal compressor.  _Information and Inference: A Journal of the IMA_ , 11(2):557–580, 04 2021.
  * Seide et al. (2014) Seide, F., Fu, H., Droppo, J., Li, G., and Yu, D.  1-bit stochastic gradient descent and application to data-parallel distributed training of speech dnns.  In _Interspeech 2014_ , September 2014.
  * Simsekli et al. (2020) Simsekli, U., Sener, O., Deligiannidis, G., and Erdogdu, M. A.  Hausdorff dimension, heavy tails, and generalization in neural networks.  _Advances in Neural Information Processing Systems_ , 33:5138–5151, 2020.
  * Sutskever et al. (2013) Sutskever, I., Martens, J., Dahl, G., and Hinton, G.  On the importance of initialization and momentum in deep learning.  In Dasgupta, S. and McAllester, D. (eds.), _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _Proceedings of Machine Learning Research_ , pp. 1139–1147, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * Toulis et al. (2020) Toulis, P., Horel, T., and Airoldi, E. M.  The proximal robbins–monro method.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 83(1):188–212, December 2020.  ISSN 1467-9868.  doi: 10.1111/rssb.12405.
  * Vardi & Zhang (2000) Vardi, Y. and Zhang, C.-H.  The multivariate ${L}_{1}$-median and associated data depth.  _Proceedings of the National Academy of Sciences_ , 97(4):1423–1426, 2000.  doi: 10.1073/pnas.97.4.1423.
  * Wang et al. (2021) Wang, H., Gurbuzbalaban, M., Zhu, L., Simsekli, U., and Erdogdu, M. A.  Convergence rates of stochastic gradient descent under infinite noise variance.  _Advances in Neural Information Processing Systems_ , 34:18866–18877, 2021.
  * Weiszfeld & Plastria (2009) Weiszfeld, E. and Plastria, F.  On the point for which the sum of the distances to _n_ given points is minimum.  _Ann. Oper. Res._ , 167(1):7–41, 2009.  doi: 10.1007/s10479-008-0352-z.  URL https://doi.org/10.1007/s10479-008-0352-z.
  * Yin et al. (2018) Yin, D., Chen, Y., Kannan, R., and Bartlett, P.  Byzantine-robust distributed learning: Towards optimal statistical rates.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 5650–5659. PMLR, 10–15 Jul 2018.
  * Zhang et al. (2020a) Zhang, J., He, T., Sra, S., and Jadbabaie, A.  Why gradient clipping accelerates training: A theoretical justification for adaptivity.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020a.  URL https://openreview.net/forum?id=BJgnXpVYwS.
  * Zhang et al. (2020b) Zhang, J., Karimireddy, S. P., Veit, A., Kim, S., Reddi, S., Kumar, S., and Sra, S.  Why are adaptive methods good for attention models?  In _Proceedings of the 34th International Conference on Neural Information Processing Systems_ , NeurIPS 2020, Red Hook, NY, USA, 2020b. Curran Associates Inc.  ISBN 9781713829546.

SGD with Clipping is Secretly Estimating the Median Gradient
APPENDIX

The Appendix is organized as follows.

  * •

In Appendix A we provide missing proofs for the sample median method gradient
descent (SMGD) presented in Section 3.

  * •

In Appendix B, we provide additional formulas and derivations for the problem
of approximating the median, studied in Section 4.

  * •

In Appendix C, we provide additional plots and supplementary information on
numerical experiments.

##  Appendix A Missing Proofs for Section 3

###  A.1 On the Assumptions on the Sample Median

We first state a simple fact about Lipschitz functions:

######  Lemma A.1.

If Assumption 3.1 holds, and if $\ell$ is Lipschitz continuous, then without
loss of generality we can assume that $\delta_{2}=\sigma_{2}=0$.

###### Proof.

If $\ell$ is $G$-Lipschitz, then for every subgradient
$\bm{u}\in\partial\ell(\bm{w})$ we have $\|\bm{u}\|\leq G$, so one can replace
$\delta_{1}^{2}$ with $\delta_{1}^{2}+\delta_{2}G^{2}$ and $\sigma_{1}^{2}$
with $\sigma_{1}^{2}+\sigma_{2}G^{2}$. ∎

In the following proof of Proposition 3.3 we will need an equivalent but more
compact formulation of Assumption 3.1.

######  Lemma A.2.

Let $\bm{m}$ be a random variable in $\mathbb{R}^{d}$, and consider the
following property:

$\mathbb{E}\left[\|\bm{m}-\nabla\ell(\bm{w})\|^{2}\right]\leq\nu_{1}^{2}+\nu_{2}\|\nabla\ell(\bm{w})\|^{2}.$
(20)

If (6) holds then (20) holds with $\nu_{1}^{2}=\delta_{1}^{2}+\sigma_{1}^{2}$
and $\nu_{2}=\delta_{2}+\sigma_{2}$. If (20) holds, then (6) holds with
$\delta_{1}=\sigma_{1}=\nu_{1}$ and $\delta_{2}=\sigma_{2}=\nu_{2}$.

###### Proof.

Expanding the squares, we get

$\mathbb{E}\left[\|\bm{m}-\nabla\ell(\bm{w})\|^{2}\right]=\mathbb{E}\left[\|\bm{m}-\mathbb{E}\left[\bm{m}\right]\|^{2}\right]+\|\mathbb{E}\left[\bm{m}\right]-\nabla\ell(\bm{w})\|^{2},$

where we used that
$\mathbb{E}\left[\langle\bm{m}-\mathbb{E}\left[\bm{m}\right],\mathbb{E}\left[\bm{m}\right]-\nabla\ell(\bm{w})\rangle\right]=0$.
Now, the first implication is trivial, and the second follows from the
positivity of the involved terms. ∎

We now verify that Assumption 3.1 holds true for the $\ell_{1}$-sample median,
provided that the noise on the subgradient is linearly depending on an
$\alpha$-stable distribution, and that the linear coefficients grow with the
norm of the subgradient.

######  Lemma A.3.

Let $n\in\mathbb{N}$ be odd and $n\geq 3$. Suppose that the subgradient
oracles $\bm{g}^{(i)}=\nabla\ell(\bm{w})+\bm{\Sigma}(\bm{w})\bm{\zeta}^{(i)}$
verify Assumption 3.5, for all $i=1,\dots,n$ and set

$\bm{m}=\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{g}^{(i)}(\bm{w})\right).$

Then, Assumption 3.1 holds with $\delta_{1}=\delta_{2}=0$, and
$\sigma_{1}^{2}=C_{n}c_{1}$, $\sigma_{2}=C_{n}c_{2}$, for some constant
$C_{n}=\mathcal{O}(1/n)$.

###### Proof.

Let us denote the noise vector by
$\bm{\xi}^{(i)}=\bm{\Sigma}(\bm{w})\bm{\zeta}^{(i)}$. By the stability
property of $\alpha$-stable distributions (Nolan, 2020, Proposition 1.3), we
have that

$\displaystyle{\xi}_{k}^{(i)}:=\bm{\xi}_{k}^{(i)}(\bm{w})=\sum_{j=1}^{d}[\bm{\Sigma}(\bm{w})]_{kj}\bm{\zeta}^{(i)}_{j}=^{\mathrm{d}}\|\bm{\Sigma}(\bm{w})_{k}\|_{\alpha}z,$
(21)

where $=^{\mathrm{d}}$ denotes equality in distribution, $z$ is a standard
symmetric $\alpha$-stable random variable in $\mathbb{R}$, and
$\bm{\Sigma}(\bm{w})_{k}$ denotes the $k$-th row of $\bm{\Sigma}(\bm{w})$.
Observe that from our Assumption 3.5 we can write that

$\sum_{i=1}^{d}\|\bm{\Sigma}(\bm{w})_{i}\|_{\alpha}^{2}\leq\rho_{\alpha}^{2}\sum_{i=1}^{d}\|\bm{\Sigma}(\bm{w})_{i}\|_{2}^{2}=\rho_{\alpha}^{2}\|\bm{\Sigma}(\bm{w})\|_{F}^{2}\leq
c_{1}^{\prime}+c_{2}^{\prime}\|\nabla\ell(\bm{w})\|^{2},\quad\mbox{for all
}\bm{w}\in\mathbb{R}^{d}.$ (22)

by setting $c_{1}^{\prime}=\rho_{\alpha}^{2}c_{1}$ and
$c_{2}=\rho_{\alpha}^{2}c_{2}$, and where $\rho_{\alpha}>0$ is the best
constant such that $\|\cdot\|_{\alpha}\leq\rho_{\alpha}\|\cdot\|_{2}$, which
is $\rho_{\alpha}=d^{\tfrac{2-\alpha}{2\alpha}}$ is $\alpha\leq 2$, and
$\rho_{\alpha}=1$ otherwise.

We start by computing the bias term that is required by Assumption 3.1.

$\displaystyle\|\mathbb{E}\left[\bm{m}\right]-\nabla\ell(\bm{w})\|^{2}=$
$\displaystyle\|\mathbb{E}\left[\texttt{median}_{\ell_{1},i\in[n]}\left(\nabla\ell(\bm{w})+\bm{\xi}^{(i)}(\bm{w})\right)\right]-\nabla\ell(\bm{w})\|^{2}$
$\displaystyle=$
$\displaystyle\|\mathbb{E}\left[\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{\xi}^{(i)}(\bm{w})\right)\right]\|^{2}$

where we used the fact that

$\texttt{median}_{\ell_{1},i\in[n]}\left(\nabla\ell(\bm{w})+\bm{\xi}^{(i)}(\bm{w})\right)=\nabla\ell(\bm{w})+\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{\xi}^{(i)}(\bm{w})\right).$

On the other hand, the $\ell_{1}$-median boils down to computing component-
wise medians, hence we have that

$\displaystyle\mathbb{E}\left[\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{\xi}^{(i)}(\bm{w})\right)\right]=\begin{bmatrix}\mathbb{E}\left[\texttt{median}_{i\in[n]}\left(\xi_{1}^{(i)}\right)\right]\\\
\vdots\\\
\mathbb{E}\left[\texttt{median}_{i\in[n]}\left(\xi_{d}^{(i)}\right)\right]\end{bmatrix},$

where
$\texttt{median}_{i\in[n]}(\xi_{k}^{(i)})=\underset{m\in\mathbb{R}}{\operatorname{argmin}}\;\frac{1}{n}\sum_{i=1}^{n}|m-\xi_{k}^{(i)}|$
is the one dimensional sample median estimator.

Now, let us compute
$\mathbb{E}\left[\texttt{median}_{i\in[n]}\left(\xi^{(i)}_{k}\right)\right]$
for some $k\in\\{1,\dots,d\\}$. Denoting $r=(n-1)/2\in\mathbb{N}$ and by using
the probability density function of the sample median (cf. (Maritz & Jarrett,
1978)), we have that:

$\displaystyle\mathbb{E}\left[\texttt{median}_{i\in[n]}\left(\xi^{(i)}_{k}\right)\right]=$
$\displaystyle\frac{n!}{r!r!}\underbrace{\int_{\mathbb{R}}x\left(F_{\bm{w},k}(x)(1-F_{\bm{w},k}(x))\right)^{r}p_{\bm{w},k}(x)\mathrm{d}x}_{=:\mathrm{A}},$

where $F_{\bm{w},k}$ and $p_{\bm{w},k}(x)$ are respectively the cumulative
distribution function (cdf) and the probability density function (pdf) of
$\xi_{k}^{(i)}\sim\mathcal{P}_{\alpha,\|\bm{\Sigma}_{k}(\bm{w})\|_{\alpha}}$
that is a symmetric $\alpha$-stable random variable with scale
$\|\bm{\Sigma}_{k}(\bm{w})\|_{\alpha}$. As $p_{\bm{w},k}$ is continuous and
symmetric around zero, we have $p_{\bm{w},k}(x)=p_{\bm{w},k}(-x)$ and
$(1-F_{\bm{w},k}(x))=F_{\bm{w},k}(-x)$ for all $x\in\mathbb{R}$. This gives
us:

$\displaystyle\mathrm{A}=$
$\displaystyle\int_{0}^{\infty}x\left(F_{\bm{w},k}(x)F_{\bm{w},k}(-x)\right)^{r}p_{\bm{w},k}(x)\mathrm{d}x+\int_{-\infty}^{0}x\left(F_{\bm{w},k}(x)F_{\bm{w},k}(-x)\right)^{r}p_{\bm{w},k}(x)\mathrm{d}x$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}x\left(F_{\bm{w},k}(x)F_{\bm{w},k}(-x)\right)^{r}p_{\bm{w},k}(x)\mathrm{d}x-\int_{0}^{\infty}x\left(F_{\bm{w},k}(x)F_{\bm{w},k}(-x)\right)^{r}p_{\bm{w},k}(x)\mathrm{d}x$
$\displaystyle=$ $\displaystyle 0.$

Hence, we obtain that

$\mathbb{E}\left[\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{\xi}^{(i)}(\bm{w})\right)\right]=0.$

This shows that the first line of Assumption 3.1 holds with
$\delta_{1}=\delta_{2}=0$.

We now proceed to the estimation of the variance. By using similar arguments,
we have that:

$\displaystyle\mathbb{E}\left[\|\bm{m}-\mathbb{E}\left[\bm{m}\right]\|^{2}\right]=\;$
$\displaystyle\mathbb{E}\left[\|\texttt{median}_{\ell_{1},i\in[n]}\left(\bm{\xi}^{(i)}(\bm{w})\right)\|^{2}\right]$
$\displaystyle=\;$
$\displaystyle\sum_{k=1}^{d}\mathbb{E}\left[\left(\texttt{median}_{i\in[n]}\left(\xi_{k}^{(i)}\right)\right)^{2}\right]$
$\displaystyle=\;$
$\displaystyle\sum_{k=1}^{d}\frac{n!}{r!r!}\mathbb{E}\left[\Bigl{(}\xi_{k}^{(i)}\Bigr{)}^{2}\Biggl{(}F_{\bm{w},k}\Bigl{(}\xi_{k}^{(i)}\Bigr{)}\Bigl{(}1-F_{\bm{w},k}\Bigl{(}\xi_{k}^{(i)}\Bigr{)}\Bigr{)}\Biggr{)}^{r}\right],$
(23)

where (23) follows from (Maritz & Jarrett, 1978, Equation 2.1). By using (21)
with (22) and denoting the cdf of the standard symmetric $\alpha$-stable
distribution by $F_{z}$, we further compute that:

$\displaystyle\mathbb{E}\left[\|\bm{m}-\mathbb{E}\left[\bm{m}\right]\|^{2}\right]=\;$
$\displaystyle\frac{n!}{r!r!}\mathbb{E}\left[z^{2}\Biggl{(}F_{z}(z)(1-F_{z}(z)\Biggr{)}^{r}\right]\sum_{k=1}^{d}\|\bm{\Sigma}_{k}(\bm{w})\|_{\alpha}^{2}$
$\displaystyle\leq\;$
$\displaystyle\frac{n!}{r!r!}\mathbb{E}\left[z^{2}\Biggl{(}F_{z}(z)(1-F_{z}(z)\Biggr{)}^{r}\right]\left(c_{1}^{\prime}+c_{2}^{\prime}\|\nabla\ell(\bm{w})\|^{2}\right)$
$\displaystyle=:\;$ $\displaystyle
C_{n}\left(c_{1}^{\prime}+c_{2}^{\prime}\|\nabla\ell(\bm{w})\|^{2}\right).$

This shows that the second line of Assumption 3.1 holds with
$\sigma_{1}^{2}=C_{n}\rho_{\alpha}^{2}c_{1}$ and
$\sigma_{2}=C_{n}\rho_{\alpha}^{2}c_{2}$.

Finally we observe that, by again (Maritz & Jarrett, 1978, Equation 2.1),
$C_{n}$ is the variance of the sample median of a set of i.i.d. standard
symmetric $\alpha$-stable variables $\\{z_{1},\dots,z_{n}\\}$, where
$z_{i}=^{\mathrm{d}}z$. As $n\to\infty$, it is well-known that the sample
median is asymptotically normal with variance $1/(4nf^{2}(0))$, where $f$
denotes the pdf of $z$ (see (Maritz & Jarrett, 1978)). Hence we conclude that
$C_{n}=\mathcal{O}(1/n)$ as $n\to\infty$. This concludes the proof. ∎

###  A.2 Rates for SMGD for Smooth Nonconvex Functions: Proof of Proposition
3.3

The following proposition is an application of biased SGD results to the SMGD
algorithm with errors. It is possible to show that Assumptions 3.1 and 3.2
both imply Assumption 9 from (Demidovich et al., 2023). So from a qualitative
point of view the result below is the same that what we would obtain by using
Theorem 3 from (Demidovich et al., 2023). But doing so would produce a bound
with worse constants, so we prefer to provide a proof exploiting directly our
assumptions.

See 3.3

###### Proof.

In this proof we note $\mathbb{E}_{t}$ for the expectation taken conditionally
to the filtration generated by $\bm{w}_{0},\dots\bm{w}_{t}$. In particular, we
will have from our Assumption 3.2 on the errors together with the tower rule
that
$\mathbb{E}_{t}\|\bm{e}_{t}\|^{2}=\mathbb{E}_{t}\mathbb{E}\left[\|\bm{e}_{t}\|^{2}\
|\ \bm{g}_{t}^{(1)},\dots,\bm{g}_{t}^{(n)}\right]\leq\varepsilon^{2}$.

Smoothness of $\ell$ implies

$\ell(\bm{w}_{t+1})-\ell(\bm{w}_{t})-\langle\nabla\ell(\bm{w}_{t}),\bm{w}_{t+1}-\bm{w}_{t}\rangle\leq\frac{L}{2}\|\bm{w}_{t+1}-\bm{w}_{t}\|^{2}.$

If we note $\bm{d}_{t}:=\bm{m}_{t}+\bm{e}_{t}$, thus
$\bm{w}_{t+1}-\bm{w}_{t}=-\eta\bm{d}_{t}$. This implies

$\ell(\bm{w}_{t+1})-\ell(\bm{w}_{t})+\eta\langle\nabla\ell(\bm{w}_{t}),\bm{d}_{t}\rangle\leq\frac{L\eta^{2}}{2}\|\bm{d}_{t}\|^{2}.$

Let us bound the two terms depending on $\bm{d}_{t}$. First, take expectation
conditioned to $\mathcal{F}_{t}$ (which we will denote $\mathbb{E}_{t}$) to
write

$\displaystyle\mathbb{E}_{t}\left[\langle\nabla\ell(\bm{w}_{t}),\bm{d}_{t}\rangle\right]$
$\displaystyle=$
$\displaystyle\langle\nabla\ell(\bm{w}_{t}),\mathbb{E}_{t}\left[\bm{d}_{t}\right]\rangle=\frac{1}{2}\|\nabla\ell(\bm{w}_{t})\|^{2}+\frac{1}{2}\|\mathbb{E}_{t}\left[\bm{d}_{t}\right]\|^{2}-\frac{1}{2}\|\mathbb{E}_{t}\left[\bm{d}_{t}\right]-\nabla\ell(\bm{w}_{t})\|^{2}$
$\displaystyle\geq$
$\displaystyle\frac{1}{2}\|\nabla\ell(\bm{w}_{t})\|^{2}-\frac{1}{2}\|\mathbb{E}_{t}\left[\bm{d}_{t}\right]-\nabla\ell(\bm{w}_{t})\|^{2}.$

Using (6) and
$\|\mathbb{E}_{t}\left[\bm{e}_{t}\right]\|^{2}\leq\mathbb{E}_{t}\left[\|\bm{e}_{t}\|^{2}\right]\leq\varepsilon^{2}$
due to Jensen’s inequality, we can bound

$\|\mathbb{E}_{t}\left[\bm{d}_{t}\right]-\nabla\ell(\bm{w}_{t})\|^{2}\leq
2\|\mathbb{E}_{t}\left[\bm{m}_{t}\right]-\nabla\ell(\bm{w}_{t})\|^{2}+2\|\mathbb{E}_{t}\left[\bm{e}_{t}\right]\|^{2}\leq
2\delta_{1}^{2}+2\delta_{2}\|\nabla\ell(\bm{w}_{t})\|^{2}+2\varepsilon^{2}.$

Second, we use Young’s inequality and our assumptions to bound
$\|\bm{d}_{t}\|^{2}$:

$\displaystyle\mathbb{E}_{t}\left[\|\bm{d}_{t}\|^{2}\right]$
$\displaystyle\leq$ $\displaystyle
2\mathbb{E}_{t}\left[\|\bm{e}_{t}+\bm{m}_{t}-\nabla\ell(\bm{w}_{t})\|^{2}\right]+2\|\nabla\ell(\bm{w}_{t})\|^{2}$
$\displaystyle\leq$ $\displaystyle
4\mathbb{E}_{t}\left[\|\bm{e}_{t}\|^{2}\right]+4\mathbb{E}_{t}\left[\|\bm{m}_{t}-\nabla\ell(\bm{w}_{t})\|^{2}\right]+2\|\nabla\ell(\bm{w}_{t})\|^{2}$
$\displaystyle\leq$ $\displaystyle
4\varepsilon^{2}+4(\delta_{1}^{2}+\sigma_{1}^{2})+(2+4\delta_{2}+4\sigma_{2})\|\nabla\ell(\bm{w}_{t})\|^{2},$

where in the last inequality we used that (6) implies (20) (cf. Lemma A.2).
Combine all the above inequalities to get

$\displaystyle\mathbb{E}_{t}\left[\ell(\bm{w}_{t+1})\right]-\ell(\bm{w}_{t})+\tfrac{\eta}{2}\|\nabla\ell(\bm{w}_{t})\|^{2}-\tfrac{\eta}{2}(2\delta_{1}^{2}+2\delta_{2}\|\nabla\ell(\bm{w}_{t})\|^{2}+\varepsilon^{2})$
$\displaystyle\hskip
8.61108pt\leq\frac{L\eta^{2}}{2}\big{(}4(\varepsilon^{2}+\delta_{1}^{2}+\sigma_{1}^{2})+(2+4\delta_{2}+4\sigma_{2})\|\nabla\ell(\bm{w}_{t})\|^{2}\big{)}.$

Rearranging yields

$\displaystyle\mathbb{E}_{t}\left[\ell(\bm{w}_{t+1})\right]-\ell(\bm{w}_{t})$
$\displaystyle\hskip
12.91663pt\leq-\eta\big{[}\tfrac{1}{2}-\delta_{2}-L\eta(1+2\sigma_{2}+2\delta_{2})\big{]}\|\nabla\ell(\bm{w}_{t})\|^{2}+\eta(\delta_{1}^{2}+\varepsilon^{2})+2L\eta^{2}(\varepsilon^{2}+\delta_{1}^{2}+\sigma_{1}^{2}).$

If $\delta_{2}\leq\tfrac{1}{8}$ and using
$\eta\leq\frac{1}{8L(1+2\sigma_{2}+2\delta_{2})}$, after taking full
expectation, we obtain

$\displaystyle\mathbb{E}\left[\ell(\bm{w}_{t+1})\right]-\mathbb{E}\left[\ell(\bm{w}_{t})\right]$
$\displaystyle\leq$
$\displaystyle-\tfrac{\eta}{4}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]+\eta(\delta_{1}^{2}+\varepsilon^{2})+2L\eta^{2}(\varepsilon^{2}+\delta_{1}^{2}+\sigma_{1}^{2}).$

Sum over $t=0,\dots,T-1$, multiply by $\frac{4}{\eta T}$ and use
$\min\limits_{t=0,\dots,T-1}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]\leq\tfrac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]$
to finally obtain

$\min\limits_{t=0,\dots,T-1}\mathbb{E}\left[\|\nabla\ell(\bm{w}_{t})\|^{2}\right]\leq\frac{4(\ell(\bm{w}_{0})-\inf\ell)}{\eta
T}+8L\eta(\varepsilon^{2}+\delta_{1}^{2}+\sigma_{1}^{2})+4(\delta_{1}^{2}+\varepsilon^{2}).$

∎

###  A.3 Rates for Biased SGD for Convex Lipschitz Functions

In this section we prove another result about the complexity of biased SGD,
that is SGD with biased estimators of the subgradient. As far as we know,
biased SGD was only studied in the case of strongly convex functions or smooth
nonconvex functions (Demidovich et al., 2023), or in the particular case of
random projections. In that regard, the next result is a new contribution to
the analysis of biased SGD.

Before stating our result, we recall a few standard definitions given that we
are working with a nonsmooth convex function. Given a convex function
$\ell:\mathbb{R}^{d}\to\mathbb{R}$, we say that $\bm{u}$ is a subgradient of
$\ell$ at $\bm{w}$ if

$(\forall\bm{w}^{\prime}\in\mathbb{R}^{d})\ell(\bm{w}^{\prime})-\ell(\bm{w})-\langle\bm{u},\bm{w}^{\prime}-\bm{w}\rangle\geq
0.$

We note $\partial\ell(\bm{w})$ the set of all subgradients of $\ell$ at
$\bm{w}$, that is called the subdifferential of $\ell$ at $\bm{w}$.

######  Theorem A.4.

Let $\ell:\mathbb{R}^{d}\to\mathbb{R}$ be a convex $G$-Lipschitz function such
that ${\rm{argmin}}~{}\ell\neq\emptyset$. Let $\eta>0$, and note
$D:=\|\bm{w}_{0}-\bm{w}_{*}\|$ for some $\bm{w}_{*}\in{\rm{argmin}}~{}\ell$.
Consider the iterates generated by a biased Stochastic Subgradient Descent
method

$\bm{w}_{t+1}=\bm{w}_{t}-\eta\bm{g}_{t},$

and assume that the estimator $\bm{g}_{t}$ has a uniformly finite bias and
variance (we denote by $\mathcal{F}_{t}$ the filtration induced by
$\bm{w}_{0},\dots,\bm{w}_{t}$):

$\|\mathbb{E}\left[\bm{g}_{t}|\mathcal{F}_{t}\right]-\bm{u}_{t}\|\leq\delta,\quad\mathbb{E}\left[\|\bm{g}_{t}-\mathbb{E}\left[\bm{g}_{t}|\mathcal{F}_{t}\right]\|^{2}|\mathcal{F}_{t}\right]\leq\sigma^{2},$

where $\bm{u}_{t}$ is any subgradient of $\ell$ at $\bm{w}_{t}$. Then, for
every $T\geq 1$

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{2D^{2}}{\eta
T}+\eta\frac{\sigma^{2}+(\delta+G)^{2}}{2}+\eta\frac{\delta^{2}T}{4}$

where $\bar{\bm{w}}_{T}$ is an average of the first $T$ iterates defined by
$\bar{\bm{w}}_{T}=\frac{1}{\sum_{t=0}^{T-1}\theta^{t}}\sum_{t=0}^{T-1}\theta^{t}\bm{w}_{t}$,
$\theta=\frac{T}{T+2}$. In particular, if we take a constant step size equal
to

$\eta=\frac{2\sqrt{2}D}{\sqrt{2(\sigma^{2}+(\delta+G)^{2})T+\delta^{2}T^{2}}}\quad\text{(resp.
$\eta=\frac{1}{\delta T+\sqrt{T}}$ )}$

then we can guarantee that

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq
2D\sqrt{\frac{\sigma^{2}+(\delta+G)^{2}}{T}+\frac{\delta^{2}}{2}}\quad\text{(resp.
$\frac{2D^{2}+\sigma^{2}+\delta^{2}+G^{2}}{\sqrt{T}}+(2D^{2}+1)\delta$ ).}$

###### Proof.

Let $w_{*}\in{\rm{argmin}}~{}\ell$ and $T\geq 1$ be fixed. Develop the squares
and use the definition of the algorithm to write

$\displaystyle\|\bm{w}_{t+1}-\bm{w}_{*}\|^{2}$ $\displaystyle=$
$\displaystyle\|\bm{w}_{t}-\bm{w}_{*}\|^{2}+2\langle\bm{w}_{t+1}-\bm{w}_{t},\bm{w}_{t}-\bm{w}_{*}\rangle+\|\bm{w}_{t+1}-\bm{w}_{t}\|^{2}$
$\displaystyle=$
$\displaystyle\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta\langle\bm{g}_{t},\bm{w}_{t}-\bm{w}_{*}\rangle+\eta^{2}\|\bm{g}_{t}\|^{2}.$

Let $\bar{\bm{g}}_{t}$ be the projection of
$\mathbb{E}\left[\bm{g}_{t}|\mathcal{F}_{t}\right]$ onto
$\partial\ell(\bm{w}_{t})$, the subdifferential of $\ell$ at $\bm{w}_{t}$.
This means that from our assumption on the bias of $\bm{g}_{t}$ we will have

$\|\mathbb{E}\left[\bm{g}_{t}|\mathcal{F}_{t}\right]-\bar{\bm{g}}_{t}\|\leq\|\mathbb{E}\left[\bm{g}_{t}|\mathcal{F}_{t}\right]-\bm{u}_{t}\|\leq\delta.$

Then, after taking conditional expectation (denoted by
$\mathbb{E}_{t}\left[\cdot\right]$ instead of
$\mathbb{E}\left[\cdot|\mathcal{F}_{t}\right]$) we can write

$\displaystyle\mathbb{E}_{t}\left[\|\bm{w}_{t+1}-\bm{w}_{*}\|^{2}\right]$
$\displaystyle=$
$\displaystyle\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta\langle\mathbb{E}_{t}\left[\bm{g}_{t}\right],\bm{w}_{t}-\bm{w}_{*}\rangle+\eta^{2}\mathbb{E}_{t}\left[\|\bm{g}_{t}\|^{2}\right]$
$\displaystyle=$
$\displaystyle\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta\langle\bar{\bm{g}}_{t},\bm{w}_{t}-\bm{w}_{*}\rangle+2\eta\langle\bar{\bm{g}}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right],\bm{w}_{t}-\bm{w}_{*}\rangle+\eta^{2}\mathbb{E}_{t}\left[\|\bm{g}_{t}\|^{2}\right]$
$\displaystyle\leq$
$\displaystyle\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta(\ell(\bm{w}_{t})-\inf\ell)+2\eta\langle\bar{\bm{g}}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right],\bm{w}_{t}-\bm{w}_{*}\rangle+\eta^{2}\mathbb{E}_{t}\left[\|\bm{g}_{t}\|^{2}\right],$

where in the last inequality we used the convexity of $\ell$ together with the
fact that $\bar{\bm{g}}_{t}\in\partial\ell(\bm{w}_{t})$. The last term can be
bounded by using our assumptions on the estimator $\bm{g}_{t}$, and the fact
that $\ell$ has bounded subgradients:

$\mathbb{E}_{t}\left[\|\bm{g}_{t}\|^{2}\right]=\mathbb{E}_{t}\left[\|\bm{g}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right]\|^{2}\right]+\|\mathbb{E}_{t}\left[\bm{g}_{t}\right]\|^{2}\leq\sigma^{2}+(\|\mathbb{E}_{t}\left[\bm{g}_{t}\right]-\bar{\bm{g}}_{t}\|+\|\bar{\bm{g}}_{t}\|)^{2}\leq\sigma^{2}+(\delta+G)^{2}.$

Injecting this in the previous inequality we obtain

$\mathbb{E}_{t}\left[\|\bm{w}_{t+1}-\bm{w}_{*}\|^{2}\right]\leq\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta(\ell(\bm{w}_{t})-\inf\ell)+2\eta\langle\bar{\bm{g}}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right],\bm{w}_{t}-\bm{w}_{*}\rangle+\eta^{2}\sigma^{2}+\eta^{2}(\delta+G)^{2}.$

Now we introduce a parameter $\eta>0$, and we use Young’s inequality together
with our assumption on the bias of $\bm{g}_{t}$ to write

$2\langle\bar{\bm{g}}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right],\bm{w}_{t}-\bm{w}_{*}\rangle\leq\frac{1}{\eta}\|\bar{\bm{g}}_{t}-\mathbb{E}_{t}\left[\bm{g}_{t}\right]\|^{2}+\eta\|\bm{w}_{t}-\bm{w}_{*}\|^{2}\leq\eta^{-1}\delta^{2}+\eta\|\bm{w}_{t}-\bm{w}_{*}\|^{2}.$

Inject this bound in our main inequality

$\mathbb{E}_{t}\left[\|\bm{w}_{t+1}-\bm{w}_{*}\|^{2}\right]\leq(1+\eta\eta)\|\bm{w}_{t}-\bm{w}_{*}\|^{2}-2\eta(\ell(\bm{w}_{t})-\inf\ell)+\eta^{2}\sigma^{2}+\eta^{2}(\delta+G)^{2}+\eta^{-1}\delta^{2},$

and after reorganizing the terms, dividing by $2\eta$ and taking expectation
we finally obtain

$\mathbb{E}\left[\ell(\bm{w}_{t})-\inf\ell\right]\leq\frac{1+\eta\eta}{2\eta}\mathbb{E}\left[\|\bm{w}_{t}-\bm{w}_{*}\|^{2}\right]-\frac{1}{2\eta}\mathbb{E}\left[\|\bm{w}_{t+1}-\bm{w}_{*}\|^{2}\right]+\eta\frac{\sigma^{2}+(\delta+G)^{2}}{2}+\eta^{-1}\frac{\delta^{2}}{2}.$

We are now ready to use a weighted telescopic sum argument. First, define
$\theta:=\frac{1}{1+\eta\eta}$ and $\rho_{t}:=\theta^{t}$. Second, multiply
our inequality by $\rho_{t+1}$ and then sum for $t=0,\dots,T-1$. Observe that
the term $(1+\eta\eta)\rho_{t}$ is equal to $\rho_{t}$ due to our definition,
which means that we have a telescoping sum where the terms
$\rho_{t}\mathbb{E}\left[\|\bm{w}_{t}-\bm{w}_{*}\|^{2}\right]$ will cancel
each other. Third, divide by $\sum_{t=0}^{T-1}\rho_{t+1}$ so to obtain

$\frac{1}{\sum_{t=0}^{T-1}\rho_{t+1}}\sum_{t=0}^{T-1}\rho_{t+1}\mathbb{E}\left[\ell(\bm{w}_{t})-\inf\ell\right]\leq\frac{\rho_{0}}{2\eta\sum_{t=0}^{T-1}\rho_{t+1}}\mathbb{E}\left[\|\bm{w}_{0}-\bm{w}_{*}\|^{2}\right]+\eta\frac{\sigma^{2}+(\delta+G)^{2}}{2}+\eta^{-1}\frac{\delta^{2}}{2}.$

Now define
$\bar{\bm{w}}_{T}:=\frac{1}{\sum_{t=0}^{T-1}\rho_{t+1}}\sum_{t=0}^{T-1}\rho_{t+1}\bm{w}_{t}=\frac{1}{\sum_{t=0}^{T-1}\theta^{t}}\sum_{t=0}^{T-1}\theta^{t}\bm{w}_{t}$,
and use Jensen’s inequality to get

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{\|\bm{w}_{0}-\bm{w}_{*}\|^{2}}{2\eta\sum_{t=0}^{T-1}\rho_{t+1}}+\eta\frac{\sigma^{2}+(\delta+G)^{2}}{2}+\eta^{-1}\frac{\delta^{2}}{2}.$

We will now simplify this upper bound, so we can later make an appropriate
choice of $\alpha$ and $\eta$ leading to the desired complexity bound. We will
focus now on the geometric sum, and start by using the definition of
$\rho_{t}$ to write

$\eta\sum_{t=0}^{T-1}\rho_{t+1}=\eta\theta\sum_{t=0}^{T-1}\theta^{t}=\eta\frac{\theta}{1-\theta}(1-\theta^{T}).$

On the one hand, we see that from the definition of $\theta$ we have

$\eta\frac{\theta}{1-\theta}=\frac{\eta}{1+\eta\eta}\frac{1}{1-\frac{1}{1+\eta\eta}}=\frac{\eta}{\eta\eta}=\eta^{-1}.$

On the other hand, we can guarantee that $(1-\theta^{T})\geq 1/2$ provided
that $\eta=\frac{2}{\eta T}$. To see this, first observe that
$(1-\theta^{T})\geq 1/2$ is equivalent to
$T\geq\frac{\ln(2)}{\ln(\theta^{-1})}$. Since we impose that
$\eta\eta=\frac{2}{T}$ then necessarily
$\theta=\frac{T}{T+2}\Leftrightarrow\theta^{-1}=1+\frac{2}{T}$. So what we
need is $T\geq\frac{\ln(2)}{\ln(1+\frac{2}{T})}$. Now, observe that
$\ln(2)\leq 1$. Moreover, it is an exercise (see Lemma A.5 which is deferred
to the end of this proof) to verify that for all $x\geq 0$,
$\frac{1}{\ln(x+1)}\leq\frac{1}{2}+\frac{1}{x}$. So for our condition to hold
it is enough to ask that $T\geq\frac{1}{2}+\frac{T}{2}$, which is equivalent
to $T\geq 1$, which is true.

Combining all those bounds together with our new definition for $\eta$, we now
have

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{2\|\bm{w}_{0}-\bm{w}_{*}\|^{2}}{\eta
T}+\eta\frac{\sigma^{2}+(\delta+G)^{2}}{2}+\frac{\eta\delta^{2}T}{4}.$

To simplify the last stage of our analysis, let us note
$a=2\|\bm{w}_{0}-\bm{w}_{*}\|^{2}$, $b=\frac{\sigma^{2}+(\delta+G)^{2}}{2}$
and $c=\frac{\delta^{2}}{4}$, so that our bound writes as

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{a}{\eta
T}+\eta b+c\eta T.$ (24)

Minimizing the right-hand side is equivalent to solve $\frac{a}{\eta T}=\eta
b+c\eta T$, where

$\frac{a}{\eta T}=\eta b+c\eta
T\Leftrightarrow\eta^{2}=\frac{a}{T(b+cT)}\Leftrightarrow\eta=\sqrt{\frac{a}{bT+cT^{2}}}\Leftrightarrow\eta=2\sqrt{\frac{2\|\bm{w}_{0}-\bm{w}_{*}\|^{2}}{2(\sigma^{2}+(\delta+G)^{2})T+\delta^{2}T^{2}}}$

With such choice of stepsize, our bound becomes

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{2a}{\eta
T}=\frac{2a}{T}\sqrt{\frac{bT+cT^{2}}{a}}=2\sqrt{a}\sqrt{\frac{b}{T}+c}=2\|\bm{w}_{0}-\bm{w}_{*}\|\sqrt{\frac{\sigma^{2}+(\delta+G)^{2}}{T}+\frac{\delta^{2}}{2}}$

and this gives us our first bound. For the second bound, let us take simply
$\eta=\frac{1}{\delta T+\sqrt{T}}$ and inject it into (24) to obtain

$\displaystyle\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]$
$\displaystyle\leq$ $\displaystyle\frac{a\delta
T+a\sqrt{T}}{T}+\frac{b}{\delta T+\sqrt{T}}+\frac{cT}{\delta T+\sqrt{T}}$
$\displaystyle\leq$ $\displaystyle
a\delta+\frac{a}{\sqrt{T}}+\frac{b}{\sqrt{T}}+\frac{c}{\delta}=\frac{a+b}{\sqrt{T}}+a\delta+\frac{c}{\delta}$
$\displaystyle\leq$
$\displaystyle\frac{2D^{2}+\sigma^{2}+\delta^{2}+G^{2}}{\sqrt{T}}+(2D^{2}+1)\delta,$

where in the last inequality we simplified some numerical constants. ∎

######  Lemma A.5.

For every $x\geq 0$, $\frac{1}{\ln(1+x)}\leq\frac{1}{2}+\frac{1}{x}$.

###### Proof.

This inequality is equivalent to $\ln(1+x)\geq\frac{2x}{x+2}$, or again
$(x+2)\ln(1+x)\geq 2x$. Define $\phi:(-1,+\infty)\to\mathbb{R}$ as
$\phi(x)=(x+2)\ln(1+x)$ and compute its derivatives:

$\phi^{\prime}(x)=\ln(1+x)+1+\frac{1}{1+x}\quad\text{ and
}\quad\phi^{\prime\prime}(x)=\frac{x}{(1+x)^{2}}.$

We see that $\phi^{\prime\prime}(x)\geq 0$ for all $x\geq 0$, so $\phi$ is
convex on $[0,+\infty)$. So we can use the tangent inequality:

$\phi(x)\geq\phi(0)+\phi^{\prime}(0)(x-0)=0+2x=2x.$ (25)

Therefore $\phi(x)\geq 2x$ for all $x\geq 0$, which is what we wanted to
prove. ∎

###  A.4 Rates for SMGD for Convex Lipschitz Functions: Proof of Proposition
3.4

See 3.4

###### Proof.

In this proof, we note $\mathcal{F}_{t}$ the filtration generated by
$\bm{w}_{0},\dots,\bm{w}_{t}$, and we will note $\mathbb{E}_{t}$ to refer to
the expectation conditioned on $\mathcal{F}_{t}$.

We want to apply Theorem A.4. Let us denote
$\bm{g}_{t}=\bm{m}_{t}+\bm{e}_{t}$, so that the iterates of SMGD verify
$\bm{w}_{t+1}=\bm{w}_{t}-\eta\bm{g}_{t}$. We are now going to show that
$\bm{g}_{t}$ has uniformly bounded bias and variance. For this we will make
use of Assumptions 3.1 and 3.2. Remember that we assumed
$\delta_{2}=\sigma_{2}=0$ without loss of generality (we can do this according
to Lemma A.1), so we note $\delta$ and $\sigma$ instead of $\delta_{1}$ and
$\sigma_{1}$. Also, Assumption 3.2 on the errors together with the tower rule
imply that
$\mathbb{E}_{t}\left[\|\bm{e}_{t}\|^{2}\right]=\mathbb{E}_{t}\left[\mathbb{E}\left[\|\bm{e}_{t}\|^{2}\
|\ \bm{g}_{t}^{(1)},\dots,\bm{g}_{t}^{(n)}\right]\right]\leq\varepsilon^{2}$.
Regarding the bias, we can write

$\|\mathbb{E}_{t}\left[\bm{g}_{t}-\nabla\ell(\bm{w}_{t})\right]\|\leq\|\mathbb{E}_{t}\left[\bm{m}_{t}-\nabla\ell(\bm{w}_{t})\right]\|+\|\mathbb{E}_{t}\left[\bm{e}_{t}\right]\|\leq\delta+\varepsilon,$

where in the last inequality we used Jensen’s inequality to write
$\|\mathbb{E}_{t}\left[\bm{e}_{t}\right]\|\leq\sqrt{\mathbb{E}_{t}\left[\|\bm{e}_{t}\|^{2}\right]}$.
As for the variance, we write

$\mathbb{E}_{t}\left[\|\bm{g}_{t}-\mathbb{E}_{t}\bm{g}_{t}\|^{2}\right]\leq
2\mathbb{E}_{t}\left[\|\bm{m}_{t}-\mathbb{E}_{t}\left[\bm{m}_{t}\right]\|^{2}\right]+2\mathbb{E}_{t}\left[\|\bm{e}_{t}-\mathbb{E}_{t}\left[\bm{e}_{t}\right]\|^{2}\right]\leq
2\sigma^{2}+8\mathbb{E}_{t}\left[\|\bm{e}_{t}\|^{2}\right]\leq
2\sigma^{2}+8\varepsilon^{2}.$

Using Theorem A.4, we obtain

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{2\|\bm{w}_{0}-\bm{w}_{*}\|^{2}}{\eta
T}+\eta\frac{2\sigma^{2}+8\varepsilon^{2}+(\delta+\varepsilon+G)^{2}}{2}+\frac{\eta(\delta+\varepsilon)^{2}T}{4}.$

Setting the step size $\eta=\frac{1}{(\delta+\varepsilon)T+\sqrt{T}}$, we
conclude that

$\mathbb{E}\left[\ell(\bar{\bm{w}}_{T})-\inf\ell\right]\leq\frac{2D^{2}+2\sigma^{2}+8\varepsilon^{2}+(\delta+\varepsilon)^{2}+G^{2}}{\sqrt{T}}+(2D^{2}+1)(\delta+\varepsilon).$

∎

###  A.5 Rates for SMGD with Stable Noise: Proof of Corollary 3.6

See 3.6

###### Proof.

First of all, remember that the $\ell_{1}$-median can be computed exactly, by
computing component-wise a one-dimensional median. This allows us to have
$\varepsilon=0$ in Assumption 3.2.

Proof of part (i). The result is a direct consequence of Lemma A.3 and
Proposition 3.3.

Proof of part (ii). The result is a direct consequence of Lemma A.3 and
Proposition 3.4.

This completes the proof. ∎

##  Appendix B Approximate Computation of Median-like Estimators

###  B.1 Median-like Estimators : General Considerations

Given a function $\mathcal{D}:\mathbb{R}^{d}\to\mathbb{R}$ we consider the
problem (4), which writes as

$\underset{\bm{m}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathbb{E}\left[\mathcal{D}(\bm{m}-\bm{z})\right].$

In what follows, we will mostly consider functions $\mathcal{D}$ satisfying
the following assumption.

######  Assumption B.1.

The function $\mathcal{D}:\mathbb{R}^{d}\to\mathbb{R}$ is convex,
${\rm{argmin}}~{}\mathcal{D}=\\{0\\}$ and the problem (4) admits a solution.

Note that assuming $\mathcal{D}$ to be convex and finite implies that it is
continuous, therefore measurable, which means that
$\mathbb{E}\left[\mathcal{D}(\bm{m}-\bm{z})\right]$ is well defined. The
assumption that ${\rm{argmin}}~{}\mathcal{D}=\\{0\\}$ guarantees that in the
trivial case when the distribution is a Dirac $\delta_{\bm{m}_{0}}$, then the
solution to (4) is exactly $\bm{m}_{0}$. The assumption that a solution to (4)
exists is easily verified in practice, under the assumption that $\mathcal{D}$
is coercive:

######  Lemma B.2.

Let $\mathcal{D}:\mathbb{R}^{d}\to\mathbb{R}$ be convex. Assume that
$\mathbb{E}\left[\mathcal{D}(\bm{z}-\mathbb{E}\left[\bm{z}\right])\right]<+\infty$.
If $\mathcal{D}$ is coercive, i.e.
$\lim\limits_{\|m\|\to\infty}\mathcal{D}(\bm{m})=+\infty$, then (4) admits a
solution.

###### Proof.

Let $\phi(\bm{m}):=\mathbb{E}\left[\mathcal{D}(\bm{m}-\bm{z})\right]$. Since
we assume that
$\phi(\mathbb{E}\left[\bm{z}\right])=\mathbb{E}\left[\mathcal{D}(\bm{z}-\mathbb{E}\left[\bm{z}\right])\right]<+\infty$
then $\phi$ is proper. Since $\mathcal{D}$ is convex then $\phi$ is convex, as
an expectation of convex functions. Moreover $\mathcal{D}$ is convex and
finite, which means that $\mathcal{D}$ is continuous on $\mathbb{R}^{d}$, see
e.g. Corollary 8.39 in (Bauschke & Combettes, 2017). Therefore $\phi$ is an
expectation of lower semicontinuous functions, which implies that $\phi$ is
lower semicontinuous.

Since $\mathcal{D}$ is coercive, convex and continuous, we know that there
exists $a\in(0,+\infty)$ and $b\in\mathbb{R}$ (see Proposition 14.16 in
(Bauschke & Combettes, 2017)) such that $\mathcal{D}(\cdot)\geq a\|\cdot\|+b$.
Therefore

$\displaystyle\phi(\bm{m})$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathcal{D}(\bm{m}-\bm{z})\right]\geq
a\mathbb{E}\left[\|\bm{m}-\bm{z}\|\right]+b$ $\displaystyle\geq$
$\displaystyle
a\|\bm{m}-\mathbb{E}\left[\bm{z}\right]\|-a\mathbb{E}\left[\|\bm{z}-\mathbb{E}\left[\bm{z}\right]\|\right]+b$
$\displaystyle\geq$ $\displaystyle
a\|\bm{m}-\mathbb{E}\left[\bm{z}\right]\|-a\mathbb{E}\left[\mathcal{D}(\bm{z}-\mathbb{E}\left[\bm{z}\right])\right]\underset{\|\bm{m}\|\to\infty}{\longrightarrow}+\infty,$

where in the inequalities we used the fact that $\|\cdot\|$ is $1$-Lipschitz,
and that $a\mathbb{E}\left[\|\bm{z}-\mathbb{E}\left[\bm{z}\right]\|\right]\leq
a\mathbb{E}\left[\mathcal{D}(\bm{z}-\mathbb{E}\left[\bm{z}\right])\right]-b<+\infty$.
So $\phi$ is coercive on top of being proper convex and lower semicontinuous.
We can then conclude that $\phi$ has a minimizer, with for instance
Proposition 11.15 from (Bauschke & Combettes, 2017). ∎

###  B.2 Approximations based on Subgradient Updates

Instead of SPP, we could solve (8) by iteratively taking steps of stochastic
subgradient descent. If $\mathcal{D}$ is convex, then in each iteration we
sample a gradient $\bm{g}_{t}$, and update our current estimate $\bm{m}_{t}$
by

$\displaystyle\bm{m}_{t+1}=\bm{m}_{t}-\tau\bm{u}_{t},\quad\bm{u}_{t}\in\partial\mathcal{D}(\bm{m}_{t}-\bm{g}_{t}).$
(26)

Here, $\tau>0$ is the learning rate, and $\bm{u}_{t}$ is a (convex)
subgradient. Since $\mathcal{D}$ could be a non-differentiable function such
as the $\ell_{2}$\- or $\ell_{1}$-norm, we use subgradients instead of
gradients.

#### Squared $\ell_{2}$-norm.

As a first simple example of (26), let
$\mathcal{D}=\tfrac{1}{2}\|\cdot\|_{2}^{2}.$ In this case (26) becomes

$\displaystyle\bm{m}_{t+1}=\bm{m}_{t}-\tau(\bm{m}_{t}-\bm{g}_{t})=(1-\tau)\bm{m}_{t}+\tau\bm{g}_{t}.$
(27)

This is again (heavy-ball) momentum (cf. (14)), now with coefficient
$\beta=1-\tau$. Note that in contrast to SPP, the coefficient $\beta$ is only
in $[0,1)$ (as typically is the case for momentum) if $\tau\in[0,\;1]$.

####  $\ell_{1}$-norm.

If we choose $\mathcal{D}=\|\cdot\|_{1}$, update (26) gives

$\displaystyle\bm{m}_{t+1}=\bm{m}_{t}+\tau\operatorname{sgn}(\bm{g}_{t}-\bm{m}_{t}),$

where the $\operatorname{sgn}$-operator can take any value in $[-1,1]$ for
coordinates $i\in[d]$ such that $(\bm{g}_{t}-\bm{m}_{t})_{i}=0$. This recovers
a version of sign-SGD (Riedmiller & Braun, 1993; Bernstein et al., 2018),
where the $\mathrm{sign}$-operation is applied to the increment
$\bm{m}_{t}-\bm{g}_{t}$. For example, if we reset $\bm{m}_{t}=0$ in every
iteration, we obtain exactly sign-SGD.

####  $\ell_{2}$-norm.

If we choose $\mathcal{D}=\|\cdot\|_{2}$, and if $\bm{g}_{t}\neq\bm{m}_{t}$,
update (26) gives

$\displaystyle\bm{m}_{t+1}=\bm{m}_{t}-\tau\frac{\bm{m}_{t}-\bm{g}_{t}}{\|\bm{m}_{t}-\bm{g}_{t}\|}=\big{(}1-\frac{\tau}{\|\bm{m}_{t}-\bm{g}_{t}\|}\big{)}\bm{m}_{t}+\frac{\tau}{\|\bm{m}_{t}-\bm{g}_{t}\|}\bm{g}_{t}.$

Otherwise, if $\bm{g}_{t}=\bm{m}_{t}$, we can choose any $\bm{u}_{t}$ with
$\|\bm{u}_{t}\|\leq 1$ and set $\bm{m}_{t+1}=\bm{m}_{t}-\tau\bm{u}_{t}$; in
particular we can set $\bm{m}_{t+1}=\bm{m}_{t}$.

###  B.3 Approximations based on the Stochastic Proximal Point: General
$\mathcal{D}$

Given a proper, closed convex function
$\phi:\mathbb{R}^{d}\to\mathbb{R}\cup\\{+\infty\\}$ and $\tau>0$, we recall
the definition of the proximal operator

$\mathrm{prox}_{\tau\phi}(\bm{m}_{0})=\underset{\bm{m}\in\mathbb{R}^{d}}{\rm{argmin}}~{}\phi(\bm{m})+\frac{1}{2\tau}\|\bm{m}-\bm{m}_{0}\|^{2}.$

Applying the SPP algorithm to (8) gives the following iterations

$\begin{cases}\text{Sample $\bm{g}_{t}$ i.i.d. from $\mathcal{G}$}\\\
\bm{m}_{t+1}=\mathrm{prox}_{\tau\mathcal{D}(\cdot-\bm{g}_{t})}(\bm{m}_{t}).\end{cases}$
(28)

As we see next, those iterations can be reformulated to simply involve the
proximal operator of $\mathcal{D}$:

######  Proposition B.3.

Let $\mathcal{D}$ verify Assumption B.1. The SPP update (28) for solving (8)
can be equivalently written as

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{g}_{t}+\mathrm{prox}_{\tau\mathcal{D}}(\bm{m}_{t}-\bm{g}_{t})$
(29)
$\displaystyle=\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\bm{m}_{t}-\cdot}{\tau})}(\bm{g}_{t}),$
(30)

where
$\mathcal{D}^{*}(\bm{m})=\sup_{\bm{y}}\left<\bm{y},\bm{m}\right>-\mathcal{D}(\bm{y})$
is the Fenchel conjugate of $\mathcal{D}$, and $\bm{g}_{t}$ is sampled i.i.d.
from $\mathcal{G}$ at each iteration.

###### Proof.

Here we make use of proximal calculus rules, which can found in most textbooks
on convex analysis, such as chapter 6 in (Beck, 2017). Starting from (9), with
$\bm{y}=\bar{\bm{y}}+\bm{g}_{t}$, we have

$\displaystyle\begin{split}&\bm{m}_{t+1}=\underset{\bm{y}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathcal{D}(\bm{y}-\bm{g}_{t})+\tfrac{1}{2\tau}\|\bm{y}-\bm{m}_{t}\|^{2}\\\
&\hskip
4.30554pt=\bm{g}_{t}+\underset{\bar{\bm{y}}\in\mathbb{R}^{d}}{\operatorname{argmin}}\;\mathcal{D}(\bar{\bm{y}})+\tfrac{1}{2\tau}\|\bar{\bm{y}}-(\bm{m}_{t}-\bm{g}_{t})\|^{2}\\\
&\hskip
4.30554pt=\bm{g}_{t}+\mathrm{prox}_{\tau\mathcal{D}}(\bm{m}_{t}-\bm{g}_{t}).\end{split}$

This proves (29). As for (30), we combine
$(\tau\mathcal{D})^{*}=\tau\mathcal{D}^{*}(\cdot/\tau)$ (Beck, 2017, Thm.
4.14) and the Moreau decomposition theorem (Beck, 2017, Thm. 6.44), to obtain

$\mathrm{prox}_{\tau\mathcal{D}}(\bm{m})=\bm{m}-\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\cdot}{\tau})}(\bm{m})$

in (29) gives

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{g}_{t}+\mathrm{prox}_{\tau\mathcal{D}}(\bm{m}_{t}-\bm{g}_{t})$
$\displaystyle=\bm{g}_{t}+\bm{m}_{t}-\bm{g}_{t}-\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\cdot}{\tau})}(\bm{m}_{t}-\bm{g}_{t})$
$\displaystyle=\bm{m}_{t}-\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\cdot}{\tau})}(\bm{m}_{t}-\bm{g}_{t})$
(31)

Finally, using that for any function $f$ it holds
$\underset{\bm{y}}{\operatorname{argmin}}\;f(\bm{y})=\bm{m}+\underset{\bar{\bm{y}}}{\operatorname{argmin}}\;f(\bar{\bm{y}}+\bm{m})$
in (31) gives

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{m}_{t}-\left(\bm{m}_{t}+\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\cdot+\bm{m}_{t}}{\tau})}(-\bm{g}_{t})\right)$
$\displaystyle=-\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\cdot+\bm{m}_{t}}{\tau})}(-\bm{g}_{t})$
$\displaystyle=\mathrm{prox}_{\tau\mathcal{D}^{*}(\frac{\bm{m}_{t}\cdot-}{\tau})}(\bm{g}_{t})$

where in the final equality we used that if $f(\bm{x})=g(-\bm{x})$ we have
that

$\mathrm{prox}_{g}(\bm{m})=-\mathrm{prox}_{f}(-\bm{m}).$

∎

###  B.4 Approximations based on the Stochastic Proximal Point: Particular
Cases of $\mathcal{D}$

####  $\ell_{p}$-norm.

We now specialize the previous results for when $\mathcal{D}$ is an
$\ell_{p}$-norm. See 4.1

###### Proof.

The proximal operator of the $\ell_{p}$-norm is given by

$\mathrm{prox}_{\tau\|\cdot\|_{p}}(\bm{m})=\bm{m}-\tau\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0})}(\bm{m}/\tau),$

see Example 6.47 in (Beck, 2017). Using this in (29) gives

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{g}_{t}+\mathrm{prox}_{\tau\|\cdot\|_{p}}(\bm{m}_{t}-\bm{g}_{t})$
$\displaystyle=\bm{m}_{t}-\tau\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},1)}\big{(}\frac{\bm{m}_{t}-\bm{g}_{t}}{\tau}\big{)}$
$\displaystyle=\bm{m}_{t}+\tau\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},1)}\big{(}\frac{\bm{g}_{t}-\bm{m}_{t}}{\tau}\big{)}$
$\displaystyle=\bm{m}_{t}+\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},\tau)}\big{(}\bm{g}_{t}-\bm{m}_{t}\big{)}$
$\displaystyle=\bm{m}_{t}+\left(\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m}_{t},\tau)}\big{(}\bm{g}_{t}\big{)}-\bm{m}_{t}\right)$
$\displaystyle=\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m}_{t},\tau)}\big{(}\bm{g}_{t}\big{)}$

where the third equality we used that

$\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},1)(\bm{m})}=-\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},1)(-\bm{m})},$

followed by

$\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m}/\tau,1)}(\bm{y}/\tau)=\frac{1}{\tau}\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m},\tau)}(\bm{y})$

in the fourth equality and

$\mathrm{Proj}_{\mathbb{B}_{q}(\bm{0},\tau)}(\bm{y})=\mathrm{Proj}_{\mathbb{B}_{q}(\bm{m},\tau)}(\bm{y}-\bm{m})-\bm{m}$

in the fifth equality. ∎

####  $\ell_{1}$-norm.

See 4.3

###### Proof.

Follows from (10) and that

$\mathrm{Proj}_{\mathbb{B}_{2}(\bm{0},\tau)}\big{(}\bm{v}\big{)}=\frac{\tau\;\bm{v}}{\max\\{\tau,\|\bm{v}\|_{2}\\}}.$

∎

####  $\ell_{2}$-norm.

See 4.2

###### Proof.

Follows directly from (10), by using that

$\;\min\\{\max\\{\bm{v},-\tau\\},\tau\\}=\mathrm{Proj}_{\mathbb{B}_{\infty}(\bm{0},\tau)}\big{(}\bm{v}\big{)}.$

∎

#### Huber function.

The Huber function is given by

$\displaystyle H_{\mu}:\mathbb{R}^{d}\to\mathbb{R},\quad
H_{\mu}(\bm{z})=\begin{cases}\frac{1}{2\mu}\|\bm{z}\|^{2}&\quad\|\bm{z}\|\leq\mu,\\\
\quad\|\bm{z}\|-\frac{\mu}{2}&\quad\text{else}.\end{cases}$

It operates like the squared loss for arguments with small norm, and like the
$\ell_{2}$-norm for arguments with large norms (potentially outliers).

######  Lemma B.4.

Update (9) with $\mathcal{D}=H_{\mu}$ is given by

$\displaystyle\bm{m}_{t+1}$
$\displaystyle=\bm{g}_{t}+\big{(}1-\tfrac{\tau}{\max\\{\|\bm{m}_{t}-\bm{g}_{t}\|,\tau+\mu\\}}\big{)}(\bm{m}_{t}-\bm{g}_{t})=\big{(}1-\tfrac{\tau}{\max\\{\|\bm{m}_{t}-\bm{g}_{t}\|,\tau+\mu\\}}\big{)}\bm{m}_{t}+\tfrac{\tau}{\max\\{\|\bm{m}_{t}-\bm{g}_{t}\|,\tau+\mu\\}}\bm{g}_{t}.$
(32)

###### Proof.

Follows from (9) and the fact that the proximal operator of the Huber function
(see (Beck, 2017, Example 6.66)) is given by $\mathrm{prox}_{\tau
H_{\mu}}(\bm{z})=\Big{(}1-\frac{\tau}{\max\\{\|\bm{z}\|,\tau+\mu\\}}\Big{)}\bm{z}.$
∎

The estimator (32) has the advantage that it is a relaxation of vectorwise
clipping (12), as
$1-\frac{\tau}{\max\\{\|\bm{m}_{t}-\bm{g}_{t}\|,\tau+\mu\\}}\geq
1-\frac{\tau}{\tau+\mu}>0$. However, it comes with the disadvantage that it
has two hyperparameters $\tau,~{}\mu$ instead of one.

###  B.5 Proof of Proposition 4.4: Convergence when $\mathcal{D}$ is a Norm

See 4.4

###### Proof.

Since we assume that $\mathcal{D}$ is a norm over $\mathbb{R}^{d}$ then it
convex and continuous. Moreover, since all norms are equivalent in
$\mathbb{R}^{d}$, we have that $\mathcal{D}(\bm{w})\leq C\|\bm{w}\|_{2}$ for
some constant $G>0$. Thus $\mathcal{D}$ is a $G$-Lipschitz function. If now we
write $\varphi_{\bm{g}}:=\mathcal{D}(\cdot-\bm{g})$, it is clear that it is
convex and Lipschitz. Therefore,
$\varphi=\mathbb{E}_{\bm{g}}\left[\varphi_{\bm{g}}\right]$ is itself a convex
$G$-Lipschitz function. We can now conclude with Theorem 4.4 from (Davis &
Drusvyatskiy, 2019). ∎

Here is a corollary stating rates on the expected distance between the
approximate sample median and the true sample median, when using the
$\ell_{1}$-norm:

######  Corollary B.5.

In the context of Proposition 4.4, assume that $\mathcal{D}=\|\cdot\|_{1}$ and
that $\mathcal{G}$ is a sum of Dirac measures
$\tfrac{1}{n}\sum_{i=1}^{n}\delta_{\bm{g}^{(i)}}$. Then we have further

$\mathbb{E}\left[\|\bar{\bm{m}}_{T}-\bm{m}_{*}\|\right]=O\left(\frac{\|\bm{m}_{0}-\bm{m}_{*}\|^{2}}{\tau
T}+\tau\right).$

###### Proof.

In this proof, we use the notation $(\bm{x}_{j})_{j=1,\ldots,d}$ to specify
the indivdual components of a vector $\bm{x}\in\mathbb{R}^{d}$. For
$\bm{x}\in\mathbb{R}^{d}$, we define the sign operator as follows: let
$\operatorname{sgn}(\bm{x})=(\operatorname{sgn}(\bm{x}_{j}))_{j=1,\ldots,d}$,
where $\bm{x}_{j}$ is the $j$-th coordinate of $\bm{x}\in\mathbb{R}^{d}$, and
where we define for a scalar $t\in\mathbb{R}$ the set-valued operator
$\operatorname{sgn}(t)=\\{+1\\}$ if $t>0$, $\operatorname{sgn}(t)=\\{-1\\}$ if
$t<0$, and $\operatorname{sgn}(t)=[-1,1]$ if $t=0$. Further, given a set
$A\subset\mathbb{R}^{d}$ we define $\|A\|_{-}:=\inf\limits_{a\in A}\|a\|$.

To obtain the desired bound, and given that we already have the conclusion of
Proposition 4.4, it is enough for us to show that there exists $\mu>0$ such
that

$(\forall\bm{m}\in\mathbb{R}^{d})\quad\mu\operatorname{dist}(\bm{m},{\rm{argmin}}~{}\varphi)\leq\varphi(\bm{m})-\inf\varphi.$

It is a standard result from variational analysis that this property555It is
sometimes referred to as (weak) sharp minima (Ferris, 1991; Burke & Ferris,
1993), (superlinear)(linear) conditioning (Lemaire, 1992; Cornejo et al.,
1997; Lemaire, 1998), or error bound (Lewis & Pang, 1998). is true if and only
if the following inequality is satisfied (see Theorem 5.1 in (Cornejo et al.,
1997) or Proposition 3.1 in (Lemaire, 1998))

$(\forall\bm{m}\in\operatorname{dom}\partial\varphi)\quad\|\partial\varphi(\bm{m})\|_{-}\geq\mu.$

Our goal now is to prove the above inequality, from which the conclusion will
follow. Let us consider $\bm{m}\notin{\rm{argmin}}~{}\varphi$, and use
standard calculus to write

$\partial\varphi(\bm{m})=\frac{1}{n}\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}-{\bm{g}}_{i})=\frac{1}{n}\sum_{i=1}^{n}\big{(}\operatorname{sgn}(\bm{m}_{j}-({\bm{g}}_{i})_{j})\big{)}_{j=1,\ldots,d}=\frac{1}{n}\Big{(}\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}_{j}-({\bm{g}}_{i})_{j})\Big{)}_{j=1,\ldots,d}.$

Our goal is to show that $\|\partial\varphi(\bm{m})\|_{-}\geq\frac{1}{n}$,
which is equivalent to show that
$\|\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}-{\bm{g}}_{i})\|_{-}\geq 1$. Because
$\bm{m}\notin{\rm{argmin}}~{}\varphi$, we know that
$0\notin\partial\varphi(\bm{m})$ so there must exist a coordinate $j$ such
that the $j$-th coordinate of $\partial\varphi(\bm{m})$ does not contain zero,
in other words such that
$0\notin\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}_{j}-({\bm{g}}_{i})_{j})$.
Using the fact that $\|\cdot\|_{2}\geq\|\cdot\|_{\infty}$, we can therefore
lower bound
$\|\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}-{\bm{g}}_{i})\|_{-}\geq|\sum_{i=1}^{n}\operatorname{sgn}(\bm{m}_{j}-({\bm{g}}_{i})_{j})|_{-}$.
Let us denote $s_{i}:=\operatorname{sgn}(\bm{m}_{j}-({\bm{g}}_{i})_{j})$, so
that we want to prove that $|\sum_{i=1}^{n}s_{i}|_{-}\geq 1$. Let us now
consider a few cases.

  * •

If $s_{i}=\\{+1\\}$ or $\\{-1\\}$ for every $i$, then $\sum_{i=1}^{n}s_{i}$ is
a singleton. Furthermore it is a sum of relative numbers, so
$\sum_{i=1}^{n}s_{i}\in\mathbb{Z}$. But we also know that
$0\notin\sum_{i=1}^{n}s_{i}$, so we conclude that
$\sum_{i=1}^{n}s_{i}\in\mathbb{Z}^{*}$, and so that
$|\sum_{i=1}^{n}s_{i}|_{-}\geq 1$.

  * •

If $s_{i}=[-1,+1]$ for every $i$, then we immediately see that
$0\in\sum_{i=1}^{n}s_{i}$ which is a contradiction.

  * •

Otherwise, the $s_{i}$ are combinations of singletons and intervals. Let us
note $I=\\{i:s_{i}=[-1,+1]\\}$ and $I^{\prime}=\\{i:s_{i}=\\{\pm 1\\}\\}$,
which are not empty by assumption. Let us also note $k\geq 1$ the cardinality
of $I$, so that $\sum_{i\in I}s_{i}=[-k,+k]$. We also introduce $s:=\sum_{i\in
I^{\prime}}s_{i}\in\mathbb{Z}$. Now we can write

$|\sum_{i=1}^{n}s_{i}|_{-}=\inf\limits_{t\in\sum_{i=1}^{n}s_{i}}|t|=\inf\limits_{t\in[-k,+k]}|t+s|=d(s;[-k,+k]),$

where the latter is the distance from the number $s$ to the interval
$[-k,+k]$. Now remember that $0\notin\sum_{i=1}^{n}s_{i}$, so this distance
cannot be zero, which means that $|s|>k$. Moreover both $s$ and $k$ are
rational numbers, so this distance must be greater or equal to $1$. This
concludes the proof.

∎

##  Appendix C Supplementary Information on Experiments

This section provides additional information for the experimental setup.

For easier interpretation, we use the first order approximation of momentum
(14), namely $\bm{m}_{t+1}=\left(1-\tau\right)\bm{m}_{t}+\tau\bm{g}_{t}$. That
is, we replace $\frac{\tau}{1+\tau}$ with its first-order Taylor approximation
around zero, which is equal to $\tau$. Note that $\tau$ is typically small,
for example for $\tau=0.01$ we have $\frac{\tau}{1+\tau}\approx 0.0099$. This
applies only to the experiments with estimating a fixed vector (or gradient)
in Section 6.1.

For the language modeling experiments, all details are identical to (Kunstner
et al., 2023), Section A.1. They also provide implementations for all tasks at
https://github.com/fKunstner/noise-sgd-adam-sign, which we use. For the
training runs for language modeling, we use batch size $256$ for PTB, $320$
for WikiText-2, and $32$ for SQuAD.

###  C.1 Additional Plots

Here, we show the results for estimating the full gradient of a transformer
architecture, with weights fixed at initialization. Fig. 7 shows the results
for PTB where we tried three different batch sizes $\\{64,256,1024\\}$. Fig. 8
shows the results for WikiText-2.

_Discussion._ From Fig. 7, we observe two phenomena: in the long run, momentum
attains the lowest error for estimating the full batch gradient. However, the
initial decrease of the error is much faster for CClip, followed by VClip.
This is important when using these estimates within a training setup such as
(17), where we only do one iteration of gradient estimation, followed by an
update of the weight (and hence a change in the full-batch gradient).
Secondly, we observe that the difference in convergence speeds is most
pronounced when the batch size increases. Hence, being robust to outliers
seems not to be solvable only by increasing the batch size and thus decreasing
the noise of the mini-batch gradient. In fact, the contrary seems to be the
case. This is similar to the observations made in (Kunstner et al., 2023).

Figure 7: Encoder transformer on PTB dataset, with weights fixed at
initialization. We use momentum, $\mathrm{clip}_{\tau,1}$ and
$\mathrm{clip}_{\tau,2}$ to estimate the full-batch gradient. Figure 8:
Encoder transformer on Wikitext-2 dataset, with weights fixed at
initialization. We use momentum, $\mathrm{clip}_{\tau,1}$ and
$\mathrm{clip}_{\tau,2}$ to estimate the full-batch gradient.

###  C.2 Learning Rate Values for Transformer Training

Here, we report for each of the language modelling tasks, the learning rate
value that is displayed in Fig. 5. For each method, we tuned the learning rate
over a set of values $10^{\frac{j}{2}}$ for a suitable interval of integer
numbers $j\in\mathbb{Z}$.

Table 1: Tuned learning rate values used for each method. All methods use constant learning rates. For SGD-M, tuning information is also reported in Section C.1 in (Kunstner et al., 2023)). Name | SGD-M | clipped-SGD | VClip | CClip | Adam
---|---|---|---|---|---
PTB | 0.1 | 0.316 | 3.16 | 0.316 | 0.001
WikiText-2 | 0.0316 | 0.316 | 1.0 | 0.316 | 0.001
SQuAD | 0.0316 | 0.1 | 0.316 | 0.1 | 0.000316