<figure><img src="image_1.png" /></figure>

HAL
open science

# Tukey-Median of Means-Gradients for Langevin Dynamics: a Robust Approach to Fitting Machine Learning Models

Stéphane Chrétien, Ben Gao, Jordan Patracone, Abdel-Rahim Mezidi, Olivier Alata

<figure><img src="image_2.png" /></figure>

## To cite this version:

Stéphane Chrétien, Ben Gao, Jordan Patracone, Abdel-Rahim Mezidi, Olivier Alata. Tukey-Median of Means-Gradients for Langevin Dynamics: a Robust Approach to Fitting Machine Learning Models. 2025. [https://arxiv.org/abs/2503.05319](https://arxiv.org/abs/2503.05319)

HAL Id: hal-05319256

https://inria.hal.science/hal-05319256v1

Preprint submitted on 17 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tukey-Median of Means-Gradients for Langevin Dynamics: a Robust Approach to Fitting Machine Learning Models

Stéphane Chrétien

Univ. Lyon2
Lab. ERIC

Ben Gao

Univ. Jean Monnet
Lab. Hubert Curien
Inria

Jordan Patracone

Univ. Jean Monnet
Lab. Hubert Curien
Inria

Abdel-Rahim Mezidi

Univ. Jean Monnet
Lab. Hubert Curien
Inria

Olivier Alata

Univ. Jean Monnet
Lab. Hubert Curien

## Abstract

Median-of-Means (MoM) estimators have emerged as powerful tools for robust learning when data are corrupted by outliers or exhibit heavy-tailed distributions. However, incorporating MoM principles to design robust loss functions typically renders the underlying optimization problem highly non-convex, posing significant challenges for the design of efficient, polynomial-time algorithms. In this work, we apply MoM principle at the gradient level using Tukey median to robustly estimate the uncorrupted gradient. More precisely, we propose a novel Tukey MoM-based variant of Stochastic Gradient Langevin Dynamics (TMoM-SGLD) which enjoys strong theoretical guarantees, including hitting-time bounds showing convergence to approximate stationary points of the underlying uncorrupted loss. Contrary to prior MoM estimators, our results hold in broad non-convex settings, making TMoM-SGLD particularly relevant to modern machine learning tasks. Experiments on corrupted linear regression tasks demonstrate that, unlike prior approaches, TMoM-SGLD remains effective even under high levels of outlier contamination. Additionally, we further validate our approach on a deep learning task, showing its capabilities to robustly train neural networks.

## 1 INTRODUCTION

Learning algorithms typically assume that training data is clean and independently drawn from a fixed distribution. However, real-world datasets often violate this assumption (they may include outliers or

corrupted labels) and may exhibit heavy-tailed distributions. These issues can significantly degrade the performance of standard stochastic optimization algorithms, sparking a growing interest in robust learning. A central objective in robust learning is to design robust estimators and robust stochastic optimization algorithms that remain statistically efficient while being resilient to corrupted or heavy-tailed data.

In response, the machine learning community has proposed several algorithmic strategies to enhance robustness. Some methods directly modify the optimization dynamics, for example by robustifying gradients (Prasad et al., 2020; Merad and Gaïffas, 2024) or momentum terms (Ilboudo et al., 2022). Others rely on heuristic adjustments such as *reweighting schemes* (Ren et al., 2018) or *curriculum learning approaches* (Jiang et al., 2018) to name a few. While these techniques often yield empirical gains, they either do not provide principled guarantees or do not exploit the underlying geometrical structure of the learning problem.

This limitation has motivated a complementary line of research, grounded in robust statistics, that seeks theoretical performance guarantees under contamination. More precisely, a prominent line of work has focused on obtaining sub-Gaussian estimators even in the presence of adversarial contamination (Audibert et al., 2011; Lugosi and Mendelson, 2016; Minsker, 2018; Lerasle, 2019).

Among the most promising tools are *Median-of-Means* (MoM) estimators. Developed in the robust statistics literature and stemming from a series of powerful ideas initiated by Nemirovsky and Yudin (1983), they have recently come to the forefront thanks to the recent theoretical achievements and their improved applicability (Kwon et al., 2018; Lecué and Lerasle, 2019; Chinot et al., 2019; Lecué and Lerasle, 2020). MoM estimators offer provable robustness to outliers and heavy-tailed noise without requiring strong assumptions on the data distribution. While effective in theory, incorporating MoM principles into optimization
pipelines introduces significant challenges. In practice, data are reshuffled at each epoch and are partitioned into blocks, then MoM principles are applied either to loss values (Lecué and Lerasle, 2020; Huang and Lederer, 2023) or to the gradients (Prasad et al., 2020; Puchkin et al., 2024). While Lecué and Lerasle (2020) and Puchkin et al. (2024) have provided theoretical guarantees for convex optimization problems, there still exist a gap to handle the non-convex case, notably encountered in deep learning tasks.

To address this, we introduce the Tukey Median-of-Means Stochastic Gradient Langevin Dynamics (TMoM-SGLD), a variant of SGLD, in which the gradient is replaced by the Tukey median (Tukey, 1975; Dyckerhoff and Mozharovskyi, 2016) computed over batch-mean gradients. Our proposal leverages two insights. First, from the robust statistics viewpoint, the Tukey median has shown promising robustness to high level of contamination (Donoho and Gasko, 1992), but its usage for learning tasks has not yet been investigated. Second, from the non-convex optimization viewpoint, SGLD permits to escape spurious local minima and benefits from polynomial-time guarantees to reach approximate stationary points (Zhang et al., 2017; Chen et al., 2020). By combining them together, our main contributions are as follows:

* We propose the Tukey-MoM as a robust gradient estimator supported by theoretical guarantees.

* We establish probabilistic theoretical guarantees in the form of hitting-time bounds, showing that TMoM-SGLD reaches approximate stationary points of the uncorrupted objective in polynomial time under non-convex settings in presence of data contaminations.

* We demonstrate that our method remains robust for higher level contaminations and heavy-tailed data, where existing baselines fall short.

Finally, our method complements and extends classical approaches in robust statistics, such as Catoni-type estimators (Catoni, 2012) and minimax estimators (Baud et al., 2017), by integrating their robustness insights into modern non-convex optimization pipelines.

In the following, we recall the robust learning setting followed by a literature review on MoM estimators in Section 3. Then, we present our contribution in detail in Section 4 and evaluate its empirical performance in Section 5. Finally, Section 6 concludes the paper with a discussion of limitations and perspectives.

## 2 PROBLEM STATEMENT

We consider the prototypical problem of supervised learning. Given a dataset of $n$ training examples

$\{X_i := (Z_i, Y_i)\}_{i=1}^n$, where each pair $(Z_i, Y_i)$ lies in $\mathbb{R}^d \times \mathbb{R}$, the goal is to learn a function that maps an input $Z_i \in \mathbb{R}^d$ to an output $Y_i \in \mathbb{R}$. The learning task is typically cast as an empirical risk minimization problem:

$$ \min_{\theta \in \Theta} \frac{1}{n} \sum_{i=1}^{n} \ell(Y_i, h_{\theta}(Z_i)), \quad (1) $$

where $h_{\theta}: \mathbb{R}^d \to \mathbb{R}$ is a parametric function indexed by a vector $\theta \in \Theta \subset \mathbb{R}^p$, with $\Theta$ being the set of admissible parameters, and $\ell: \mathbb{R} \times \mathbb{R} \to \mathbb{R}_+$ is a prescribed loss function. This formulation is general enough to cover most supervised learning tasks, including modern deep learning applications where $h_{\theta}$ is implemented as a neural network.

A standard assumption underpinning most popular learning algorithms is that the training samples are drawn i.i.d. from some unknown distribution over $\mathbb{R}^d \times \mathbb{R}$. However, in many real-world settings, this assumption is overly idealized: training data may be subject to heavy-tailed noise or contaminations that can impair learning. To better reflect such scenarios, we consider that the dataset is a mixture of clean and corrupted samples, though the identity of each is unknown. Formally, we assume the existence of an unknown partition of the index set:

$$ \{1, \dots, n\} = \mathcal{I}_{\text{uncorr}} \cup \mathcal{I}_{\text{corr}}, \quad \mathcal{I}_{\text{uncorr}} \cap \mathcal{I}_{\text{corr}} = \emptyset, \quad (2) $$

where $\mathcal{I}_{\text{uncorr}}$ indexes the clean examples and $\mathcal{I}_{\text{corr}}$ the corrupted ones. Ideally, one would like to minimize the empirical risk on the clean examples only, i.e.,

$$ \min_{\theta \in \Theta} \left\{ \hat{F}_{\text{uncorr}}(\theta) := \frac{1}{|\mathcal{I}_{\text{uncorr}}|} \sum_{i \in \mathcal{I}_{\text{uncorr}}} f_{\theta}(X_i) \right\}, \quad (3) $$

where the notation $f_{\theta}(X_i) := \ell(Y_i, h_{\theta}(Z_i))$ will be used as a shorthand hereafter. However, since the partition $\mathcal{I}_{\text{uncorr}} \cup \mathcal{I}_{\text{corr}}$ is unknown, directly optimizing this objective is infeasible. Standard empirical risk minimization over the full dataset may be heavily biased by even a small fraction of corrupted points. This motivates the need for robust learning strategies that can approximate the clean empirical risk without prior knowledge of which samples are corrupted. In the following, we recall one such approach based on the MoM principle, which offers provable robustness under minimal assumptions.

## 3 Robust Stochastic Optimization via the Median of Means Principle

Median of Means (MoM) estimators, originally introduced as robust mean estimators in the univariate setting (Donoho and Gasko, 1992), achieve sub-Gaussian deviation bounds while remaining resilient
to a small fraction of outliers. Formally, given n i.i.d. samples {$S_i$}$_{i=1}^n$ with $n = K b$, partitioned into K disjoint blocks $I_k = (k-1)b+1, \dots, kb$, the MoM estimator is defined as $MoM_K[S] = \text{median}(\frac{K}{n}\sum_{i \in I_1} S_i, \dots, \frac{K}{n}\sum_{i \in I_K} S_i)$.

Beyond mean estimation, the MoM principle has re-
cently been repurposed for robust supervised learning,
as described in Section 2. Two main strategies have
emerged: applying MoM at the level of loss values or
at the level of gradients.

Lecué and Lerasle (2020) applied the MoM princi-
ple at the loss values level using minmax strategy,
and achieved sub-Gaussian deviation bound in the
presence of a small portion of outliers. By partition-
ing {$X_i := (Z_i, Y_i)$}$_{i=1}^n$ into K equally sized disjoint
groups $I_1, \dots, I_K$, the minmax MoM estimator is given
by the solution of

$$
\min_{\theta \in \Theta} \max_{\theta' \in \Theta} \mathrm{MoM}_K[f_{\theta} - f_{\theta'}], \quad (4)
$$

