Bayesian Robustness: A Nonasymptotic Viewpoint

TL;DR

This paper introduces Rob-ULA, a robust sampling algorithm for Bayesian inference that effectively handles adversarial outliers, providing finite-sample guarantees and demonstrating practical performance on various tasks.

Contribution

The paper proposes Rob-ULA, a novel robust variant of ULA with finite-sample analysis, addressing adversarial contamination in Bayesian posterior sampling.

Findings

Rob-ULA achieves accurate sampling within . .  iterations.

Theoretical bounds relate sample complexity to data dimension and contamination level.

Experimental results validate robustness and effectiveness on real-world datasets.

Abstract

We study the problem of robustly estimating the posterior distribution for the setting where observed data can be contaminated with potentially adversarial outliers. We propose Rob-ULA, a robust variant of the Unadjusted Langevin Algorithm (ULA), and provide a finite-sample analysis of its sampling distribution. In particular, we show that after $T= \tilde{\mathcal{O}}(d/\varepsilon_{\textsf{acc}})$ iterations, we can sample from $p_T$ such that $\text{dist}(p_T, p^*) \leq \varepsilon_{\textsf{acc}} + \tilde{\mathcal{O}}(\epsilon)$, where $\epsilon$ is the fraction of corruptions. We corroborate our theoretical analysis with experiments on both synthetic and real-world data sets for mean estimation, regression and binary classification.

Equations

$$p(\theta |\mathcal{D},\alpha )\propto p_{\theta }(\theta |\alpha )\cdot \prod _{i=1}^{n}p(z_i|\theta ).$$

$$q(h|\mathcal{D},\alpha ):={\int }_{{\mathbb{R}}^d}h(\theta )p(\theta |\mathcal{D},\alpha )\mathrm{d}\theta .$$

$$\hat{q}(h|\mathcal{D},\alpha )=\frac{1}{T}\sum _{t=1}^{T}h({\theta }_t).$$

$$p^{\ast }(\theta )\propto \exp (-f(\theta )),$$

$${\theta }_{k+1,\eta }={\theta }_{k,\eta }-{\eta }_{k+1}\nabla f({\theta }_{k,\eta })+\sqrt{2{\eta }_{k+1}}{\xi }_{k+1},$$

$$\mathrm{d}{\theta }_t=-\nabla f({\theta }_t)\mathrm{d}t+\sqrt{2}\mathrm{d}{B}_t,t\ge 0,$$

$$p(\theta |\mathcal{D},\alpha )\propto \exp \left(\log (p_{\theta }(\theta |\alpha ))-\sum _{i=1}^n{g}_i(\theta )\right),$$

$${\theta }_{k+1,\eta }={\theta }_{k,\eta }-{\eta }_{k+1}\left(-\nabla \log (p_{\theta }({\theta }_{k,\eta }|\alpha )+\sum _{i=1}^n\nabla g_i({\theta }_{k,\eta })\right)+\sqrt{2{\eta }_{k+1}}{\xi }_{t+1},$$

$$z_i\sim (1-ϵ)P+ϵQ,$$

$$p(\theta |{\mathcal{D}}_c,\alpha )\propto p_{\theta }(\theta |\alpha )\cdot \prod _{i\in {\mathcal{D}}_c}(z_i|\theta ).$$

$$q(h|\mathcal{D},\alpha ):={\int }_{{\mathbb{R}}^d}h(\theta )p(\theta |{\mathcal{D}}_c,\alpha )\mathrm{d}\theta ,$$

$$\hat{q}(h|{\mathcal{D}}_c,\alpha )=\frac{1}{T}\sum _{t=1}^{T}h({\theta }_t^c).$$

$$p(\theta |{\mathcal{D}}_c,\alpha )\propto p_{\theta }(\theta |\alpha )\cdot \prod _{i\in {\mathcal{D}}_c}(z_i|\theta ).$$

$$ϵ-\sqrt{\frac{2}{n}\log \left(\frac{1}{\delta }\right)}\le {ϵ}_n\le ϵ+\underset{e_n}{\underbrace{\sqrt{\frac{2}{n}\log \left(\frac{1}{\delta }\right)}}}.$$

$${\Vert \hat{\nabla }{U}_{\theta }-\frac{1}{|{\mathcal{D}}_c|}\sum _{i\in {\mathcal{D}}_c}\nabla g_i(\theta )\Vert }_2,$$

$$\Vert \hat{\nabla }{U}_{\theta }-\nabla U_{\theta }{\Vert }_2\le C_1{C}_4^{\frac{1}{4}}\sqrt{\gamma \log (d)\Vert {\Sigma }_{\theta }{\Vert }_2}.$$

$$\text{KL}(p(\theta )\Vert p^{\ast }(\theta ))=\int p(\theta )\ln \left(\frac{p(\theta )}{p^{\ast }(\theta )}\right)\mathrm{d}\theta .$$

$${\text{W}}_2^2(p,q):=\underset{\zeta \in \Gamma (p,q)}{\inf }\int \Vert x-y{\Vert }_2^2\mathrm{d}\zeta (x,y),$$

$$\Vert \nabla f(\theta )-\nabla f(\nu )\Vert \le L\Vert \theta -\nu \Vert ,\forall \theta ,\nu \in {\mathbb{R}}^d\text{and}{\nabla }^{2}f(\theta )\text{ exists for all }\theta \in {\mathbb{R}}^d.$$

$$mI⪯{\nabla }^{2}f(\theta ),\forall \theta \in {\mathbb{R}}^d.$$

$${\Vert \nabla f({\theta }_{k,\eta })-\hat{\nabla }f({\theta }_{k,\eta })\Vert }^2\le n^2ϵC_R{\Vert {\Sigma }_{\theta }\Vert }_2\log d\text{and}{\Vert {\Sigma }_{\theta }\Vert }_2\le C_{\Sigma ,1}{\Vert \theta -\tilde{\theta }\Vert }^2+C_{\Sigma ,2},$$

$${\text{W}}_2^2(p_{k\eta },p^{\ast })$$

$$+ϵ\left(\frac{4C_R{C}_{\Sigma ,2}}{{\bar{m}}^2}\log d+\frac{8C_R{C}_{\Sigma ,1}}{{\bar{m}}^2}\frac{d\log d}{n}\right)$$

$$\le \frac{1}{\bar{m}}\log \frac{\bar{L}}{\bar{m}}\frac{d}{n}{e}^{-\bar{m}kh}+8\left(C_R{C}_{\Sigma ,2}\frac{{\bar{L}}^4}{{\bar{m}}^4}ϵ\log d+\frac{{\bar{L}}^4}{{\bar{m}}^3}\frac{d}{n}\right)h^2+4\frac{{\bar{L}}^2}{{\bar{m}}^2}\frac{d}{n}h$$

$$+ϵ\left(\frac{4C_R{C}_{\Sigma ,2}}{{\bar{m}}^2}\log d+\frac{8C_R{C}_{\Sigma ,1}}{{\bar{m}}^2}\frac{d\log d}{n}\right),$$

$$\underset{(I)}{\underbrace{\frac{2e^{-n\bar{m}k\eta }}{n\bar{m}}\text{KL}(p_0\Vert p^{\ast })}}+\underset{(II)}{\underbrace{C\left(\frac{{\bar{L}}^4}{{\bar{m}}^4}ϵ\log d+\frac{{\bar{L}}^4}{{\bar{m}}^3}\frac{d}{n}\right)h^2+4\frac{{\bar{L}}^2}{{\bar{m}}^2}\frac{d}{n}h}}+\underset{(III)}{\underbrace{\frac{C}{{\bar{m}}^2}ϵ\log d}},$$

$$\eta =\mathcal{O}\left(\frac{{\epsilon }_{\text{acc}}}{n\kappa \bar{L}d}\right)\text{and}T\ge \mathcal{O}\left(\frac{\bar{L}}{\bar{m}}\log \left(\frac{\text{KL}(p_0\Vert p^{\ast })}{n\bar{m}{\epsilon }_{\text{acc}}}\right)\right),$$

$${\text{W}}_2^2(p_{T,\eta },p^{\ast })\le {\epsilon }_{\text{acc}}+\tilde{\mathcal{O}}\left(\frac{ϵ}{{\bar{m}}^2}\right).$$

$$\text{KL}(p_0\Vert p^{\ast })=\int p_0(\mathbf{x})\log \left(\frac{p_0(\mathbf{x})}{p^{\ast }(\mathbf{x})}\right)\mathrm{d}\mathbf{x}\le \frac{d}{2}\log \frac{\bar{L}}{\bar{m}}.$$

$${\Theta }_{(k+1)\eta }={\Theta }_{k\eta }-\eta \hat{\nabla }f({\Theta }_{k\eta })+\sqrt{2}(B_{(k+1)\eta }-B_{k\eta }),$$

$$\mathrm{d}{\Theta }_t=-\hat{\nabla }f({\Theta }_{kh})\mathrm{d}t+\sqrt{2}\mathrm{d}{B}_t,\text{ for }k\eta <t\le (k+1)\eta .$$

$$\mathbb{E}\left[{\Vert {\Theta }_{k\eta }-\tilde{\theta }\Vert }_2^2\right]\le \frac{4}{{\bar{m}}^2}{C}^{\prime }ϵ\log d+\frac{4d}{n\bar{m}},$$

Full text

1 Introduction

Robustness has been of ongoing interest in both the Bayesian [8, 3] and frequentist setting [36, 15, 16] since being introduced by George Box in 1953 [4]. The goal is to capture the sensitivity of inferential procedures to the presence of outliers in the data and misspecifications in the modelling assumptions, and to mitigate overly large sensitivity. The Bayesian approach has been focused on capturing possible anomalies in the observed data via the model and in choosing priors that have minimal effect on inferences. The frequentist approach, on the other hand, has focused on the development of estimators that identify and guard against outliers in the data. We refer the reader to [18, Chap 15] for a comprehensive discussion.

The focus on model robustness in Bayesian statistics is implemented via sensitivity studies to understand effects of misspecification of the prior distribution [3, 24] and its propagation towards the posterior [17]. There is, however, little in the way of a comprehensive formal finite-sample framework for Bayesian robustness. Huber asked “Why there is no finite sample Bayesian robustness theory?” and Kong suggested that such a theory is infeasible in full generality, arguing that it is computationally infeasible to carry out the necessary calculations even in finite spaces.

We address this issue by providing a formal framework for studying Bayesian robustness and by proposing a robust inferential procedure with finite-sample guarantees. We address issues of computational infeasibility by refraining from modelling outlier data explicitly. Instead, we posit that the collected data contains a small fraction of observations which are not explained by the modelling assumptions. This corruption model, termed an $\epsilon$-contamination model, was first studied by Huber in 1964 [15] and has been the subject of recent computationally-focused work in the frequentist setting [10, 20, 30].

Given data corrupted according to an $\epsilon$-contamination model, our goal is to sample from the clean posterior distribution $p^{*}$: the posterior distribution conditioning only on the uncorrupted (“clean”) data. Our key idea is to leverage a robust approach for estimating the mean in the context of gradient-based sampling techniques. Our overall procedure is a robust variant of the Unadjusted Langevin Algorithm (ULA) that we refer to as “Rob-ULA.” The underlying ULA algorithm and its variants have been used for efficient large-scale Bayesian posterior sampling [32, 39] and their convergence analysis has been a recent topic of interest [6, 12, 5, 7, 22]; see Section 2.2 for a detailed overview. Informally, our main result shows that after $T=\tilde{\mathcal{O}}(d/\varepsilon_{\textsf{acc}})$ iterations of Rob-ULA, the iterate $\theta_{T}$ has a distribution $p_{T}$ such that $\text{dist}(p_{T},p^{*})\leq\varepsilon_{\textsf{acc}}+\tilde{\mathcal{O}}(\epsilon)$, where $\epsilon$ is the fraction of corrupted points in the data set.

The remainder of the paper is organized as follows: Section 2 contains a discussion of the related literature, Section 3 discusses relevant background as well as the formal problem setup, and Section 4 describes the proposed algorithm Rob-ULA and states our main theorem regarding its convergence. In Section 5 we discuss the fundamental problems of Bayesian mean estimation and linear regression in the robust setting. Section 6 consists of experimental evaluation of Rob-ULA on synthetic as well as real world data sets and we conclude with Section 7.

2 Related Work

In this section, we review related work on robust estimation procedures in both the Bayesian and frequentist settings. We also discuss work on using Langevin dynamics to sample from distributions over continuous spaces.

2.1 Robust statistical procedures

There are many threads in the literature on robust estimation and outlier detection [16, 4, 8]. In the frequentist parameter estimation setting, the most commonly studied model is Huber’s classical $\epsilon$-contamination model. There has also been a recent focus on an adversarial paradigm that is devoted to developing computationally efficient problem-dependent estimators for mean estimation  [20, 10], linear regression [19, 1, 2, 33], and general risk minimization problems [30, 11]. Particular relevant to our setup are [30] and [11] which utilize the robust mean estimators of  [20, 10] to robustify gradient-based procedures for empirical risk minimization.

The study of robustness in the Bayesian framework has focused primarily on developing models and priors that have minimal impact on inference. An important line of work has focused on the sensitivity of the posterior distribution to and has led to the development of noninformative priors [3, 24, 23]. These methods are orthogonal to those considered in the current paper, as they do not aim to robustify inferential procedures against corruptions in the observed data set. In principle a complete Bayesian model would have the capacity for explaining the outliers present in the data, but this would require performing an integral over all models with a proper prior. Such an approach would generally be computationally intractable.

An important class of procedures for handling outliers in the data set focuses on reweighing the data points via a transformation of the likelihood function. Huber [18] considers assigning binary weights to data points and identifies model-dependent procedures to identify outliers. In contrast, Wang et al. [38] consider these weights as latent variables and infers these weight variables along with the parameters of the model. These methods are susceptible to the choice of priors over these weighting variables. An alternate set of robust procedures are based on the idea of localization [8, 37]: each data point is allowed to depend on its own version of the latent variable and these individual latent variables are tied together with the use of a hyperparameter, often fitted using empirical Bayes estimation. Although these methods build on intuitive notions of outlier removal, there is little theoretical understanding of the kind of outliers that these methods can tolerate.

2.2 Sampling methods

Turning to sampling methods for posterior inference, there have been various zeroth-order [21, 25] and first-order methods [13, 31, 27] proposed for sampling from distributions over continuous spaces. Our focus in this paper is on the overdamped Langevin MCMC method which was first proposed by Ermak [13] in the context of molecular dynamics applications. Its nonasymptotic convergence (in total variation distance) was first studied by Dalalyan [6] for log-smooth and log-strongly concave distributions. Cheng and Bartlett [5] extended the analysis to obtain a similar convergence result in Kullback-Leibler divergence. Such nonasymptotic convergence guarantees are essential to understanding the robustness of computational procedures as they simultaneously capture the dependence of the sampling error on the number of iterations $T$, the dimensionality $d$, and the contamination level $\epsilon$.

3 Background

In the section we briefly review relevant background on Bayesian computation and we formally describe our problem setup.

3.1 Bayesian modelling

Given parameters $\theta\in\mathbb{R}^{d}$ and a data set $\mathcal{D}=\{z_{1},z_{2},\ldots,z_{n}\}$, we assume that the statistician has specified a prior distribution, $p_{\theta}(\theta|\alpha)$, and a likelihood, $p(z|\theta)$. We can then form the posterior distribution, $p(\theta|\mathcal{D},\alpha)$, as follows:

[TABLE]

We are generally concerned with the estimation of some test function $h(\theta)$ under the posterior distribution, which is accomplished in the Bayesian setting by computing a posterior mean:

[TABLE]

In practice, one way of computing this posterior mean is to use a Monte Carlo algorithm to generate a sequence of samples $\{\theta_{t}\}_{t=1}^{T}$ and form an estimate $\hat{q}(h|\mathcal{D},\alpha)$:

[TABLE]

3.2 Unadjusted Langevin Algorithm

We consider a specific Monte Carlo algorithm, the Unadjusted Langevin Algorithm (ULA), for sampling from probability distributions over continuous spaces. Generically, we consider distributions over $\mathbb{R}^{d}$ of the form

[TABLE]

for a class of functions $f$ which are square integrable. The ULA algorithm starts from some initial point $\theta_{0}\in\mathbb{R}^{d}$ and defines a sequence of parameters, $\{\theta_{k,\eta}\}_{k=1}^{T}$, according to the following update equation:

[TABLE]

where $\eta=\{\eta_{i}>0\}$ denotes a sequence of step sizes and $\{\xi_{i}\}_{i\in\mathbb{N}}\sim\mathcal{N}(0,\text{I}_{d\times d})$ are i.i.d. Gaussian vectors. The Markov chain in Equation (1) is the Euler discretization of a continuous-time diffusion process $\{\theta_{t}\}_{t\geq 0}$ known as the Langevin diffusion. The stochastic differential equation governing the Langevin diffusion is given by

[TABLE]

where $\{B_{t}\}_{t\geq 0}$ represents a $d$-dimensional Brownian motion. Denoting the distribution of $\theta_{k,\eta}$ by $p_{k,\eta}$, Cheng and Bartlett [5] showed that $\text{KL}(p_{k,\eta}\;\|\;p^{*})\leq\epsilon$ after $t=\tilde{\mathcal{O}}(\frac{d}{\epsilon})$ steps for functions $f$ which are smooth and strongly-convex. Specializing to the Bayesian modelling setup we rewrite the posterior distribution as:

[TABLE]

where $g_{i}(\theta):=-\log(p(z_{i}|\theta))$ is the negative log-likelihood corresponding to the $i^{th}$ data point. The ULA algorithm can then be used to form an approximation to the posterior as follows:

[TABLE]

where $\eta$ and $\xi_{t+1}$ are the step-size sequence and independent Gaussian noise respectively.

3.3 Problem Setup

We turn to a formal treatment of the robustness problem in the Bayesian setting. We consider the $\epsilon$-contamination model introduced by Huber [15] and let the collection of $n$ data points be obtained from the following mixture distribution:

[TABLE]

where $P$ denotes the true underlying generative distribution while $Q$ is any arbitrary distribution. A data set $\mathcal{D}$ drawn from such a mixture distribution has each data point $z_{i}$ corrupted adversarially with probability $\epsilon$. We denote by $\mathcal{D}_{c}$ the subset of data points in $\mathcal{D}$ sampled from the true distribution $P$ and similarly let $\mathcal{D}_{a}$ denote the subset of data sampled from $Q$. Given data points $\mathcal{D}=\mathcal{D}_{c}\cup\mathcal{D}_{a}$, the likelihood function $p(z|\theta)$ and the prior $p_{\theta}(\theta|\alpha)$, we aim to form a clean posterior distribution given by:

[TABLE]

Accordingly, as in Section 3.1, we would like to robustly estimate the mean of the test function $h(\theta)$ under the uncorrupted posterior $p(\theta|\mathcal{D}_{c},\alpha)$:

[TABLE]

which we approximate via an estimate:

[TABLE]

In the following section we present an algorithm, Rob-ULA, for generating the sequence of samples $\{\theta^{c}_{t}\}_{t=1}^{T}$ and provide theoretical guarantees on its convergence properties. The main idea is to exploit gradient-based iterative sampling methods, and to leverage a robust mean estimation procedure to compute a robust estimate of the gradient at each iteration.

4 Rob-ULA: Robust Unadjusted Langevin Algorithm

We turn to our proposed algorithm, Rob-ULA (Algorithm 1), which aims to solve the robust posterior inference problem defined in Section 3.3. Rob-ULA is a simple modification of the ULA algorithm, described in Section 3.2, where in each iteration instead of using the complete set of data points for computing the gradient, we construct a robust estimator of the gradient and update the parameter using this estimate. This robust estimator ensures that the outlier data points do not exert too much influence on the gradient and allow Rob-ULA to obtain samples from a distribution close to the clean posterior distribution:

[TABLE]

Before proceeding to establish the convergence guarantees for Rob-ULA, we present the robust gradient estimation procedure.

4.1 Robust Gradient Estimation

Algorithm 2 describes our robust gradient estimation procedure. Based on the robust mean estimator of Lai et al. [20], it takes as input the gradients of the negative log-likelihoods $\nabla g_{i}(\theta)$ and outputs an estimate of the robust mean of the gradient vectors ($\widehat{\nabla}U_{\theta}$ in Algorithm 1), assuming a fraction $\epsilon$ of them are arbitrarily corrupted. Algorithm 1 then scales this gradient estimate by the number of samples $n$, to obtain a robust estimate of gradients of the likelihood $\sum_{i=1}^{n}\nabla g_{i}(\theta)$.

Note that the model described in Section 3.3 assumes that each data point $z$ is sampled i.i.d. from the mixture distribution $(1-\epsilon)P+\epsilon Q$, where $P$ represents the true generative distribution and $Q$ can be any arbitrary distribution. An application of the Hoeffding bound for Bernoulli random variables shows that with probability at least $1-\delta$, the fraction of corrupted points $\epsilon_{n}$ in the sampled data set $\mathcal{D}$ satisfy

[TABLE]

For the remainder of the paper, we condition on this high probability event and state our results assuming this event holds. Following the proof strategy of [20] and [30], we derive a bound on the estimation error of the true average log-likelihood gradient,

[TABLE]

uniformly for any value of the iterate $\theta$ in the following lemma. We let $\nabla U_{\theta}:=\frac{1}{|\mathcal{D}_{c}|}\sum_{i\in\mathcal{D}_{c}}\nabla g_{i}(\theta)$ denote the true value of this average log-likelihood gradient.

Lemma 1** (Robust Gradient Estimation)**

Let $P$ denote the uniform distribution over $\mathcal{D}_{c}$ and let $P_{\theta}$ denote the corresponding distribution over $\nabla g_{i}(\theta)$ with mean given by $\nabla U_{\theta}$, covariance $\Sigma_{\theta}$ and fourth moment given by $C_{4}$. There exists a positive constant $C_{1}>0$ for which the robust mean estimator when instantiated with the contamination level $\gamma:=\epsilon+e_{n}$, returns, with probability $1-\delta$, an estimate $\widehat{\nabla}U_{\theta}$ such that for all $\theta\in^{d}$, we have that,

[TABLE]

Remark.

Note that Proposition 1 of Prasad et al. [30] presents a high-probability bound similar to ours which is applicable for a fixed parameter $\theta$. Such a bound, however, does not suffice to ensure convergence of Rob-ULA because the additive Gaussian noise at every iterate requires us to obtain a uniform high-probability recovery error bound (see Section 4.2 for details). Lemma 1 establishes a uniform bound for the specific distribution $P$ which is uniform over the clean data $\mathcal{D}_{c}$. This restriction of the distribution $P$ also allows us to avoid sample-splitting at every iteration of Algorithm 1 which was essential for both [20] and [30].

In addition, Lemma 1 indicates that irrespective of the sample size $n$, one can estimate the mean of the gradient robustly up to error $\tilde{\mathcal{O}}(\sqrt{\epsilon\|\Sigma_{\theta}\|_{2}})$. This implies that at each iteration of Rob-ULA, we incur an error of $\tilde{\mathcal{O}}(n\sqrt{\epsilon\|\Sigma_{\theta}\|_{2}})$ since we scale the average gradient estimate by $n$ during the update. In Theorem 2 we show how with an appropriate choice of step size $\eta=\mathcal{O}(1/n)$, one can control the propagation of this bias in the convergence analysis for Rob-ULA.

The detailed proof of Lemma 1 can be found in Appendix A.

4.2 Convergence Analysis

In this section, we study the convergence of the proposed algorithm Rob-ULA. For ease of notation, we let $f(\theta;\mathcal{D})=\sum_{i\in\mathcal{D}}g_{i}(\theta)-\log(p_{\theta}(\theta|\alpha))$ and similarly denote the clean and corrupted versions of the function $f(\theta;\mathcal{D}_{c})$ and $f(\theta;\mathcal{D}_{a})$. The objective of the robust Bayesian posterior inference problem is then to obtain samples from the clean posterior distribution given by $p^{*}(\theta|\mathcal{D},\alpha)\propto\exp\left(-f(\theta;\mathcal{D}_{c})\right)$. For clarity of exposition, we drop the dependence of the posterior distribution on the data set $\mathcal{D}$ as well as the hyperparameters $\alpha$ and let $f(\theta):=f(\theta;\mathcal{D}_{c})$.

We quantify the convergence of distribution $p(\theta)$ following a stochastic process to the stationary distribution $p^{*}(\theta)$ through the Kullback-Leibler divergence, $\text{KL}(p(\theta)\;\|\;p^{*}(\theta))$:

[TABLE]

We define the Wasserstein distance $\text{W}_{2}^{2}(p,q)$ between a pair of distributions $(p,q)$ as:

[TABLE]

where $\Gamma(p,q)$ denotes the set of joint distributions such that the first set of coordinates has marginal $p$ and the second set has marginal $q$.

We begin by making the following assumptions on the function $f(\theta)$:

• **Assumption 1 **(Lipschitz smoothness). The function $f(\theta)$ is $L$-Lipschitz smooth and its Hessian exists for all $\theta\in\mathbb{R}^{d}$. That is,

[TABLE]

• **Assumption 2 **(Strong convexity). The function $f(\theta)$ is $m$-strongly convex for all $\theta\in^{d}$. That is,

[TABLE]

We further denote the condition number of the function $f$ as $\kappa=L/m$.

The assumptions of Lipschitz smoothness and strong convexity are standard in both the sampling and optimization literatures. In addition to the assumptions, we define the average Lipschitz constant $\bar{L}=L/n$ and the average strong convexity of $f$ as $\bar{m}=m/n$. We now state our main theorem concerning the convergence guarantees for Rob-ULA.

Theorem 2** (Main Result)**

Let $p^{*}(\theta)\propto\exp(-f(\theta))$, where $f$ satisfies Assumptions 4.2 and 4.2. Further, assume that the gradient estimates $\widehat{\nabla}f(\theta)$ satisfy

[TABLE]

*where $\Sigma_{\theta}$ is the covariance of uniform distribution on $\nabla g_{i}(\theta)$ induced by the clean data set $\mathcal{D}_{c}$, $\tilde{\theta}$ satisfies $\nabla F(\tilde{\theta})=0$ and $\epsilon$ is the fraction of corrupted points, satisfying $\epsilon\leq\dfrac{\bar{m}^{2}}{4C_{R}C_{\small{\Sigma,1}}\log d}$. Then the iterates of Rob-ULA, when initialized with $\theta_{0}\sim\mathcal{N}\left(0,\frac{1}{L}I_{d}\right)$ (with corresponding density $p_{0}$) and step size $\eta\leq\frac{1}{n\bar{L}}$ (and define $h:=n\eta\leq\frac{1}{\bar{L}}$), satisfy:

---

[TABLE]

where $p_{k\eta}$ represents the distribution of the iterate $\theta_{k,\eta}$.

Remarks.

Before proceeding to the proof of this theorem, a few comments are in order. First observe that the error term consists of three different components:

[TABLE]

where a) term (I) comprises an exponentially decaying dependence (with the number of time-steps $t$) on the initial error $\text{KL}(p_{0}\;\|\;p^{*}(\theta))$, b) term (II) is a discretization error term and c) term (III) captures the dependence on the fraction of corrupted points $\epsilon$ and vanishes as $\epsilon$ goes to zero.