where $\text{MoM}_K[f_\theta - f_{\theta'}] = \text{median}(\frac{K}{n}\sum_{i \in I_1}(f_\theta - f_{\theta'})(X_i), \dots, \frac{K}{n}\sum_{i \in I_K}(f_\theta - f_{\theta'})(X_i))$.

The authors demonstrated that the minmax MoM
estimators can weed out outliers during the train-
ing process when the contamination ratio is small.
Since the median of means is operated on loss val-
ues {($f_\theta - f_{\theta'}$)($X_i$)}$^{n}_{i=1}$, we refer to the minmax MoM
estimators as LMoM in the reminder of the paper.
Although the corresponding concentration bounds de-
vised in Lecué and Lerasle (2020) are statistically pow-
erful, they rely on computing the Rademacher com-
plexity, which limits their practical applicability. Most
importantly, the objective (4) is highly non-convex and
often exhibits numerous spurious local minima due to
the nature of the MoM functional.

Another line of research, aiming at heavy-tailed
data/noise, proposed to apply the MoM principle
at the gradient level, which amounts to perform-
ing a robust mean estimation of random vectors,
namely {$\nabla_{\theta} f_{\theta}(X_i)$}$_{i=1}^n$. In this framework, standard
optimization algorithms are modified by using the
median of mean-gradients instead of mean gradi-
ents. Given a sequence of step sizes {$\rho_l$)$_l$, the iterates
are defined for all $l \in \mathbb{N}_+$ as

$$
\theta^{(l+1)} = \theta^{(l)} - \rho_{l+1} \text{ median} (G_{\pi^{(l)},1}, \dots, G_{\pi^{(l)},K}), (5)
$$

where the batch-mean gradients are given by

$$
G_{\pi^{(l)}, k} = \frac{K}{n} \sum_{i \in I_{\pi^{(l)}, k}} \nabla_{\theta} f_{\theta^{(l)}}(X_i). \quad (6)
$$

Here, {$\pi^{(l)}$}$_{l \in \mathbb{N}_+}$ denotes a sequence of independent
random permutations, each drawn uniformly at ran-
dom. These permutations are used to reshuffle the

dataset at every iteration and to define the index
sets $I_{\pi^{(l)},k} = \pi^{(l)}(\{(k-1)\frac{n}{K} + 1, \dots, k\frac{n}{K}\})$.

In contrast to the univariate case, there are mul-
tiple ways to define a median in high-dimensional
spaces that give different statistical properties and ro-
bustness. Minsker (2015) gives an exponential con-
centration bound for the geometric median of mean
(GMoM). Prasad et al. (2020) used GMoM to esti-
mate gradient, at iteration l the gradient estimation is
defined as the solution of

$$
\arg \min_{\mu \in \mathbb{R}^p} \frac{1}{K} \sum_{k=1}^{K} \| G_{\pi^{(l)}, k} - \mu \|_2. \quad (7)
$$

More recently, Puchkin et al. (2024) proposed the Smoothed Median-of-Means (SMoM) for smooth convex objectives. They proposed to use component-wise MoM and integrated it into clipped variants of stochastic gradient descent. Component-wise median ignores the geometric structure when the covariance matrix of gradients is not axis-aligned. In addition, the breakpoint of the component-wise median decreases drastically in high dimensions.

Overall, existing MoM-based gradient estimators, including SMoM, lack guarantees against adversarial contamination: in particular, no result establishes robustness when the fraction of outliers exceeds about $\epsilon = 0.05$. By contrast, our use of the Tukey median preserves geometric structure and yields robustness to much higher contamination levels, up to $\epsilon = 1/3$, as formalized in Theorem 1.

4 Tukey-Median of Mean-Gradients
for Langevin Dyamics

In this section, we detail our Tukey Median-of-Means Stochastic Gradient Langevin Dynamics (TMoM-SGLD) method designed for solving (3), without any access to the index of clean examples. First, we introduce and study the robustness of Tukey-Median of Mean-Gradients. Second, we integrate these robustness results into a non-convex optimization scheme.

4.1 Tukey-Median of Mean-Gradients

Among the multiple ways to define a median of a set
of vectors $\mathcal{G} = \{G_1, \dots, G_K\} \subset \mathbb{R}^p$, the Tukey median
is the point in $\mathbb{R}^p$ of maximal Tukey depth.

Definition 1 (Tukey median (Tukey, 1975)). For any point $\eta \in \mathbb{R}^p$, the Tukey depth at $\eta$ denoted depth($\mathcal{G}$, $\eta$) is the minimum fraction of points lying on either side of any hyperplane passing through $\eta$:

$$
\text{depth}(\mathcal{G}, \eta) = \inf_{\|v\|_2=1} \frac{|\{G \in \mathcal{G} : \langle G - \eta, v \rangle \ge 0\}|}{K}. \quad (8)
$$
The Tukey median is then defined as the point that maximizes this depth:

$$ \text{Tukey}(\mathcal{G}) = \underset{\eta \in \mathbb{R}^p}{\operatorname{argmax}} \operatorname{depth}(\mathcal{G}, \eta). \quad (9) $$

One of our key contributions is to show that the Tukey median, when applied to batch-mean gradients (the Tukey-MoM gradient), provides a robust estimator of the uncorrupted mean gradient. In other words, we prove that the Tukey median of $\mathcal{G} = \{G_{\pi^{(l)},1}, \dots, G_{\pi^{(l)},K}\} = \mathcal{G}_{uncorr} \cup \mathcal{G}_{corr}$ (as in (6)), formed by both uncorrupted and corrupted batch-mean gradients, remains close to the true uncorrupted mean gradient even under high level of contamination. To establish this, we rely on the following assumptions.

(A1) *Data distribution.* The samples $\{X_i\}_{i \in \mathcal{I}}$ are i.i.d.

(A2) *Contamination bound.* The fraction of corrupted gradients $|\mathcal{G}_{corr}|/K$ is at most $\varepsilon < 1/3$.

Let $\mu = \mathbb{E}[G]$ and $\Sigma = \text{Cov}(G) = BB^\top$ denote the mean and covariance of uncorrupted gradients. We let $\sigma_1$ and $\sigma_{\text{rank}(\Sigma)}$ denote its largest and smallest singular values. The following assumptions apply to the distribution of these uncorrupted vectors $G \in \mathcal{G}_{uncorr}$.

(A3) *Concentration condition.* Along any direction, the distribution of the uncorrupted gradients is sub Gaussian. Formally, for all $u \ge 0$ and all unit vectors $v \in \mathbb{R}^p$, the distribution of $G \in \mathcal{G}_{uncorr}$ satisfies:

$$ \begin{aligned} & \mathbb{P}(\langle B^\dagger(G - \mu), v \rangle > \sqrt{\frac{K}{n}}u) \quad (10) \\ & \le \Phi^c(u) + \frac{c_1RK}{\sqrt{n}}\exp\left(-c_2\frac{u^2}{R^2}\right), \end{aligned} $$

for some appropriate constants $c_1$, $c_2$ and $R$ (Austern and Mackey, 2022). Here $\Phi^c(u) := 1 - \Phi(u)$ denotes the tail of the standard normal distribution. For convenience, we introduce $\epsilon_{res} > 0$ as a residual term so that the inequality may be written as $\mathbb{P}(\langle B^\dagger(G - \mu), v \rangle > \sqrt{K/n}u) \le \Phi^c(u) + \epsilon_{res}$.

(A4) *Bounded anisotropy.* The maximal directional imbalance of the marginal distributions of $G$ is bounded, i.e., there exists a finite

$$ \epsilon_{\text{iso}} = \max_{\|v\|_2=1, \lambda \in \mathbb{R}} \mathbb{P}[\langle G - \mu, v \rangle \ge \lambda] - \min_{\|v\|_2=1, \lambda \in \mathbb{R}} \mathbb{P}[\langle G - \mu, v \rangle \ge \lambda]. \quad (11) $$

We are now well-equipped to state our guarantees regarding Tukey-Median of Mean-Gradients, inspired by Li (2019).

**Theorem 1** (Robustness of the Tukey-MoM gradient). Suppose assumptions (A1)–(A4) hold and let $t \ge 0$. Define the critical separation radius

$$ r_{\star} = \sqrt{\frac{K}{n}} \Phi^{-1}\left(\frac{1}{2} + 2\epsilon_{\text{iso}} + \epsilon_{\text{res}} + 2t + \frac{\varepsilon}{1-\varepsilon}\right). \quad (12) $$

Then, with probability at least $1 - 2e^{-Kt^2}$, provided that the batch-size satisfies

$$ n/K \ge C \log(1/\epsilon_{\text{res}})/r_{\star}^2 \quad (13) $$

for some constant $C$ depending only on $\{c_1, c_2, R\}$ from (A3), the Tukey-MoM gradient yields

$$ \|\text{Tukey}(\mathcal{G}) - \mu\|_2 \le \left[ \varepsilon_{\mu} := \frac{\sigma_1}{\sigma_{\text{rank}(\Sigma)}} r_{\star} \right], \quad (14) $$

where $\mu = \mathbb{E}[G]$ and $\Sigma = \text{Cov}(G) = BB^\top$ denote the mean and covariance of the uncorrupted gradients.

**Theorem 1** establishes a probabilistic non-asymptotic bound on the error of the Tukey median estimator when applied to $K$ batch-mean gradients formed from $n$ samples and subject to an $\varepsilon$-fraction of arbitrary outliers. The deviation is controlled in $\ell_2$-norm and bounded by the critical separation radius $r_\star$, whose leading term scales as $\sqrt{K/n}$. This rate is modulated by the spectral condition number of the covariance $\Sigma$, the anisotropy parameter $\epsilon_{iso}$, the residual tail mass $\epsilon_{res}$ stemming from mild concentration assumptions, and a contamination term. In essence, the estimator achieves the expected $\sqrt{K/n}$ accuracy up to condition-number and residual corrections.

**Remark 1.** The critical separation radius $r_\star$ characterizes the largest deviation from the true mean $\mu$ that the Tukey-MoM gradient can attain under $\varepsilon$-contamination. In particular, any candidate point farther than $r_\star$ from $\mu$ (in the covariance-normalized geometry) has insufficient Tukey depth to be selected, which guarantees robustness of the estimator.

For comparison, our closest competitors SMoM (Puchkin et al., 2024) and GMoM (Prasad et al., 2020) bound the error of uncorrupted mean gradient estimation by a term depending on Tr($\Sigma$) while ours depends on the condition number of $\Sigma$. Finally, while their analysis are restricted to convex objectives, our robustification serves as the basis for extending the study of median-of-means methods to non-convex regimes, as developed in Section 4.3.

## 4.2 Proof of Theorem 1

We show that the Tukey median lies within a covariance-normalized radius of the true mean uncorrupted
gradient $\mu$. The proof proceeds in three steps: (i) establish a lower bound on the depth at $\mu$; (ii) upper bound the depth at a candidate $\eta$; (iii) compare the two to obtain the critical separation radius $r_\star$.

Let $\Sigma = BB^\top$ denote the Cholesky factorisation of $\Sigma$. Let $B^\dagger := (B^\top B)^{-1}B^\top$ denote the Moore-Penrose pseudo-inverse of $B$.

**Step 1.** Lower bound on the depth at $\mu$. For any $v$ of unit norm in $\mathbb{R}^d$, we have

$$ \begin{aligned} & |\{G \in \mathcal{G} : \langle G - \mu, v \rangle \ge 0\}| \\ & = |\{G \in \mathcal{G}_{\text{uncorr}} : \langle G - \mu, v \rangle \ge 0\}| \\ & \quad + |\{G \in \mathcal{G}_{\text{corr}} : \langle G - \mu, v \rangle \ge 0\}|. \end{aligned} \quad (15) $$

Hence,

$$ \begin{aligned} & |\{G \in \mathcal{G} : \langle G - \mu, v \rangle \ge 0\}| \\ & \ge |\{G \in \mathcal{G}_{\text{uncorr}} : \langle G - \mu, v \rangle \ge 0\}|. \end{aligned} \quad (16) $$

For any $i \in \mathcal{I}_{\text{uncorr}}$, define the binary variable

$$ Y_i = \begin{cases} 1 & \text{if } \langle G_i - \mu, v \rangle \ge 0, \\ 0 & \text{otherwise.} \end{cases} \quad (17) $$

That is, $Y_i = 1$ if the $i$-th uncorrupted gradient lies in the halfspace $\{g : \langle g - \mu, v \rangle \ge 0\}$, and $Y_i = 0$ otherwise. By (A4), the imbalance across halfspaces is bounded, hence $\mathbb{P}(Y_i = 1) \ge \frac{1}{2} - \epsilon_{\text{iso}}$. By construction, the number of uncorrupted gradients in this halfspace is given by

$$ |\{G \in \mathcal{G}_{\text{uncorr}} : \langle G - \mu, v \rangle \ge 0\}| = \sum_{i \in \mathcal{I}_{\text{uncorr}}} Y_i. \quad (18) $$

Applying Hoeffding's inequality yields, with probability at least $1 - e^{-(1-\epsilon)Kt^2}$,

$$ \frac{1}{(1 - \epsilon)K} |\{G \in \mathcal{G} : \langle G - \mu, v \rangle \ge 0\}| \ge \frac{1}{2} - \epsilon_{\text{iso}} - t. \quad (19) $$

Hence, we deduce that the depth at $\mu$ is, with high probability, larger than $1/2 - \epsilon_{\text{iso}} - t$, from which we deduce that the Tukey median depth is itself larger than or equal to $1/2 - \epsilon_{\text{iso}} - t$.

**Step 2.** Upper bound on the depth at a candidate $\eta$. Let us consider any candidate point $\eta \in \mathbb{R}^d$. If $\|\eta - \mu\|_2$ is large, we show that its depth is too small for it to be a Tukey median. For this purpose, suppose that the unit norm vector $v \in \mathbb{R}^d$ is optimal for $\eta$. Arguing as above,

$$ |\{G \in \mathcal{G} : \langle G - \eta, v \rangle \ge 0\}| \le |\{G \in \mathcal{G}_{\text{uncorr}} : \langle G - \mu, v \rangle \ge \langle \eta - \mu, v \rangle\}| + \varepsilon K. $$

Equivalently

$$ \begin{aligned} & |\{G \in \mathcal{G} : \langle G - \eta, v \rangle \ge 0\}| \\ & \le \left| \left\{ G \in \mathcal{G}_{\text{uncorr}} : \left\langle B^{\dagger}(G - \mu), \frac{\eta - \mu}{\|B^{\top}(\eta - \mu)\|_2} \right\rangle \ge r(\eta, \mu) \right\} \right| + \varepsilon K. \end{aligned} \quad (20) $$

where $r(\eta, \mu) = \|B^{\dagger/2}(\eta - \mu)\|_2^2 / \|B^{\dagger}(\eta - \mu)\|_2$ denotes its covariance-normalized separation from the mean $\mu$, hereafter denoted $r$.

Applying Assumption (A3) and Hoeffding's inequality gives, with probability at least $1 - e^{-Kt^2}$,

$$ \frac{1}{(1-\epsilon)K} |\{G \in \mathcal{G}_{\text{uncorr}} : \langle B^{\dagger}(G - \mu), \frac{\eta - \mu}{\|B^{\top}(\eta - \mu)\|_2} \rangle \ge r\}| \le \Phi^c(\sqrt{n/K} r) + \frac{c_1 R \sqrt{K}}{\sqrt{n}} \exp \left(-c_2 \frac{n r^2}{K R^2}\right) + \epsilon_{\text{iso}} + t. $$

Including the contribution of corrupted gradients then gives

$$ \begin{aligned} & \frac{1}{(1 - \epsilon)K} |\{G \in \mathcal{G} : \langle G - \mu, v \rangle \ge 0\}| \\ & \le \Phi^c(\sqrt{n/K} r) + \frac{c_1 R \sqrt{K}}{\sqrt{n}} \exp \left(-c_2 \frac{n r^2}{K R^2}\right) + \epsilon_{\text{iso}} + t \\ & \quad + \epsilon_{\text{iso}} + t + \frac{\epsilon}{1 - \epsilon}. \end{aligned} \quad (21) $$

This shows that, as the separation $r$ from $\mu$ grows, the fraction of gradients supporting $\eta$ in direction $v$ decays rapidly, limiting the possible depth of distant candidates.

**Step 3.** Critical separation radius and conclusion. We now examine what could make the upper bound (21) on $|\{G \in \mathcal{G} : \langle G - \eta, v \rangle \ge 0\}|$ become lower than the lower bound (19) on $|\{G \in \mathcal{G} : \langle G - \mu, v \rangle \ge 0\}|$, i.e.

$$ \begin{aligned} & \Phi^c(\sqrt{n/K} r) + \frac{c_1 R \sqrt{K}}{\sqrt{n}} \exp \left(-c_2 \frac{n r^2}{K R^2}\right) + \epsilon_{\text{iso}} + t \\ & \quad + \frac{\epsilon}{(1 - \epsilon)} \le \frac{1}{2} - \epsilon_{\text{iso}} - t \end{aligned} \quad (22) $$

This inequality is implied by

$$ \begin{aligned} \Phi^c(\sqrt{n/K} r) \le \frac{1}{2} - 2\epsilon_{\text{iso}} - \frac{c_1 R \sqrt{K}}{\sqrt{n}} \exp \left(-c_2 \frac{n r^2}{K R^2}\right) \\ - 2t - \frac{\epsilon}{(1 - \epsilon)}. \end{aligned} \quad (23) $$

By symmetry of the Gaussian distribution, this is equivalent to

$$ -\sqrt{\frac{n}{K}} r \le \Phi^{-1} \left( \frac{1}{2} - 2\epsilon_{\text{iso}} - \epsilon_{\text{res}} - 2t - \frac{\epsilon}{(1 - \epsilon)} \right), \quad (24) $$
as long as the exponential remainder is bounded $\epsilon_{res}$, i.e., $\frac{c_1 R \sqrt{K}}{\sqrt{n}} \exp(-c_2 \frac{n}{K} r^2) \le \epsilon_{res}$. Equivalently,

$$r \ge r_{\star} := \sqrt{\frac{K}{n}} \Phi^{-1}(\frac{1}{2} + 2\epsilon_{iso} + \epsilon_{res} + 2t + \frac{\epsilon}{1-\epsilon}). \quad (25)$$

Since $r \ge \frac{\sigma_{\text{rank}(\Sigma)}}{\sigma_1} \|\eta - \mu\|_2$, we get

$$\|\eta - \mu\|_2 \ge \sqrt{\frac{K}{n}} \frac{\sigma_1}{\sigma_{\text{rank}(\Sigma)}} r_{\star} \quad (26)$$

implies (25), and thus implies that $\eta$ cannot be a Tukey median. Finally, (25) implies that (24) is satisfied for as long as $n/K \ge C \log(1/\epsilon_{res})/r^2$, for some constant depending on $c_1$, $c_2$, and $R$. This completes the proof.

## 4.3 Hitting Time Guarantees for TMoM-SGLD

We next present our result on TMoM-SGLD, which integrates Tukey-MoM gradients into the Stochastic Gradient Langevin Dynamics (SGLD) framework.

SGLD is particularly appealing as it augments stochastic gradient descent with Gaussian noise, allowing efficient exploration of non-convex loss landscapes and ensuring polynomial-time convergence to neighborhoods of stationary points (Zhang et al., 2017). In TMoM-SGLD, we retain these favorable properties while enhancing robustness by replacing the empirical mean gradient with the Tukey-MoM gradient, which mitigates the impact of corrupted samples. Our algorithm is defined below.

**Definition 2 (TMoM-SGLD).** Given a sequence of step sizes $\{\rho_l\}_{l \in \mathbb{N}}$, the iterates are defined for all $l \in \mathbb{N}_+$ as

$$\theta^{(l+1)} = \theta^{(l)} - \rho_{l+1} \text{ Tukey } [G_{\pi^{(l)},1}, \dots, G_{\pi^{(l)},K}] + \sqrt{2\rho_{l+1}}\zeta^{(l+1)}, \quad (27)$$

where $\zeta^{(l+1)} \sim \mathcal{N}(0, I)$ is a standard Gaussian noise.

Our result below establishes non-asymptotic probabilistic convergence guarantees: with high probability, we bound the number of iterations required for TMoM-SGLD to reach a region where the expected gradient of the uncorrupted loss is small.

**Theorem 2 (Hitting time of TMoM-SGLD iterates).** Let $\epsilon > 0$ be a desired accuracy and consider the iterates $\{\theta^{(l)}\}_{l \in \mathbb{N}^+}$ generated by TMoM-SGLD with step-sizes of the form $\rho_l = \rho_0 l^{-1/2}$. Define the hitting time

$$\tau_{\epsilon} = \inf \{ l \ge 0 : \| \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)}) \| \le \epsilon \}.$$

Suppose assumptions (A1)–(A4) hold and that $n/K \ge C \log(1/\epsilon_{res})/r_{\star}^2$, as required in Theorem 1. Then, for

all $L \ge 1$, the hitting time of the TMoM-SGLD iterates satisfies

$$\mathbb{P} (\tau_{\epsilon+\epsilon_{\mu}} \ge L) \le 2 \frac{\|\nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(0)})\|_2}{\rho_0 (\epsilon + \epsilon_{\mu})^2 \sqrt{L}} + 6\rho_0 \frac{\log(L) (\epsilon_{\mu} + B_1)}{(\epsilon + \epsilon_{\mu})^2 \sqrt{L}} + 6 \frac{\delta_0^2 p}{(\epsilon + \epsilon_{\mu})^2} + 2Le^{-Kt^2}. \quad (28)$$