For any given accuracy $\varepsilon_{\textsf{acc}}$, if the step size and the number of iterations satisfy:

[TABLE]

then the error in convergence can be bounded as

[TABLE]

As we show in Section 5, for problems such as Bayesian linear regression and Bayesian mean estimation, the average strong convexity parameter $\bar{m}$ scales independently of the sample size $n$. This implies that the resulting error can be bounded by $\varepsilon_{\textsf{acc}}+\tilde{\mathcal{O}}(\epsilon)$. While the accuracy can be set to arbitrarily small values which would result in a corresponding increase in the number of time steps, there is a bias term depending on the contamination level $\mathcal{O}(\epsilon)$ which cannot be reduced by either increasing the sample size or by increasing the number of iterations. This is consistent with results in the frequentist literature [1, 10, 20, 30], which show that such inconsistency is a result of the adversarial corruptions and in general cannot be avoided.

Lemma B.20 in Appendix B presents the following bound on the initial error:

[TABLE]

Proof of Theorem 2

We proceed to a proof of our main convergence result, Theorem 2. We begin by considering the process described by Rob-ULA as a discretization of the Langevin dynamics given by Equation 2, with the following gradient estimate:

[TABLE]

where ${{\Theta}}_{k\eta}$ represents the random variable describing the process at the $k^{th}$ iterate using step size $\eta$. This is equivalent to defining the following stochastic differential equation

[TABLE]

We next state a lemma which provides a bound on the variance of the distribution of the $k^{th}$ iterate.

Lemma 3

For $\Theta_{t}$ following Eq. (7), if $\Theta_{0}\sim\mathcal{N}\left(0,\dfrac{1}{n\bar{L}}\mathrm{I}\right)$, $\epsilon\leq\dfrac{\bar{m}^{2}}{4C\log d}$, and $h:=n\eta\leq\dfrac{1}{\bar{L}}$, then for all $k\in\mathbb{N}^{+}$,

[TABLE]

where $C$ is a universal constant and $\tilde{\theta}$ is the fixed point satisfying $\nabla f(\tilde{\theta})=0$.

The proof of this lemma is deferred to Appendix B. Treating Lemma 3 as given, we proceed to the proof of Theorem 2.

We consider the dynamics in Equation 7 within the time range $k\eta<t\leq(k+1)\eta$. From the Girsanov theorem [28] we have that ${{\Theta}}_{t}$ admits a density function $p_{t}$ with respect to the Lebesgue measure. This density function can also be represented as

[TABLE]

where $p({{\theta}},t|{{\zeta}},k\eta)$ is the weak solution to the following Kolmogorov forward equation:

[TABLE]

where $p({{\theta}},t|{{\zeta}},k\eta)$ and its derivatives are defined via $P_{t}(f)=\int f({{\theta}})p({{\theta}},t|{{\zeta}},k\eta)\mathrm{d}{{\theta}}$ as a functional over the space of smooth bounded functions on $\mathbb{R}^{d}$ (we refer the readers to  [34] for more details). As shown by Cheng and Bartlett [5], the time derivative of the KL divergence along $p_{t}$ is given by:

[TABLE]

where the expectation is taken with respect to the joint distribution of ${{\Theta}}_{t}$ and ${{\Theta}}_{k\eta}$. Hence

[TABLE]

where ${\sf(i)}$ follows from the definition of $p^{*}(\theta)\propto\exp(-f(\theta))$. We first focus on the second term in the above expression which can be bounded as:

[TABLE]

where ${\sf(i)}$ follows by an application of Young’s inequality $2a^{\top}b\leq\|a\|_{2}^{2}+\|b\|_{2}^{2}$ and ${\sf(ii)}$ follows from the point-wise assumption that $\left\|\nabla f({{\theta}}_{k,\eta})-\widehat{\nabla}{f}({{\theta}}_{k,\eta})\right\|^{2}\leq n^{2}\epsilon C_{R}\left\|\Sigma_{{\theta}}\right\|_{2}\log d$ and that $\left\|\Sigma_{{\theta}}\right\|_{2}\leq C_{\small{\Sigma,1}}\left\|{{\theta}}-\tilde{\theta}\right\|^{2}+C_{\small{\Sigma,2}}$. Let us define new constant $C_{13}:=C_{R}\cdot C_{\small{\Sigma,1}}$ and $C_{14}:=C_{R}\cdot C_{\small{\Sigma,2}}$.

Next, we proceed to bound the term $\mathbb{E}\left[\left\|{{\Theta}}_{t}-{{\Theta}}_{k\eta}\right\|^{2}\right]$ using Lemma 3. Let us define the variable $\tau:=t-k\eta\in(0,\eta]$ and bound the term as:

[TABLE]

where ${\sf(i)}$ follows from Assumption 4.2 and we define $\nu=n\tau$ ($h=n\eta$). Plugging in the bounds obtained in Equations (4) and (4) into Equation (4), we get for $k\eta<t\leq(k+1)\eta$:

[TABLE]

where ${\sf(i)}$ follows from an application of the log-Sobolev inequality with $m$ being the log-Sobolev constant and ${\sf(ii)}$ follows from the fact that $\epsilon\leq\dfrac{\bar{m}^{2}}{4C_{13}\log d}$ and the substitution

[TABLE]

Finally, using Gronwall’s inequality we have the following one-step progress equation:

[TABLE]

Repeated application of this progress inequality leads us to

[TABLE]

The final result of the theorem can then be obtained by using Talagrand’s inequality [29] which states that for the probability distributions $p_{k\eta}$ and ${p^{*}}$, we have that

[TABLE]

which concludes the proof of the theorem. $\blacksquare$

5 Consequences for Mean Estimation and Regression

In this section, we study the fundamental problems of Bayesian mean estimation (Section 5.1) and Bayesian linear regression (Section 5.2) under the Huber $\epsilon$-contamination model.

5.1 Robust Bayesian mean estimation

We begin with the robust Bayesian mean estimation (RBME) problem and instantiate the convergence guarantees for Rob-ULA (Algorithm 1) for this problem. For simplicity, we study the setup in which the likelihood is Gaussian:

[TABLE]

for a mean vector ${{\theta}}\in\mathbb{R}^{d}$ and a fixed positive definite covariance matrix $\Sigma\in\mathbb{R}^{d\times d}$. We consider the corresponding conjugate prior over ${{\theta}}$ given by

[TABLE]

where $\theta_{0}$ is the mean and $\Sigma_{0}\succ 0$ is the covariance matrix. The set of parameters $(\theta_{0},\Sigma_{0})$ are the hyperparameters $\alpha$. Given data points $\mathcal{D}=\{z_{i}\}_{i=1}^{n}$ sampled from the Huber contamination model, where $Z_{i}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}(1-\epsilon)P+\epsilon Q$, with $Q$ being an arbitrary adversarially chosen distribution, the objective of the RBME problem is to sample from the posterior induced by the uncorrupted data points,

[TABLE]