Our result is an extension of (Chen et al., 2020, Theorem 2) where a different treatment of the main results from Zhang et al. (2017) is presented as well as interesting generalizations. The proof follows the lines of the proof given in Chen et al. (2020) and the details are given in the supplementary material.

Theorem 2 provides a non-asymptotic high-probability upper bound on the hitting time $\tau_{\epsilon+\epsilon_{\mu}}$ of a small-gradient region for the iterates of TMoM-SGLD. Specifically, it bounds the number of iterations $L$ needed before the expected robustified gradient norm—estimated via a Tukey-MoM procedure—falls below $\epsilon+\epsilon_{\mu}$, where $\epsilon_{\mu}$ is the statistical error from Theorem 1. The bound depends on the initial energy gap, the decaying step-size $\rho_l = \rho_0 l^{-1/2}$, $\epsilon_{\mu}$, and problem parameters $(d, \delta_0, B_1)$ such as local smoothness constants. Unlike prior works such as SMoM (Puchkin et al., 2024), which do not track iteration-wise dynamics, our result explicitly links statistical robustness and optimization dynamics, showing how the combination governs the time to reach near-stationary points under contamination.

Moreover, this analysis captures the explicit trade-off between robustness and efficiency: more conservative (i.e., robust) estimators may lead to larger $\epsilon_{\mu}$, which slows convergence, but protect against catastrophic gradient corruption. Compared to classical convergence guarantees for SGLD under smooth and bounded gradient noise (Zhang et al., 2017; Raginsky et al., 2017), our result accounts for structured corrections and gradient bias introduced by robustification. This is particularly relevant in high-dimensional or federated settings, where malicious updates or heavy perturbations are common and traditional moment assumptions fail to hold (Diakonikolas et al., 2020).