where $\mathcal{D}_{c}$ represents the subset of points in $\mathcal{D}$ sampled from the distribution $P$. We note that for data $X_{i}$ sampled from the Huber contamination model, we have from Equation (5) that with probability at least $1-\delta$,

[TABLE]

Let us denote by $\bar{\epsilon}:=\epsilon+e_{n}$ and by $\underline{\epsilon}:=\epsilon-e_{n}$ the corresponding upper and lower bounds on the number of corrupted data points. Following the notation in Section 4, the function $f(\theta)$ for this problem is given by:

[TABLE]

Also, we define the corresponding function $g_{i}(\theta)$ from Section 3 as well as the sample mean and covariance for the clean data points.

[TABLE]

The following corollary instantiates the guarantees of Theorem 2 for the specific function $f({{\theta}})$ defined for the RBME problem.

Corollary 4

*Consider the RBME problem with posterior given by Equation (11) and data sampled from the Huber model. Then the iterates of Rob-ULA with step size $\eta\leq\frac{1}{n\bar{L}}$ and $h:=n\eta$ satisfy:

---

[TABLE]

where $C$ is a constant depending on the fourth moment of the clean data $\mathcal{D}_{c}$ and

[TABLE]

with probability at least $1-\delta$ for $\epsilon<1/2$.

Remark.

Observe that the step-size parameter $h$ is independent of $n$ while both $\bar{L}$ and $\bar{m}$ are asymptotically independent of the sample size $n$. The above bound shows that for an appropriately chosen step size one can obtain samples from a distribution which is $\tilde{\mathcal{O}}(\epsilon)$ away from the true distribution $p^{*}$. The number of iterations required to obtain such a sample scales linearly with the number of samples $n$, the dimension $d$ and the condition number $\kappa=\frac{\bar{L}}{\bar{m}}$.

5.2 Robust Bayesian linear regression

We turn to the robust Bayesian Linear Regression (RBLR) problem. For this problem, we let the data set $\mathcal{D}=\{z_{i}=(x_{i},y_{i})\}_{i=1}^{n}$ be such that $x_{i}\in^{d}$ and $y_{i}\in$ be the covariate vectors and response variables sampled from the Huber contamination model. Note that the Huber contamination model is on the variable $z_{i}$ and hence allows for corruption in both the features $x_{i}$ as well as the response variables $y_{i}$. In addition, we assume that there exists a vector ${{\theta}}^{*}$ such that

[TABLE]

where $z_{i}\sim\mathcal{N}(0,\sigma^{2})$ are sampled independently of $x_{i}$. This assumption is for simplifying the presentation and in general one can work with $\theta^{*}$ which is the best linear approximation to the data. For the RBLR problem, we consider likelihood functions of the form:

[TABLE]

for a fixed variance parameter $\sigma^{2}$. Also, we consider a Gaussian prior over the parameter ${{\theta}}$,

[TABLE]

for some fixed mean vector ${{\theta}}_{0}$ and positive-definite covariance matrix $\Sigma_{0}$ which form the set of hyperparameters $\alpha$. Given a data set $\mathcal{D}$ sampled from the Huber $\epsilon$-contamination model, the objective of the RBLR problem is to sample from the posterior induced by the uncorrupted set of data points,

[TABLE]

Following a similar calculation to that in Section 5.1, we have that with probability at least $1-\delta$,

[TABLE]

The corresponding function $f({{\theta}})$ for the RBLR problem is then defined to be

[TABLE]

We denote by ${{\theta}}_{\text{reg}}$ the estimator which minimizes the function $U({{\theta}})$, and is given by:

[TABLE]

where $X_{c}\in^{n_{c}\times d}$ represents the set of covariate vectors of the clean data points and $y_{c}\in^{n_{c}}$ represents the corresponding response values. In addition, we define the following functions required for the analysis of the robust gradient estimator,

[TABLE]

We make the following moment assumptions on the covariates in the clean data set $\mathcal{D}_{c}$.

• **Assumption 3 **(Positive-definite data covariance). The unnormalized data covariance matrix is positive definite: $\widetilde{\Sigma}_{x}\succ 0$.

• **Assumption 4 **(Bounded fourth moment). The data satisfies a bounded fourth moment condition, i.e., for every unit vector $v$, we have that

[TABLE]

for come constant $C_{x,4}$.

Note that these assumptions are satisfied with high probability if say, each $x_{i}\in\mathcal{D}_{c}$ is sampled i.i.d. from the standard normal distribution. The following corollary then instantiates the guarantees of Theorem 2 for the specific function $f({{\theta}})$ defined for the RBLR problem in Equation (16).

Corollary 5

*Consider the RBLR problem described above with posterior given by Equation (14) and data sampled from the Huber model. Then the iterates of Rob-ULA with step size $\eta\leq\frac{1}{n\bar{L}}$ and $h=n\eta$ satisfy:

---

[TABLE]

*with probability at least $1-\delta$ for $\epsilon\leq\frac{\bar{m}^{2}}{4CC_{\small{\Sigma,1}}\log d}$. The constant $C$ depends on $C_{x,4}$ from Assumption 5.2 and the remaining parameters are defined as:

---

[TABLE]

Remark.

As in the RBME problem, the quantities $h$, $\bar{L}$ and $\bar{m}$ are asymptotically independent of the sample size $n$. However, the guarantees above hold only for a value of $\epsilon\leq\tilde{\mathcal{O}}(\frac{1}{\kappa(\Sigma_{x})})$, that is, they depend on the condition number of the covariate distribution. Such a dependence seems inherent to the problem of linear regression since the adversary is allowed to corrupt the covariate vector arbitrarily.

6 Experiments

In this section, we compare the performance of the proposed Rob-ULA with the non-robust variant ULA. We first compare them on synthetic data sets for the problem of Bayesian mean estimation and Bayesian linear regression in order to understand the variation in performance as a function of the problem parameters. In Section 6.2, we perform experiments comparing the algorithms on some real-world binary classification data sets using logistic regression.

6.1 Synthetic data sets

6.1.1 Robust Bayesian Mean Estimation

In this section, we focus on experiments related to the robust Bayesian mean estimation problem described in Section 5.1. We begin by describing the experimental setup before proceeding to a discussion of the experimental findings.

Experiment setup.

The mean vector $\theta\in^{d}$ was sampled as a uniform distribution independently in each coordinate over the interval $[0,1]$. The clean samples $z_{i}$ were then obtained independently from $\mathcal{N}(\theta,I)$. The corrupted distribution $Q$ was chosen as the Gaussian distribution $\mathcal{N}(\theta_{c},I)$ with mean given by $\theta_{c}=\theta+\theta_{\sf cor}$ where each entry of $\theta_{\sf cor}$ is sampled i.i.d. from the uniform distribution over $[0,10]$. The default parameters were set as follows: number of clean samples $n_{c}=1000$, dimension $d=200$ and fraction of corruption $\epsilon=0.2$. In every experiment, one of the parameters was varied keeping the others fixed. For both Rob-ULA and ULA, the burn-in period was set to $300$ samples and a total of $n_{\sf samp}=1000$ samples were collected following the burn-in period. Each experiment was repeated for 10 runs and we report the mean performance of the methods across these runs.

Recovery guarantees.

Figures 1 and 2 compare the performance of the algorithms for the mean estimation problem. Figure 1 shows the variation in parameter recovery error, $\|\hat{\theta}-\theta\|_{2}$ where $\hat{\theta}=\frac{1}{n_{\sf samp}}\sum_{i}\theta_{i}$ is the average of the collected samples. Figure 1(a) shows that with increasing number of data points, the error in Rob-ULA’s estimate decreases until it starts to saturate. In addition, with an increasing dimension and fraction of outliers, the error in estimation for both Rob-ULA and ULA increases. This is consistent with Theorem 2. Figure 2 studies the variation in average log-likelihood of the average estimate $\hat{\theta}$ on a held-out test set. We also plot the likelihood values obtained by plugging in the true parameter (dotted line). The samples output by Rob-ULA have likelihood values identical to the true underlying parameter while those of ULA are much lower for all experiments. Note that ULA fails to have finite likelihood values (up to Matlab precision) for dimensions $d>200$ and fraction of corruption $\epsilon>0.2$ and hence have been omitted from the curves.

6.1.2 Robust Bayesian Linear Regression

In this section, we discuss the robust Bayesian linear regression problem described in Section 5.2.

Experiment setup.