We conclude this section by mentioning that when training Deep Neural Networks, the Polyak-Lojasiewicz (PL) condition justifies that we only consider reaching region of the weight space where the gradient is small. Notice further that the landscape around global minimisers in Deep Learning may not be convex, which justifies not imposing the condition that the Hessian be positive semi-definite in the study of the hitting time.

**Remark 2 (Approximate Tukey Median).** Computing
<figure><img src="image_4.png" /><figcaption>Figure 1: Robustness of different methods under increasing contamination in linear regression on synthetic data. (Left) Test error (log-scale) as a function of contamination level $\epsilon$. (Right) Average number of times contaminated samples are used during training</figcaption></figure>

the Tukey depth is notoriously difficult Bremner and Chen (2009) and known methods are of exponential complexity in the dimension of the data. In practice, will compute an approximate version of the Tukey median to reduce computational complexity. Instead of searching over all $\eta \in \mathbb{R}^p$, we restrict the candidate set to $G = \{G_1, \dots, G_K\}$, and select the mean gradient with the highest Tukey depth. To train neural networks, we accelerate the computations by first applying principal-component analysis (PCA) to the matrix of mean gradients. In practice, we decided to keep, at least 67% of the total variance, thereby preserving the essential information while substantially reducing the computational load. A thorough theoretical analysis of efficient approximations of the Tukey depth was recently proposed in Briënd et al. (2023) and there is a wide avenue for further practical and theoretical improvements of the computation of the Tukey depth and the associated Tukey median.

# 5 Numerical Experiments

We now turn to the numerical study of our approach. We first consider a convex setting to assess the robustness of Tukey MoM (TMoM) in the presence of heavy-tailed noise and adversarial outliers. Then we move to a more realistic non-convex scenario.

## 5.1 Baselines and Optimizers

We compare our method against several baselines: standard optimizers (SGD and Adam); a robust Adam variant for heavy-tailed noise and outliers, TAdam (Ilboudo et al., 2022); MoM baselines at the gradient level, GMoM (Prasad et al., 2020) and SMoM (Puchkin et al., 2024); and a MoM baseline at the loss level, LMoM (Lecué and Lerasle, 2020).

For all baselines, we resort to SGD optimizer for the

convex case (Section 5.2) and Adam for the non-convex case (Section 5.3), since Adam is widely acknowledged to outperform SGD/SGLD in deep learning tasks through its adaptive learning-rate scaling, which reduces the need for manual tuning.

## 5.2 Convex Case: Linear Regression

We begin with a simulated linear regression task, a standard benchmark in the robust learning literature, to validate our method in a convex setting.

**Dataset.** We generate $n = 2000$ uncorrupted samples $X_i = (Z_i, Y_i)_{i=1}^n$ from the linear model

$$Y = Z\theta + \xi, \quad \theta = [1, 1, 1]^T, \quad \xi \sim \mathcal{N}(0, 0.1^2), \quad (29)$$

where $Z \in \mathbb{R}^3$. Unlike prior works (Lecué and Lerasle, 2020; Prasad et al., 2020), which assume $Z \sim \mathcal{N}(0, I)$, we sample $Z$ from a Gaussian distribution with covariance matrix $\Sigma$ whose eigenvalues are $(100, 10, 1)$ and whose principal axes are chosen at random.

To introduce contamination, following Lecué and Lerasle (2020), an $\epsilon$-fraction of training labels is replaced by a fixed outlier value $Y = 1000$, with $\epsilon \in \{0, 0.02, 0.04, 0.06, 0.08, 0.10\}$. The dataset is then split evenly into training and test sets.

This design serves two purposes. First, it provides a setting that is closer to heavy-tailed assumptions, moving beyond the standard isotropic benchmarks. Second, Assumption (A4) in Theorem 1 explicitly allows for anisotropic gradient distributions. In fact, we show in the appendix that under the linear regression model (29), choosing $Z \sim \mathcal{N}(0, I)$ results in a quasi-identity covariance for the gradients $G$, whereas an anisotropic $\Sigma$ produces a fully anisotropic gradient covariance, thereby offering a more challenging and representative testbed.

**Optimization.** All methods are tested with learning
rates $\gamma \in \{10^{-2}, 1/3 \times 10^{-2}, 10^{-3}, 1/3 \times 10^{-3}, 10^{-4}\}$ and trained for 2000 epochs. For MoM-based methods, we set $K = 25$, corresponding to a batch size of 40 = 1000/25. At each epoch, the training set is reshuffled. We additionally track the indices of corrupted samples, following Lecué and Lerasle (2020), to count how often contaminated samples contribute to updates in our method and in LMoM.

Results are averaged over five independent runs for all contamination levels, each corresponding to a new train/test split, parameter initialization, covariance matrix generation, and random permutation sequence. We report the best mean test Mean Squared Error (MSE) across learning rates, along with the corresponding standard deviations of the means.

**Results.** Figure 1 (left) reports the log-scale test errors as a function of the contamination level $\epsilon$. As expected, standard optimizers (SGD, Adam) and even the robust variant TAdam quickly diverge as soon as corrupted samples are introduced (i.e., $\epsilon > 0$), showing a steep increase in error. Among the MoM baselines, LMoM, GMoM, and SMoM offer some robustness but still degrade significantly as $\epsilon$ grows. In contrast, TMoM maintains the same performance as in the clean setting, remaining stable up to $\epsilon = 0.06$. This highlights TMoM's stronger robustness against outliers compared to the baselines.