The true parameter $\theta^{*}$ was selected similarly to the mean vector in the RBME experiments. The $x$ were sampled i.i.d. from $\mathcal{N}(0,I)$ and the corresponding response variable were set as $y=x^{\top}\theta^{*}+\mathcal{N}(0,1)$. For the corrupted distribution, each coordinate of the feature vector $x$ was sampled i.i.d. from a $\chi^{2}$ distribution and the corresponding response variables were set as $y=x^{\top}\theta^{*}+\text{Unif}[0,10]$. The default parameters were set as follows: number of clean samples $n_{c}=1000$, dimension $d=200$ and fraction of corruption $\epsilon=0.2$. In every experiment, one of the parameters was varied keeping the others fixed. For both Rob-ULA and ULA, the burn-in period was set to $100$ samples and a total of $n_{\sf samp}=300$ samples were collected. Each experiment was repeated for 10 runs and we report the mean performance of the methods across these runs.

Performance guarantees.

Figure 3 shows the performance of Rob-ULA and ULA on the linear regression problem in terms of parameter recovery. Note that similar to the mean estimation setup, the error curves for Rob-ULA are lower than those for ULA as we vary number of data points, the dimensionality of the problem and the fraction of corruptions. The error shows a decreasing trend with increasing number of samples but increases with increasing dimension and fraction of outliers, showing that the robust problems indeed become harder in higher-dimensional spaces and with a larger fraction of outliers.

6.2 Real-world data sets

In this section, we explore the performance of Rob-ULA on several real-world binary classification data sets obtained from the UCI repository [9]. We use a logistic regression model. While technically this model does not fall within the scope of Theorem 2 (primarily because of the strong-convexity assumption), we find that the experimental results are nonetheless consistent with the theory. We begin by providing details on the data sets used and then proceed to the experimental observations. For all the experiments in the section, the standard normal distribution was chosen as the prior.

Data sets.

The logistic regression experiments were carried out with the following publicly available binary classification data sets: a) Astro [14], b) Phishing [26], c) Breast-Cancer [25], d) Diabetes [9] and e) German Credit [9]. We normalized the features to scale between $[-1,1]$ and used $70\%$ of the available data for training purposes and the remaining $30\%$ for testing purposes. For cases where we were required to tune hyperparameters, we used $20\%$ of the train data for validation purposes. Once the hyperparameters were fixed, we retrained the model with the complete training set. In order to understand the effects of corruptions in these data sets, we manually add corruptions to the training subset in the form of label flips (for randomly chosen data points) for the Breast-Cancer, Diabetes and German Credit data set. The experiments with Astro and Phishing data sets are with the original uncorrupted data sets.

Evaluation Metric.

For the above data sets, since the true underlying logistic parameter is unknown, we evaluate the algorithms using the log-likelihoods on the test set. For all data sets, we show two plots: plot (a) displays log-likelihood per data point in the test set, sorted in descending order for both algorithms and plot (b) shows a histogram of log-likelihoods. Plot (a) helps understand trends in prediction likelihoods by showing how the prediction quality degrades while Plot (b) provides an understanding of how the likelihoods concentrate.

6.2.1 Binary Classification

Figure 4 compares the performance of Rob-ULA on a binary classification task with no corruptions in the training set for the Astro and Phishing data sets. In both figures, Rob-ULA is seen to perform better than vanilla ULA: the histogram of likelihoods is more concentrated towards the origin. These data sets are not linearly separable and hence the logistic model may provide a poor fit to the data; Rob-ULA exploits this fact and focuses on a subset of points which it can fit well. This allows it to perform better for a larger range of points as compared to vanilla ULA. These experiments show that if there is model misspecification and the chosen model doesn’t fit the complete data, the robust model might fit data selectively and give better confidence bounds for those data points.

6.2.2 Binary Classification with Label Flips

Figure 5 shows the performance of Rob-ULA for the Breast-Cancer ($\epsilon=0.10$), German Credit ($\epsilon=0.10$) and Diabetes ($\epsilon=0.15$) data sets respectively. For these data sets, we manually added corruptions via label flips. For all three data sets, we see a similar trend in likelihood plots: Rob-ULA performs better than vanilla ULA for a majority of data points but its prediction quality decreases for the tail points. This can also be seen in the three histogram plots (part (b) of the respective figures) wherein Rob-ULA has likelihoods extending to larger negative values. This particular behavior can be attributed to the fact that Rob-ULA is unable to learn effective representations in the space where label flips were added since it chooses to ignore those data points. Hence, in the region corresponding to the uncorrupted points, Rob-ULA achieves higher likelihood values than vanilla ULA but for the corrupted regions, the performance of Rob-ULA degrades slightly. This behavior was consistently seen for varying levels of $\epsilon$ with minor shifts in the likelihood curves for different corruption levels.

7 Conclusions

We have discussed the problem of robustness to adversarial outliers in a Bayesian framework and proposed Rob-ULA, a robust extension of the classical Unadjusted Langevin algorithm. We obtain nonasymptotic convergence guarantees for Rob-ULA.

We identify multiple directions for future work. On the statistical side, it would be interesting to extend the robustness guarantees of Rob-ULA for statistical models which do not fall within the scope of our current assumptions, notably the case of nonconvex likelihood functions. On the computational side, an important question to understand is whether one can accelerate the convergence of Rob-ULA in the presence of outliers.

Acknowledgments

This work was done in part while KB, YM and PLB were visiting the Simons Institute for the Theory of Computing.

Appendix A Proof of Lemma 1

We begin by obtaining a bound on the performance of Algorithm 2 in the one-dimensional setting and use this as a building block towards the proof for the general $d$-dimensional setting. For ease of exposition, we denote by $\mu_{\theta}$ and $\hat{\mu}$ the mean gradient $\nabla U_{\theta}$ and its estimate $\widehat{\nabla}U_{\theta}$ respectively.

A.1 Proof for 1 dimensional setting

We begin by analyzing Algorithm 2 for the case when $d=1$.

Lemma 6

Let $P_{\theta}$ denote the empirical distribution on $\mathcal{D}_{c}$ in $\mathbb{R}$ with mean $\mu_{\theta}$, variance $\sigma$ and with fourth moment constant $C_{4}$. Let $\eta$ be the fraction of corruption in the samples in $\mathcal{D}$. Then Algorithm 2 returns an estimate of the mean $\mu$ such that, for a universal constant $C$,

[TABLE]

Proof A.7**.**

Let $I_{1-\eta}$ be the interval around the true sample mean $\mu_{\theta}$ containing $1-\eta$ fraction of $\mathcal{D}_{c}$. Using the bounded fourth moment assumption on $P_{\theta}$, we have that,

[TABLE]

Let $\tilde{S}$ be the set of smallest interval containing $(1-\eta)^{2}$ fraction of all the points and let this interval be denoted by $\tilde{I}$. Now, we have that $\text{length}(\tilde{I})\leq\text{length}(I_{1-\eta})$ since the chosen interval has the smallest length. Further, $\tilde{S}$ contains at least $1-3\eta$ fraction of $\mathcal{D}_{c}$.

For any value of $\eta<\nicefrac{{1}}{{6}}$, we have that the intervals $\tilde{I}$ and $I$ must overlap. Therefore, adversarially corrupted points in the set $\tilde{S}$ are within distance $2\cdot\text{length}(I_{1-\eta})$ of the true parameter $\mu_{\theta}$. We now bound the deviation of the estimate $\hat{\mu}$ from $\mu_{\theta}$ by controlling the sources of error:

a. Error due to points in $\tilde{S}$ from $\mathcal{D}_{a}$:

Since there can be at most $\eta$ fraction of corrupted points in our selected sample and each of them within distance $2\cdot\text{length}(I_{1-\eta})$, their total contribution to the deviation is upper bounded by $2\eta\cdot\text{length}(I_{1-\eta})$.

b. Error due to points in $\tilde{S}$ from $\mathcal{D}_{c}$:

Define $\mathcal{A}$ to be the event that a point of $\mathcal{D}_{c}$ is present in the set $\tilde{S}$. From our discussion above, we have that $P(A)>1-3\eta>1/2$. Using Lemma 3.11 **[20]**, we have that there exits a constant $C$ such that

[TABLE]

Combining the analysis of parts (a) and (b) above, we get that,

[TABLE]

*Plugging in the bound for $\text{length}(I_{1-\eta})$ completes the proof. *

A.2 Proof for d-dimensional setting

We now proceed to prove the robustness properties of Algorithm 2 for the general $d$-dimensional setting. Through the course of this section, we let $\tilde{S}$ be the set of points returned after the outlier truncation procedure with $\tilde{S}_{c}$ being the clean points and $\tilde{S}_{a}$ being the adversarially corrupted points. Also, we denote by $\mu_{\tilde{S}}:=\text{mean}(\tilde{S})$, $\mu_{\tilde{S}_{c}}:=\text{mean}(\tilde{S}_{c})$ and $\mu_{\tilde{S}_{a}}:=\text{mean}(\tilde{S}_{a})$ the corresponding mean vectors of the relevant subsets.