**Mitigating corrupted gradients.** Figure 1 (right) shows that TMoM effectively weeds out corrupted samples during training. For standard optimizers such as SGD and Adam, gradients from all samples are always used. Similarly, GMoM and SMoM incorporate all samples unless the gradient aggregation itself is robust. In contrast, LMoM and TMoM select batches for optimization, enabling them to discard corrupted batches whenever the contamination ratio is not too high and $K$ is not too small. To assess this ability, we plot the total number of times corrupted samples are used during training, averaged across five runs for all contamination levels. TMoM completely avoids corrupted gradients at low contamination levels. Even when contamination is unavoidable at higher $\epsilon$, TMoM consistently selects less corrupted batches than LMoM. This robustness directly stems from the high breakdown point of the Tukey median, which underpins the superior resilience of TMoM to contamination.

## 5.3 Non-Convex Case: Image Classification

We now turn to the non-convex setting. We consider a deep learning task on a real image classification dataset, providing the first empirical study of MoM-based robustification in non-convex optimization.

**Dataset.** We consider a binary classification task derived

<figure><img src="image_5.png" /><figcaption>Figure 2: Robustness of MoM estimators on real-data.</figcaption></figure>

from CIFAR-10 by selecting only the *cat* (class 3) and *dog* (class 5) categories. To simulate corrupted samples, we introduce contamination levels $\epsilon \in \{0.05, 0.10, 0.15, 0.20, 0.25\}$ by flipping the corresponding proportion of labels in the training set. The test set remains clean. We use a convolutional neural network, whose implementation details can be found in the supplementary material.

**Results.** As standard methods and GMoM diverge, we report the final error rate (1 - accuracy) of all other MoM-based methods. Results are reported in Figure 2 and show that all MoM variants somewhat perform comparably in terms of robustness and accuracy. These early findings suggest that MoM-based strategies may extend to non-convex settings, but more extensive experiments are required to draw definitive conclusions.

# 6 Conclusion and Perspectives

We introduced a Median-of-Means approach for robust gradient estimation based on the Tukey median, and showed that it provides a key ingredient for robustifying empirical risk minimization in the presence of outliers. Building on this estimator, we analyzed a robustified version of Stochastic Gradient Langevin Dynamics and established non-asymptotic hitting-time guarantees. Numerical experiments further demonstrated that our method outperforms baselines when the portion of outliers increases.

Future research directions include addressing the computational challenges of approximating the Tukey depth in high dimensions, which may affect both scalability and accuracy, and developing adaptive strategies, such as bandit-based procedures, for selecting minibatch sizes, along with a theoretical study of their performance. Another promising direction is to extend our robustification beyond gradients to momentum terms, which would be particularly relevant for adaptive optimizers such as Adam in non-convex deep learning.
## References

Audibert, J.-Y., Catoni, O., et al. (2011). Robust linear least squares regression. *The Annals of Statistics*, 39(5):2766-2794.

Austern, M. and Mackey, L. (2022). Efficient concentration with gaussian approximation. *arXiv preprint arXiv:2208.09922*.

Baraud, Y., Birgé, L., and Sart, M. (2017). A new method for estimation and model selection: *p*-estimation. *Inventiones mathematicae*, 207(2):425-517.

Bremner, D. and Chen, D. (2009). A branch and cut algorithm for the halfspace depth problem. *arXiv preprint arXiv:0910.1923*.

Briend, S., Lugosi, G., and Oliveira, R. I. (2023). On the quality of randomized approximations of tukey's depth. *arXiv preprint arXiv:2309.05657*.

Catoni, O. (2012). Challenging the empirical mean and empirical variance: A deviation study. *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques*, 48(4):1148-1185.

Chen, X., Du, S. S., and Tong, X. T. (2020). On stationary-point hitting time and ergodicity of stochastic gradient Langevin dynamics. *J. Mach. Learn. Res.*, 21(1).

Chinot, G., Lecué, G., and Lerasle, M. (2019). Robust statistical learning with lipschitz and convex loss functions. *Probability Theory and related fields*, pages 1-44.

Diakonikolas, I., Kane, D. M., and Stewart, A. (2020). Outlier-robust high-dimensional gradient descent for stochastic optimization. In *Advances in Neural Information Processing Systems (NeurIPS)*, volume 33, pages 4399-4411.

Donoho, D. L. and Gasko, M. (1992). Breakdown properties of location estimates based on halfspace depth and projected outlyingness. *The Annals of Statistics*, pages 1803-1827.

Dyckerhoff, R. and Mozharovskyi, P. (2016). Exact computation of the halfspace depth. *Computational Statistics & Data Analysis*, 98:19-30.

Huang, S.-T. and Lederer, J. (2023). Deepmom: Robust deep learning with median-of-means. *Journal of Computational and Graphical Statistics*, 32(1):181-195.

Ilboudo, W. E. L., Kobayashi, T., and Sugimoto, K. (2022). Robust stochastic gradient descent with student-t distribution based first-order momentum. *IEEE Transactions on Neural Networks and Learning Systems*, 33(3):1324-1337.

Jiang, L., Zhou, Z., Leung, T., Li, L.-J., and Fei-Fei, L. (2018). Mentornet: Learning data-driven curriculum for very deep neural networks on corrupted labels. In *International conference on machine learning*, pages 2304-2313. PMLR.

Kwon, J., Lecué, G., and Lerasle, M. (2018). A mom-based ensemble method for robustness, subsampling and hyperparameter tuning. *arXiv preprint arXiv:1812.02435*.

Lecué, G. and Lerasle, M. (2019). Learning from machine's principles: Le cam's approach. *Stochastic Processes and their applications*, 129(11):4385-4410.

Lecué, G. and Lerasle, M. (2020). Robust machine learning by median-of-means: Theory and practice. *The Annals of Statistics*, 48(2):906-931.

Lerasle, M. (2019). Lecture notes: Selected topics on robust statistical learning theory. *arXiv preprint arXiv:1908.10761*.

Li, J. (2019). Mean estimation in high dimensions.

Lugosi, G. and Mendelson, S. (2016). Risk minimization by median-of-means tournaments. *arXiv preprint arXiv:1608.00757*.

Merad, I. and Gaiffas, S. (2024). Robust stochastic optimization via gradient quantile clipping. *Transactions on Machine Learning Research*.

Minsker, S. (2015). Geometric median and robust estimation in Banach spaces. *Bernoulli*, 21(4):2308-2335.

Minsker, S. (2018). Uniform bounds for robust mean estimators. *arXiv preprint arXiv:1812.03523*.

Nemirovsky, A. S. and Yudin, D. B. (1983). *Problem Complexity and Method Efficiency in Optimization*. Wiley-Interscience Series in Discrete Mathematics. John Wiley & Sons, Inc., New York. Translated from the Russian and with a preface by E. R. Dawson.

Prasad, A., Suggala, A. S., Balakrishnan, S., and Ravikumar, P. (2020). Robust estimation via robust gradient estimation. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 82(3):601-627.

Puchkin, N., Gorbunov, E., Kutuzov, N., and Gasnikov, A. (2024). Breaking the heavy-tailed noise barrier in stochastic optimization problems. In *International Conference on Artificial Intelligence and Statistics*, pages 856-864. PMLR.

Raginsky, M., Rakhlin, A., and Telgarsky, M. (2017). Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis. In *Proceedings of the 30th Conference on Learning Theory (COLT)*, pages 1674-1703. PMLR.
Ren, M., Zeng, W., Yang, B., and Urtasun, R. (2018). Learning to reweight examples for robust deep learning. In *International conference on machine learning*, pages 4334–4343. PMLR.

Tukey, J. W. (1975). Mathematics and the picturing of data. In *Proceedings of the international congress of mathematicians*, volume 2, pages 523–531. Vancouver.

Zhang, Y., Liang, P., and Charikar, M. (2017). A hitting time analysis of stochastic gradient Langevin dynamics. In *Conference on Learning Theory*, pages 1980–2022. PMLR.
Supplementary Materials

In this appendix, we first present the proof of Theorem 2.<fnref n="1" /> We then show that in our experimental design of section 5.2, the anisotropy of data covariance affects the anisotropy of gradient covariance. Finally, we include the implementation details of our real experiments.

A PROOF OF THEOREM 2

Recall that, given a sequence of step sizes {$ρ_l$}$_{l∈N}$, the iterates are defined for all $l ∈ N_+$ as

$$ \theta^{(l+1)} = \theta^{(l)} - \rho_{l+1} \text{ Tukey } [G_{\pi^{(l)},1}, \dots, G_{\pi^{(l)},K}] + \sqrt{2\rho_{l+1}} \zeta^{(l+1)}, \quad (30) $$

where $\zeta^{(l+1)} \sim \mathcal{N}(0, I)$ is a standard Gaussian noise. Consider the sequence of gradient sets {$G^{(l)}$}$_{l=1}^L$, where each $G^{(l)} = \{G_{\pi^{(l)},1}, \dots, G_{\pi^{(l)},K}\}$ is computed from the iterates in (30). Let $\mathcal{E}_L$ denote the event

$$ \|\text{Tukey}(G^{(l)}) - \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2 \le \left[ \varepsilon_{\mu} := \frac{\sigma_1}{\sigma_{\text{rank}(\Sigma)}} r_{\star} \right], $$

for $l = 1, \dots, L$, with

$$ r_{\star} = \sqrt{\frac{K}{n}} \Phi^{-1}\left(\frac{1}{2} + 2\epsilon_{\text{iso}} + \epsilon_{\text{res}} + 2t + \frac{\epsilon}{1-\epsilon}\right). $$

as in defined in Theorem 1. Recall from Theorem 1, that the probability of $\mathcal{E}_L^c$ at most $2Le^{-Kt^2}$.

Let $h^{(l)} := \theta_{l+1} - \theta_l$ denote the step at iteration $l$. Then, on the event $\mathcal{E}_L$ and by the construction of the algorithm,

$$ h^{(l)} = -\eta_{l+1} \left( \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)}) + \epsilon_{\mu}^{(l+1)} + \xi^{(l+1)} \right) + \sqrt{\eta_{l+1}} \zeta^{(l+1)} \quad (31) $$

Notice further that,

$$ \hat{F}_{\text{uncorr}}(\theta^{(l+1)}) = \hat{F}_{\text{uncorr}}(\theta^{(l)}) + \langle \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)}), h^{(l)} \rangle + \frac{1}{2} \int_{0}^{1} \int_{0}^{s} h^{(l)t} \nabla_{\theta}^{2} \hat{F}_{\text{uncorr}}(\theta^{(l)} + r h^{(l)}) h^{(l)} dr ds, $$

and combining this last equation with (31), we get

$$ \hat{F}_{\text{uncorr}}(\theta^{(l+1)}) = \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)}) + \langle \nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l)}), -\eta_{l+1} (\nabla_{\theta} \hat{F}_{\text{uncorr}}(\theta^{(l+1)}) + \epsilon_{\mu}^{(l+1)} + \xi^{(l+1)}) + \sqrt{\eta_{l+1}} \zeta^{(l+1)} \rangle \\ + \frac{1}{2} h^{(l)t} \nabla_{\theta}^{2} \hat{F}_{\text{uncorr}}(\theta^{(l)} + r_{l+1} h^{(l)}) h^{(l)} \quad (32) $$

for some $r_{l+1} \in (0, 1)$. Notice that, by the definition of $C_2$, the last term satisfies

$$ \frac{1}{2} h^{(l)t} \nabla^2 \hat{F}_{\text{uncorr}}(\theta^{(l)} + r_{l+1} h^{(l)}) h^{(l)} \le \frac{C_2}{2} \|h^{(l)}\|_2^2 $$

and using $\|a + b + c\|_2^2 \le 3 (\|a\|_2^2 + \|b\|_2^2 + \|c\|_2^2)$, we get

$$ \|h^{(l)}\|_2^2 \le 3\eta_{l+1}^2 \|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l+1)})\|_2^2 + 3\eta_{l+1}^2 \|\epsilon_{\mu}^{(l+1)}\|_2^2 + 3\eta_{l+1}^2 \|\xi^{(l+1)}\|_2^2 + 3\eta_{l+1} \|\zeta^{(l+1)}\|_2^2. $$

Taking the expectation, conditional on $\mathcal{E}_L$ in this last inequality, we get

$$ \mathbb{E}[\|h^{(l)}\|_2^2 | \mathcal{E}_L] \le 3\eta_{l+1}^2 \mathbb{E}[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l+1)})\|_2^2 | \mathcal{E}_L] + 3\eta_{l+1}^2 \|\epsilon_{\mu}^{(l+1)}\|_2^2 + 3\eta_{l+1}^2 B_1 + 3\eta_{l+1} \delta_0^2 d. $$

<footnote n="1">¹A similar analysis for Adam is a challenging direction that we plan to investigate in future work.</footnote>
Plugging this into (32) after taking the expectation, we get

$$
\begin{aligned}
\mathbb{E}[\nabla \hat{F}_{\text{uncorr}}(\theta^{(l+1)}) | \mathcal{E}_L] &\le \mathbb{E}[\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) | \mathcal{E}_L] \\
&\quad - \eta_{l+1} \mathbb{E}[\langle \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}), \epsilon_{\mu}^{(l+1)} \rangle | \mathcal{E}_L] - \eta_{l+1} \mathbb{E}[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2^2 | \mathcal{E}_L] \\
&\quad + \frac{3C_2}{2} \eta_{l+1}^2 \mathbb{E}[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l+1)})\|_2^2 | \mathcal{E}_L] + 3\eta_{l+1}^2 \|\epsilon_{\mu}^{(l+1)}\|_2^2 + 3\eta_{l+1}^2 B_1 + 3\eta_{l+1}\delta_0^2 d.
\end{aligned}
$$

which gives

$$
\begin{aligned}
\mathbb{E}[\nabla \hat{F}_{\text{uncorr}}(\theta^{(l+1)}) | \mathcal{E}_L] &\le \mathbb{E}_n[\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) | \mathcal{E}_L] - \eta_{l+1} \langle \mathbb{E}[\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) | \mathcal{E}_L], \epsilon_{\mu}^{(l+1)} \rangle \\
&\quad - (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2) \mathbb{E}[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2^2 | \mathcal{E}_L] + 3\eta_{l+1}^2 \|\mu^{(l+1)}\|_2^2 + 3\eta_{l+1}^2 B_1 + 3\eta_{l+1}\delta_0^2 d.
\end{aligned}
$$

On the other hand,

$$
\begin{aligned}
\eta_{l+1} \mathbb{E}[\langle \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}), \mu^{(l+1)} \rangle | \mathcal{E}_L] &+ (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2) \mathbb{E}[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2^2 | \mathcal{E}_L] \\
&= (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2) \mathbb{E}\left[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2^2 + \frac{\eta_{l+1}}{\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2} \epsilon_{\mu}^{(l+1)}\|_2^2 | \mathcal{E}_L\right]
\end{aligned}
$$

Therefore, we get

$$
\begin{aligned}
\mathbb{E}[\hat{F}_{\text{uncorr}}(\theta^{(l+1)}) | \mathcal{E}_L] &\le \mathbb{E}[\hat{F}_{\text{uncorr}}(\theta^{(l)}) | \mathcal{E}_L] - (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2) \mathbb{E}\left[\|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\right. \\
&\qquad \left. + \frac{\eta_{l+1}}{\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2} \epsilon_{\mu}^{(l+1)}\right] \|_2^2 | \mathcal{E}_L \right] \\
&\quad + 3\eta_{l+1}^2 \|\epsilon_{\mu}^{(l+1)}\|_2^2 + 3\eta_{l+1}^2 B_1 + 3\eta_{l+1}\delta_0^2 d.
\end{aligned}
$$

Rearranging the terms, and summing over $l = 0, \dots, L-1$, we obtain