Lemma A.8**.**

Let $P_{\theta}$ be the empirical distribution over $\mathcal{D}_{c}$ in $\mathbb{R}^{d}$ with mean $\mu_{\theta}$, covariance matrix $\Sigma_{\theta}$ and with fourth moment constant $C_{4}$. Let $\eta$ be the fraction of corrupted data points in $\mathcal{D}$. Then, we can obtain a vector $a\in\mathbb{R}^{d}$ such that for a constant $C$,

[TABLE]

Proof A.9**.**

Let $e_{1},\ldots,e_{d}$ be the $d$ canonical basis vectors. Projecting the problem onto these vectors and solving in each direction independently using the method for one dimension, we obtain the bound above by using Lemma 6 separately in each dimension.

Lemma A.10**.**

Let $\eta$ denote the fraction of outliers in $\mathcal{D}$. After the outlier truncation procedure, for every point in the returned set $\tilde{S}$, we have that for a constant $C$,

[TABLE]

Proof A.11**.**

Let $B^{*}=\mathcal{B}(\mu_{\theta},r_{1}^{*})$ for $r_{1}^{*}=\frac{CC_{4}^{\frac{1}{4}}}{\eta^{\frac{1}{4}}}\sqrt{d\|\Sigma_{\theta}\|_{2}}$ be the $\ell_{2}$ ball of radius $r_{1}^{*}$ around $\mu_{\theta}$. In order to bound the fraction of clean points in $B^{*}$, observe that,

[TABLE]

*where the probability is with respect to the empirical distribution over $\mathcal{D}$. Now, $\mathbb{E}\left[(\right]\|x-\mu_{\theta}\|_{2}^{4}\leq d^{2}\max_{i}\mathbb{E}\left[(\right](x-\mu_{\theta})_{i}^{4})\leq C_{4}d^{2}\|\Sigma_{\theta}\|_{2}^{2}$. Plugging this value in Equation 18, we get that at least $1-\eta$ fraction of the points of $\mathcal{D}_{c}$ lie in $B^{*}$.

Using Lemma A.8, we have that for $r_{2}^{*}=CC_{4}^{\frac{1}{4}}\sqrt{\text{tr}(\Sigma_{\theta})}\eta^{\frac{3}{4}}$, there are at least $(1-\eta)$ fraction of good points at a distance $r_{1}^{*}+r_{2}^{*}$ away from $a$. For $\eta<\nicefrac{{1}}{{6}}$, we have that the chosen ball of points around $\tilde{S}$ and $B^{*}$ must overlap. Therefore, minimum radius ball has radius at most $r_{1}^{*}+r_{2}^{*}$ and when combined with the bound from Lemma A.8 and triangle inequality, we get that,

[TABLE]

which completes the proof.

Lemma A.12**.**

Let $\eta$ be the fraction of corrupted points in $\mathcal{D}$. Then we have that after the outlier truncation step of Algorithm 2, for a constant $C$,

[TABLE]

Proof A.13**.**

We first consider the bound on the mean shift and then proceed with the bound on the covariance matrix.

Mean Shift Bound:

Let $\mathcal{A}$ be the event that a point $x\in\mathcal{D}_{c}$ is not removed by the outlier truncation procedure. Then using Lemma 3.11 [**[20]**] for $\eta<\nicefrac{{1}}{{6}}$ for the random variable $X=x^{\top}\frac{\mu_{\tilde{S}_{c}}-\mu_{\theta}}{\|\mu_{\tilde{S}_{c}}-\mu_{\theta}\|_{2}}$ for $x\sim\mathcal{D}_{c}$, we have that,

[TABLE]

Covariance Matrix Bound:

Consider the following decomposition for bounding the spectral norm of $\Sigma_{\tilde{S}_{c}}$:

[TABLE]

where $\zeta_{1}$ follows by using Corollary 3.13 **[20*]** for the same event $\mathcal{A}$ as above in the mean shift bound. *

Lemma A.14**.**

$P_{W}$* is the projection operator on the bottom $d/2$ eigenvectors of the matrix $\Sigma_{\tilde{S}}$. Then, for a constant $C$, we have that,*

[TABLE]

where $\delta_{\mu}:=\mu_{\tilde{S}_{a}}-\mu_{\tilde{S}_{c}}$.

Proof A.15**.**

Consider the matrix $\Sigma_{\tilde{S}}$. It can be decomposed as:

[TABLE]

By Weyl’s inequality, we have that,

[TABLE]

We begin by first controlling the term $\lambda_{d/2}(\Sigma_{2})$. We have that,

[TABLE]

where $\zeta_{1}$ follows by using the fact that all selected points are in a ball of radius $r_{1}^{*}+r_{2}^{*}$ where $r_{1}^{*}$ and $r_{2}^{*}$ are as defined in Lemma A.10. Next we consider the term $\lambda_{1}(\Sigma_{1})$ as follows,

[TABLE]

where $\zeta_{1}$ follows from using the bound in Lemma A.12. Combining Equations (19) and (20), we have that,

[TABLE]

Using the fact that $P_{W}$ is the projection operator on the bottom $d/2$ eigenvectors of the matrix $\Sigma_{\tilde{S}}$, we have that,

[TABLE]

Following some algebraic manipulation as in [20], we obtain that

[TABLE]

which completes the proof.

Proof of Lemma 1

Let $\tilde{S}$ be the subset of samples returned by the outlier truncation procedure and let $\tilde{S}_{c}$ be the set of clean points contained in $\tilde{S}$. Then, we have,

[TABLE]

where $\zeta_{1}$ follows from the orthogonality of the spaces $V$ and $W$, $\zeta_{2}$ follows from using triangle inequality and $\zeta_{3}$ follows from contraction of projection operators. Note that $(I)$ is a problem defined on the subspace $V$ which is of ambient dimension $d/2$ and is solved recursively by Algorithm 2. Thus, one can recursively bound the overall error of the algorithm as,

[TABLE]

Using Lemma A.12 and Lemma A.14, we can bound the above error as,

[TABLE]

which completes the proof of the lemma. ∎

Appendix B Convergence of

Algorithm 1**.**

: Proofs for Auxiliary Lemmas

Lemma B.16**.**

For $\Theta_{t}$ following Eq. (7), if the initial iterate $\Theta_{0}\sim\mathcal{N}\left(0,\dfrac{1}{n\bar{L}}\mathrm{I}\right)$, the fraction of corruption $\epsilon\leq\dfrac{\bar{m}^{2}}{4C_{R}C_{\small{\Sigma,1}}\log d}$, and scaled step-size $h=n\eta\leq\dfrac{1}{\bar{L}}$, then for all $k\in\mathbb{N}^{+}$,

[TABLE]

Proof B.17**.**

Consider first the initial iterate for $k=0$. The distribution $p_{0}$ satisfy $\mathbb{E}_{\theta\sim p_{0}}\left[\left\|\theta\right\|_{2}^{2}\right]=\dfrac{d}{n\bar{L}}\leq\dfrac{4}{\bar{m}^{2}}C_{14}\epsilon\log d+\dfrac{4d}{n\bar{m}}$. We will prove the lemma statement by strong induction.

In the induction hypothesis step, assume that for some $k\geq 0$, for all $t=0,\eta,\ldots,k\eta$, $\mathbb{E}_{\theta\sim p_{t}}\left[\left\|\theta\right\|_{2}^{2}\right]\leq\dfrac{4}{\bar{m}^{2}}C_{R}C_{\small{\Sigma,2}}\epsilon\log d+\dfrac{4d}{n\bar{m}}$. We consider obtaining a bound on $\mathbb{E}_{\theta\sim p_{(k+1)\eta}}\left[\left\|\theta\right\|_{2}^{2}\right]$, where $p_{t}$ follows Equation 7, for $t\in(k\eta,(k+1)\eta]$ (denote $\tau=t-k\eta\in(0,h]$):

[TABLE]

Given the above equation, we consider the bound on $\mathbb{E}\left[\left\|\Theta_{t}-\tilde{\theta}\right\|_{2}^{2}\right]$ for some $t\in(k\eta,(k+1)\eta]$ as follows:

[TABLE]

We next define $\nu=n\tau$ and obtain a bound on the term $\mathbb{E}\left[\left\|\left(\Theta_{k\eta}-\tilde{\theta}\right)-\widehat{\nabla}f(\Theta_{k\eta})\tau\right\|_{2}^{2}\right]$:

[TABLE]

where ${\sf(i)}$ follows by an application of Cauchy–Schwarz inequality. Next, using the assumption on the robust estimation of the gradient from Theorem 2, we have that

[TABLE]

and further simplifying the above using Lemma B.18, we obtain

[TABLE]

Therefore, since $\nu\leq\dfrac{1}{\bar{L}}$ and the corruption factor$\epsilon\leq\dfrac{\bar{m}^{2}}{4C_{R}C_{\small{\Sigma,1}}\log d}$, we have that

[TABLE]

where we have defined the constants $C_{13}:=C_{R}C_{\small{\Sigma,1}}$ and $C_{14}:=C_{R}C_{\small{\Sigma,2}}$. Using the above bounds, we have that

[TABLE]

Note that $\nu\leq\dfrac{1}{\bar{L}}\leq\dfrac{1}{\bar{m}}$ and $\mathbb{E}\left[\left\|\Theta_{k\eta}-\tilde{\theta}\right\|^{2}\right]\leq\dfrac{4C_{14}}{\bar{m}^{2}}\epsilon\log d+\dfrac{4d}{n\bar{m}}=\dfrac{2}{\bar{m}}\left(\dfrac{2}{\bar{m}}C_{14}\epsilon\log d+\dfrac{2}{n}d\right)$. Combining these, we can obtain the final bound stated in the lemma as follows:

[TABLE]

for any $t\in(k\eta,(k+1)\eta]$. This concludes the proof.

We now prove the Lemma B.18 which was used in the proof of Lemma B.16.

Lemma B.18**.**

For $\tau\leq\dfrac{1}{\bar{L}}$, and $\dfrac{1}{n}\nabla f(\theta)$ being $\bar{m}$ strongly convex and $\bar{L}$ Lipschitz smooth, we have

[TABLE]

Proof B.19**.**

To bound $\left\|\theta-\dfrac{1}{n}\nabla f(\theta)\nu\right\|_{2}^{2}$, we consider the following function: $F(\theta)=\dfrac{1}{2}\left\|\theta\right\|_{2}^{2}-\dfrac{1}{n}f(\theta)\nu$. First note strong convexity and Lipschitz smoothness of $\dfrac{1}{n}\nabla f(\theta)$ implies that $\bar{m}\mathrm{I}\preceq\dfrac{1}{n}\nabla^{2}f(\theta)\preceq\bar{L}\mathrm{I}$. Thus with $\nu\leq\dfrac{1}{\bar{L}}$, we have that

[TABLE]

Note that the point $\tilde{\theta}$ satisfies $\nabla f(\tilde{\theta})=0$. Using this we have that:

[TABLE]

where ${\sf(i)}$ follows from the Lipschitz-smoothness of $F$.

Lemma B.20** (Bound on Initial Error).**

If we let the initial iterate $\Theta_{0}$ have distribution given by

[TABLE]

*and $p^{*}$ following Assumptions 4.2–4.2, then the initial error $\text{KL}(p_{0}\;\|\;p^{*})$ is bounded as:

---

[TABLE]

Proof B.21**.**

We want to bound ${\text{KL}(p_{0}\;\|\;p^{*})}=\displaystyle\int p_{0}(\theta)\ln\left(\dfrac{p_{0}(\theta)}{p^{*}(\theta)}\right)\mathrm{d}\theta$, where $p^{*}(\theta)\propto e^{-f(\theta)}$. First define $\bar{f}(\theta)=f(\theta)-f(\theta^{*})$, where $\theta^{*}$ is the minimum of $f$. Then

[TABLE]

By Assumptions 4.2 and 4.2, we have that $\dfrac{m}{2}\|\theta\|^{2}\leq\bar{f}(\theta)\leq\dfrac{L}{2}\|\theta\|^{2}$, $\forall\theta\in\mathbb{R}^{d}$. Therefore,

[TABLE]

where ${\sf(i)}$ follows from using the Lipschitz-smoothness of $\bar{f}$. Hence

[TABLE]

We can also calculate similarly that

[TABLE]

Combining the above, we get that

[TABLE]

This concludes the proof of the statement.

Appendix C Proofs for Mean Estimation and Regression

C.1 Proof of Corollary 4

Throughout this proof, we condition on the high probability event described by Equation (5). We proceed to obtain a bound on the strong-convexity parameter $m$ and Lipschitz smoothness parameter $L$ for the function $f({{\theta}})$ defined above in Equation (13). The gradient $\nabla f({{\theta}})$ is given by

[TABLE]

and the corresponding hessian $\nabla^{2}f({{\theta}})$ is given by,

[TABLE]

Since both the matrices $\Sigma_{0}$ and $\Sigma$ are positive-definite, we have the following bounds for the parameters $m$ and $L$:

[TABLE]

We now obtain a bound on the covariance $\Sigma_{\theta}$ of the gradients $\nabla g_{i}(\theta)$ using the empirical distribution of the points in $\mathcal{D}_{c}$ as follows:

[TABLE]

where $z\sim_{u}\mathcal{D}_{c}$ denotes data sampled uniformly from the data set $\mathcal{D}_{c}$. Thus, from the above equation, we get that $C_{\small{\Sigma,1}}=0$ and $C_{\small{\Sigma,2}}=\frac{\lambda_{\text{max}}(\Sigma_{z})}{\lambda_{\text{min}}(\Sigma)}$. Plugging in these values in Theorem 2 gives us the desired bound. ∎

Appendix D Proof of Corollary 5

We condition on the high probability event described by Equation (15). We begin by obtaining a bound on the strong-convexity parameter $m$ and Lipschitz smoothness parameter $L$ for the function $f({{\theta}})$ defined above in Equation (16). The gradient $\nabla f({{\theta}})$ is given by

[TABLE]

and the corresponding hessian $\nabla^{2}f({{\theta}})$ is given by,

[TABLE]

Since both the matrices $\Sigma_{0}$ and $\widetilde{\Sigma}_{x}$ are positive-definite, we have the following bounds for the parameters $m$ and $L$:

[TABLE]

We now proceed to obtain the bound on the spectral norm of the covariance matrix $\Sigma_{{\theta}}$ of the gradients $\nabla g_{i}(\theta)=x_{i}(\langle x_{i},{{\theta}}\rangle-y_{i})$ using the empirical distribution of the points in $\mathcal{D}_{c}$. We use the notation $\mathbb{E}_{\mathcal{D}_{c}}$ to denote $\mathbb{E}_{(x,y)\sim_{u}\mathcal{D}_{c}}$, the sampling of pairs $({\bf x},y)$ uniformly from the clean data set $\mathcal{D}_{c}$.

[TABLE]

where ${\sf(i)}$ follows from the fact that $\mathbb{E}_{\mathcal{D}_{c}}\left[\nabla g_{i}(\theta)\right]\mathbb{E}_{\mathcal{D}_{c}}\left[\nabla g_{i}(\theta)\right]^{\top}\succeq 0$ and ${\sf(ii)}$ follows from the Cauchy-Schwarz inequality. We now obtain a bound on the two expectations in the above equation.

Bound on $\mathbb{E}_{\mathcal{D}_{c}}\left[(v^{\top}x)^{4}\right]$: This term can be bounded using bounded fourth moment assumption (Assumption 5.2) as follows:

[TABLE]

Bound on $\mathbb{E}_{\mathcal{D}_{c}}\vphantom{\sum}\left[(\langle x,{{\theta}}\rangle-y)^{4}\right]$: We simplify this term by using the modelling assumption on the data $(x,y)$ and then proceed to bound this using the $c_{r}$-inequality.

[TABLE]

where $\Delta_{{\theta}}:={{\theta}}-{{\theta}}^{*}$, ${\sf(i)}$ follows from using the $c_{r}$-inequality $\mathbb{E}|X+Y|^{r}\leq 2^{r-1}\left(\mathbb{E}|X|^{r}+\mathbb{E}|Y|^{r}\right)$ and ${\sf(ii)}$ follows from Assumption 2 on bounded fourth moment of the covariates. Using Lemma D.22 for bounding the fourth moment of the noise variables, we have with probability at least $1-\delta$,

[TABLE]

Substituting the bounds obtained in Equations (29) and (30) in Equation (D), along with an application of triangle inequality, we have that with probability at least $1-\delta$,

[TABLE]

One can now use the above bounds in conjunction with the values of $\bar{L}$ and $\bar{m}$ from Equation (27) to obtain the final result. ∎

The following lemma obtains a concentration bound for the fourth moment of a Guassian random variable and can be obtained by appropriate instantiation of the Hypercontractivity Concentration Inequality (Theorem 1.9) by Schudy and Sviridenko [35].

Lemma D.22** (Concentration Bound for Gaussian Fourth Moment).**

Let $z_{1},z_{2},\ldots,z_{n}$ be i.i.d. random variable sampled from $\mathcal{N}(0,1)$. Then, there exists a universal constant $C_{z,r}$ such that for any $\epsilon>0$, we have that,

[TABLE]