$$
\begin{aligned}
\sum_{l=0}^{L-1} (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2) \mathbb{E} \left[ \left\| \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) + \frac{\eta_{l+1}}{\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2} \epsilon_{\mu}^{(l+1)} \right\|_2^2 \middle| \mathcal{E}_L \right] \\
\le \mathbb{E} \left[ \hat{F}_{\text{uncorr}}(\theta^{(0)}) \middle| \mathcal{E}_L \right] - \mathbb{E} \left[ \hat{F}_{\text{uncorr}}(\theta^{(L)}) \middle| \mathcal{E}_L \right] + 3\eta_{1:L}^2 \epsilon_{\mu} + 3\eta_{1:L}^2 B_1 + 3\eta_{l+1}\delta_0^2 d.
\end{aligned}
$$

Define

$$
\tilde{\eta}_{l+1} = (\eta_{l+1} - \frac{3C_2}{2} \eta_{l+1}^2).
$$

Assuming that $\eta_0$ is sufficiently small, we have $\tilde{\eta}_{l+1} \ge \frac{1}{2} \eta_{l+1}$, $l = 0, \dots, L$. Therefore,

$$
\mathbb{P} \left( \bigcap_{l=1,\dots,L} \|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2 \ge \epsilon + \epsilon_{\mu} \right) \le \mathbb{P} \left( \sum_{l=0}^{L-1} \tilde{\eta}_{l+1} \|\nabla \hat{F}_{\text{uncorr}}(\theta^{(l)})\|_2^2 \ge \frac{1}{2} \eta_{1:L} (\epsilon + \epsilon_{\mu})^2 \right).
$$

Notice that

$$
\mathbb{P} \left( \bigcap_{l=1, \dots, L} \| \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) \|_2 \ge \epsilon + \epsilon_{\mu}, \right) \le \mathbb{P} \left( \bigcap_{l=1, \dots, L} \| \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) \|_2 \ge \epsilon + \epsilon_{\mu}, | \mathcal{E}_L \right) + 2Le^{-Kt^2} \quad (33)
$$

and combining Markov's inequality with (33), we have

$$
\mathbb{P} \left( \bigcap_{l=1, \dots, L} \| \nabla \hat{F}_{\text{uncorr}}(\theta^{(l)}) \|_2 \ge \epsilon + \epsilon_{\mu}, \right) \le 2 \frac{\mathbb{E} \left[ \| \nabla \hat{F}_{\text{uncorr}}(\theta^{(0)}) \|_2 | \mathcal{E}_L \right] + 3\eta_{1:L}^2 (\epsilon_{\mu} + B_1) + 3\eta_{1:L}\delta_0^2 d}{\eta_{1:L}(\epsilon + \epsilon_{\mu})^2} + 2Le^{-Kt^2}
$$
Since $\nabla \hat{F}_{uncorr}(\theta^{(0)})$ is independent of $\mathcal{E}_L$ and since

$$ \eta_{1:L} \geq \eta_0 L^{1/2} $$

and

$$ \eta_{1:L} \leq \log(L) \eta_0^2 $$

we get that

$$ \mathbb{P} \left( \bigcap_{l=1,...,L} \| \nabla \hat{F}_{uncorr}(\theta^{(l)}) \|_2 \geq \epsilon + \epsilon_{\mu}, \right) \leq 2 \frac{\mathbb{E}[\| \nabla \hat{F}_{uncorr}(\theta^{(0)}) \|_2]}{\eta_0 (\epsilon + \epsilon_{\mu})^2 \sqrt{L}} + 6\eta_0 \frac{\log(L) (\epsilon_{\mu} + B_1)}{(\epsilon + \epsilon_{\mu})^2 \sqrt{L}} + 6 \frac{\delta_0^2 d}{(\epsilon + \epsilon_{\mu})^2} + 2Le^{-Kt^2} $$

as announced.

## B ANISOTROPIC GRADIENT COVARIANCE

Consider the setting where data $X_i = (Z_i, Y_i)_{i=1}^n$ are generated from the linear model

$$ Y = Z\theta^* + \xi, \text{ where } Z \sim \mathcal{N}(0, \Sigma) \text{ and } \xi \sim \mathcal{N}(0, \sigma^2). \quad (34) $$

The gradient of the per-sample $\ell_i(\theta) = \frac{1}{2}(x_i^\top \theta - y_i)^2$ with respect to $\theta$ for sample $i$ is

$$ g_i = \nabla_{\theta} \ell_i(\theta) = (x_i^\top \theta - y_i)x_i = (x_i^\top \theta - x_i^\top \theta^* - \epsilon_i)x_i = x_i x_i^\top (\theta - \theta^*) - \epsilon_i x_i \quad (35) $$

Let $w = \theta - \theta^*$, then the covariance matrix of $g_i$ is

$$ \begin{aligned} \mathrm{Cov}[g_i] &= \mathrm{Cov}[x_i x_i^\top w - \epsilon_i x_i] \\ &= \mathrm{Cov}[x_i x_i^\top w] + \mathrm{Cov}[\epsilon_i x_i] \\ &= (w^\top \Sigma w)\Sigma + \Sigma w w^\top \Sigma + \sigma^2 \Sigma. \end{aligned} \quad (36) $$

Therefore, if $\Sigma = I_d$, then

$$ \mathrm{Cov}[g_i] = (\|w\|^2 + \sigma^2)I_d + ww^\top \quad (37) $$

which is mostly isotropic, up to a rank-one perturbation, which can be considered negligible when $\theta$ is close to the true parameter $\theta^*$. On the other hand, when $\Sigma$ is a general positive definite matrix, then it is easily seen from (36) that $\mathrm{Cov}[g_i]$ can be severely anisotropic.

## C IMPLEMENTATION DETAILS FOR REAL DATASET

**Binary CIFAR-10.** We design a binary classification task from CIFAR-10 by restricting the dataset to the classes “cat” (class id 3) and “dog” (class id 5). Images are normalized channel-wise using fixed values computed on the training set: mean $\mu = (0.4977, 0.4605, 0.4160)$ and standard deviation $\sigma = (0.2537, 0.2481, 0.2535)$.

**Contamination.** We simulate symmetric label noise at rate $\epsilon \in [0, 1]$. Let $Y \in \{0, 1\}$ denote labels with $Y = 0$ for cat and $Y = 1$ for dog. For each class, we uniformly select $[\epsilon \cdot |\mathcal{D}_{train}|/2]$ indices and flip their labels: $\tilde{Y}_i \leftarrow 1 - Y_i$. This yields a total contamination rate approximately $\epsilon$ and preserves class balance of flips. We test on $\epsilon \in \{0, 0.05, 0.1, 0.15, 0.2, 0.25\}$.

**Training.** We use the following parameters: learning rate $3 \times 10^{-4}$, 5000 maximum epochs, dataset batch size 2000, and mini-batch size 50 (used for all MoM-based methods). Early stopping is enabled with a validation split of 20% taken from the training set and a patience of 750 steps. We obtain 8000 training, 2000 validation and 2000 test samples. **Remark:** Although theoretically using early stopping from a contaminated dataset could lead to unexpected behavior, empirical results show that the validation loss follows closely the evolution of the test loss.
**Model.** We use a CNN for binary classification: three convolutional blocks with ReLU and 2 × 2 max pooling, followed by global average pooling and a two-layer MLP. Concretely: Conv(3→24, 3 × 3, pad 1) → ReLU → MaxPool → Conv(24→48, 3×3, pad 1) → ReLU → MaxPool → Conv(48→96, 3×3, pad 1) → ReLU → MaxPool → GAP → Linear(96→48) → ReLU → Linear(48→1). The model outputs a logit $z \in \mathbb{R}$ for the positive class (dog).

**Loss and metrics.** Training uses the binary cross-entropy with logits,

$$ \mathcal{L}_{\text{BCE}}(z, y) = -y \log \sigma(z) - (1-y) \log(1-\sigma(z)) $$

with $\sigma(z) = 1/(1 + e^{-z})$. At evaluation, we report error rate on the test set,

$$ \text{Error} = 1 - \frac{1}{N} \sum_{i=1}^{N} 1\{\sigma(z_i) > 0.5 \Leftrightarrow y_i = 1\}, $$

**TMoM minibatching.** For TMoM, we separate the samples of each batch based on their classes and form minibatches from each class. Then we compute the TMoM update on these minibatches.

For other MoM methods, we split each batch in mini-batches directly, as per the original implementations.