On Empirical Risk Minimization with Dependent and Heavy-Tailed Data
Year: 2022

Abstract:
In this work, we establish risk bounds for the Empirical Risk Minimization (ERM) with both dependent and heavy-tailed data-generating processes. We do so by extending the seminal works [Men15, Men18] on the analysis of ERM with heavytailed but independent and identically distributed observations, to the strictly stationary exponentially β-mixing case. Our analysis is based on explicitly controlling the multiplier process arising from the interaction between the noise and the function evaluations on inputs. It allows for the interaction to be even polynomially heavy-tailed, which covers a significantly large class of heavy-tailed models beyond what is analyzed in the learning theory literature. We illustrate our results by deriving rates of convergence for the high-dimensional linear regression problem with dependent and heavy-tailed data.

heading: 1 Introduction
text: Given a random vector (X,Y ) ∈ Rd × R, with joint distribution (X,Y ) ∼ π, and a class of closed, convex set of functions F ⊂ L2(π), the objective in statistical learning theory is to find the best function in the set F that maps the input X to the target Y . The quality of this mapping is measured by a user-defined loss function ℓ : R → R+ ∪ {0}. The most well-studied approach for the above task is that of risk minimization, where the best function is defined as the one that minimizes the expected loss over the set F :
f∗ = argmin f∈F Pℓf := argmin f∈F Eπ [ℓ (f(X)− Y )] .
The above problem requires the knowledge of the distribution π which is typically unknown in practice. However, we are usually given observations Zi = (Xi, Yi) for i = 1, . . . , N , from the distribution π which leads to the Empirical Risk Minimization (ERM) procedure defined as
f̂ = argmin f∈F PN ℓf := argmin f∈F
1
N N∑ i=1 ℓ (f(Xi)− Yi) .
The convergence of the empirical risk minimizer f̂ to the true risk minimizer f∗ is typically analyzed by considering the underlying empirical process, a topic which dates back to the seminal work of [VC71]; see also [VDVW96, vdG00, BBM05, Kol06, Kol11]. In a representative analysis in this
∗Department of Statistics, University of California, Davis.abroy@ucdavis.edu. Research of this author was supported in part by NSF TRIPODS grant CCF-1934568
†Department of Statistics, University of California, Davis. kbala@ucdavis.edu. Research of this author was supported in part by UC Davis CeDAR (Center for Data Science and Artificial Intelligence Research) Innovative Data Science Seed Funding Program.
‡Department of Computer Science and Department of Statistical Sciences at the University of Toronto, and Vector Institute. erdogdu@cs.toronto.edu. Research of this author was supported in part by NSERC Grant [2019-06167], Connaught New Researcher Award, CIFAR AI Chairs program, and CIFAR AI Catalyst grant
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
setting, a majority of the works assume the observations Zi are generated independent and identically distributed (iid) from π, and the analysis is based on uniform concentration. However, there are important limitations associated with this approach, particularly due to the (Talagrand’s) contraction principle which naturally requires a Lipschitz loss function (see, for example, [LT13, Corollary 3.17] or [Kol11, Theorem 2.3]). As a result, in order to work with standard (unbounded) loss functions such as squared-error loss or Huber loss, it is generally assumed that the range of f ∈ F is uniformly bounded and/or the noise ξ := Y − f(X) is also uniformly bounded π-almost surely. Several attempts have been made in the literature to overcome the limitations of the standard ERM analysis. A significant progress was made by Mendelson [Men15, Men18], who proposed the socalled learning without concentration framework for analyzing ERM procedures with unbounded noise or loss functions. The approach is based on a combination of small-ball type assumption on the input samples Xi, along with developing multiplier empirical process inequalities under weaker moment assumptions. We refer the interested reader, for example, to [Men17b, Men17a, LM18, LRS15, GM20] for details. The aforementioned works, while relaxing the prior analysis of ERM to handle heavy-tailed data-generating process (DGP), still require the more stringent iid assumption for their analysis. This restricts the practical applicability of the developed theoretical results significantly. Indeed, heavy-tailed and dependent data appear naturally in various practical learning scenarios [BF89, JM01, DKBR07]; however, theoretical guarantees are still missing.
Our Contributions: Aiming to fill the above gap, we analyze ERM with convex loss functions (that are locally strongly-convex around the origin) when the DGP is both heavy-tailed and non-iid. We do so by extending the small-ball technique of [Men15, Men18] to the strictly stationary exponentially β-mixing data. In the iid case, the interaction between the noise and the inputs is handled by an analysis based on multiplier empirical process. However, developing similar techniques in the noniid case is fundamentally restrictive due to the limitations of the analysis based on empirical process. We side-step this issue for the non-iid case by directly making assumptions on the interaction, which allows for it to be either exponentially or polynomially heavy-tailed. For the exponentially heavy-tailed interactions, we leverage the concentration inequalities developed by [MPR11]. For the polynomially heavy-tailed case, we develop new concentration inequalities extending the recent work [BMdlP20] to β-mixing random variables. We illustrate our results in the context of ERM with sparse linear function class and stationary β-mixing DGP under both squared and Huber loss.
Motivation: A natural question arises in this context: Why study ERM with convex loss functions when the DGP is heavy-tailed? Firstly, convex loss functions cover a large class of robust loss function that are tailored to deal with the heavy-tailed behavior present in the noise and/or input data. Some examples include the Huber loss [Hub92], conditional value-at-risk [RU02, RS06, MGW20, SY20] and the so-called spectral risk measures [Ace02, HH21]. While there exist studies for nonconvex loss functions suited for heavy-tailed input data (for example, [Loh17]), such analyses are mostly in a model-based setting and focus on estimation error. Secondly, while alternatives to ERM have also been proposed and analyzed in the literature for the iid case (with the most prominent one being the median-of-means framework and its variants [MM19, LM19, LL20, BM21]), it is not immediately clear how to extend such methods to the dependent DGP that we consider in this paper. We view our work as taking the first step in developing risk bounds for statistical learning when the DGP is both heavy-tailed and dependent.
Related Works: The seminal work [Yu94] extended the analysis based on empirical process to the stationary mixing process using a blocking technique. [Irl97] and [BR97] studied consistency of non-parametric regression methods under mixing and exchangeability conditions on the DGP, respectively. [Nob99] established lower bounds to achieving consistency when learning from dependent data. [SHS09] studied consistency of ERM with Zi being an α-mixing (not necessarily stationary) process, when F is a reproducing kernel Hilbert space. More recently, [Han21] and [DT20] studied learnability under a general stochastic process setup. The works [AV90, BL97, Pes10, Gam03] extend Valiant’s Probably Approximately Correct (PAC) learning model to Markovian and related drifting DGP, assuming bounded loss function and/or noise to obtain rates of convergence. Furthermore, [ZCY12] and [HS14] analyzed ERM for least-squares regression (with bounded noise) with clipped loss functions and an α-mixing DGP. Rademacher complexity results for predominantly stationary dependent processes were developed by [MS11] and [MR08]. [RSS10] and [ALW13] developed PACBayes bounds in the non-iid setting. [RST15, RS14] developed notions of sequential Rademacher complexity to characterize the complexity in online nonparametric learning in the worst-case. More recently [DDDJ19, KDD+21] considered learning under weakly-dependent data for specific models.
However, such works mainly rely on bounded loss functions in their analysis. Furthermore, [Mei00] and [AW12] studied model selection for time series forecasting in a possibly unbounded setup but do not consider conditional prediction and is limited to light-tailed cases.
Apart from the aforementioned works, the recent works [KM17, HW19, WLT20] are closely related to our setup as they consider rates of convergence of ERM under heavy-tailed and dependent DGP. In [KM17, Section 8], generalization bounds are developed when Zi is an asymptotically stationary β-mixing sequence. However, their conditions on the function class F are rather opaque and it is not clear if their method actually handles the heavy-tailed DGP that we focus on. [HW19] considered a setup based on a statistical model: for i = 1, . . . , N , Yi = f∗(Xi)+ϵi, with the following conditions: (i) ϵi being independent of Xi, (ii) Xi is iid, and (iii) ϵi is iid. For this setting, they assumed that the noise has a bounded p-th moment (with p ≥ 2) and F satisfies the standard entropy condition (see, for example [Kol06, Example 4]) with exponent α ∈ (0, 2) and obtained convergence rates of the order O(N− 1 2+α +N− 1 2+ 1 2p ). Finally, [WLT20] provides an analysis of L1-regularized ERM with quadratic loss and linear function class where the heavy-tailed behavior is induced by a sub-Weibull assumption, and data dependency is characterized by a stationary β-mixing condition. However, their analysis specializes to sparse linear function classes and their focus is on parameter estimation error and in-sample prediction accuracy.

heading: 2 Assumptions and Preliminaries
text: We assume that there exists a function f∗ ∈ F that minimizes the population risk E[ℓ(f(X)− Y )]. In what follows, we provide the conditions that we require on the DGP, specifically, the exponentially β-mixing condition, for characterizing the dependency among data points.
Definition 1 ([Yu94]). Suppose that {Zi}∞i=−∞ is a strictly stationary sequence of random variables. For any i, j ∈ Z ∪ {−∞,∞}, let σji denote the σ-algebra generated by {Zb} j b=i. Then for any positive integer b, the β-mixing coefficient of the stochastic process {Zi}∞i=−∞ is defined as
β(b) = sup n E B∈σn−∞
[ sup
A∈σ∞n+b |P(A |B)− P(A)|
] .
The sequence {Zi}∞i=−∞ is said to be β-mixing if β(b) → 0 as b→ ∞. Furthermore, it is said to be exponentially β-mixing if there exist β0, β1, r > 0 such that β(b) ≤ β0exp(−β1br) for all b.
The β-mixing condition is frequently used when studying non-iid DGP, and imposes a dependence structure between data samples that weakens over time. The coefficient β(b) is a measure of the dependence between events that occur within b units in time. Indeed, β-mixing is often used in the analysis of non-iid data in statistics and machine learning, [Vid13]. Before stating our assumptions formally, we present the following decomposition of empirical risk for a convex loss using Taylor’s expansion:
PNℓf ≥ 1
16N N∑ i=1 ℓ′′(ξ̃i)(f − f∗)2(Xi) + 1 N N∑ i=1 ℓ′(ξi)(f − f∗)(Xi), (1)
where ξ̃i is a suitably chosen midpoint between f(Xi)− Yi and f∗(Xi)− Yi := ξi. For quadratic loss functions ℓ(t) = t2, ℓ′(ξi) = 2ξi and ℓ′′(ξ̃i) = 2, ∀i. At a high level, establishing risk bounds boils down to proving a positive lower bound on the second term on the Right Hand Side (RHS) of (1), and a concentration result for the first term on the RHS with high probability. We now introduce the precise assumptions we make on the DGP and the function class to formalize the above strategy. For the sake of clearer exposition, we first introduce our assumptions in the context of quadratic loss and then indicate the changes required to handle more general locally strongly-convex loss functions.
Assumption 2.1 (Squared loss). The DGP {Zi}∞i=−∞ and the function class F satisfy the following:
(a) β-mixing data. The process {Zi}∞i=−∞ is a strictly stationary exponentially β-mixing sequence, i.e., β(k) ≤ exp(−ckη1), for some c, η1 > 0, with strict stationary distribution π.
(b) Small ball condition. Let F ⊂ L2(π) be closed, convex class of functions and define F − F := {f − h : f, h ∈ F}. Then, the function class F is such that there is a τ > 0 for which QF−F (2τ) > 0, where QH(u) = infh∈H P (|h| ≥ u∥h∥L2) .
(c) Deviations of interaction. The stationary noise ξ1 and the error (f − f∗)(X1) satisfy either:
(i) For all f ∈ F , and Z1 ∼ π, for some η2 > 0 we have
P (|ξ1(f − f∗)(X1)− E [ξ1(f − f∗)(X1)] | ≥ t) ≤ exp(1− tη2). (2)
or, (ii) For all f ∈ F , and Z1 ∼ π, for some η2 > 2, and c > 0 we have
P (|ξ1(f − f∗)(X1)− E [ξ1(f − f∗)(X1)] | ≥ t) ≤ ct−η2 . (3)
(d) Heavy-tail/data-mixing trade-off. Under condition (c)-(i), 1/η := 1/η1 + 1/η2 > 1.
Condition (a) above, exponentially β-mixing data, has been assumed in various works, see for example [Vid13, WLT20, KM17], to obtain rates of convergence for ERM procedures in general. Indeed, (exponential) mixing assumption holds in several time-series applications. For example, [Mok88] showed that certain ARMA processes can be modeled as an exponentially β-mixing stochastic process. Furthermore [VK06] showed that globally exponentially stable unforced dynamical systems subjected to finite-variance continuous density input noise give rise to exponentially mixing stochastic process; see also [FS12]. Condition (b), referred to as the well-known small-ball condition, has been previously employed in the iid case [Men15]. Intuitively, it models heavy-tailedness by restricting the mass allowed near any small neighborhoods of zero; thus, forcing the tails to be necessarily heavy. To our knowledge, the small-ball condition has not been used under dependent DGP assumptions. Condition (c) is proposed in this work as a way to model the interaction between the stationary noise ξ1 and the stationary error (f − f∗)(X1). For the iid setting, [Men15] modeled the interaction between ξ1 and (f − f∗)(X1) uniformly over the class of F via the multiplier empirical process and captured the complexity through a parameter αN (see (108) for the definition). This requires using symmetrization argument in the proof which is not applicable in the non-iid setting that we consider in this work. In addition, how different tail conditions on the data and noise affect the high-probability statement on the learning rate is not apparent from the parameter αN . We revisit the relationship between our condition (in the context of iid observations) and the multiplier empirical process approach used in [Men15] in Section F. Finally, condition (d) models the relationship between the allowed degree of dependency and the allowed degree of interaction between ξ1 and (f − f∗)(X1). Next, we modify condition (c) in Assumption 2.1 for the case of general locally strongly-convex loss functions because now the interaction part involves ℓ′(ξ) instead of ξ (recall the decomposition (1)).
Assumption 2.2 (Convex loss). When the loss function is locally strongly-convex around the origin and globally convex, condition (c) in Assumption 2.1 is modified as:
(c) Deviations of interaction. The stationary noise ξ1 and the error (f − f∗)(X1) satisfy either, (i) For all f ∈ F , and Z1 ∼ π, for some η2 > 0 we have
P (|ℓ′(ξ1)(f − f∗)(X1)− E [ℓ′(ξ1)(f − f∗)(X1)] | ≥ t) ≤ exp(1− tη2). (4)
or, (ii) For all f ∈ F , and Z1 ∼ π, for some η2 > 2, and c > 0 we have
P (|ℓ′(ξ1)(f − f∗)(X1)− E [ℓ′(ξ1)(f − f∗)(X1)] | ≥ t) ≤ ct−η2 . (5)
Complexity Measures: We now introduce the complexity measures that play a crucial role in characterizing the rates of convergence. The use of β-mixing assumption enables us to define complexity measures based on the blocking technique proposed by [Yu94], also utilized by the works of [MR08, KM17, WLT20]. We partition the training sample of size N , S := {Zi}Ni=1, into two sequences of blocks Sa and Sb. Each block in Sa, and Sb is of length a, and b respectively. Sa and Sb, both are of length µ, i.e., µ(a+ b) = N . Formally, Sa and Sb are given by
Sa = ( Z (a) 1 , Z (a) 2 , · · · , Z (a) µ ) with Z(a)i = {z(i−1)(a+b)+1, · · · , z(i−1)(a+b)+a},
Sb = ( Z (b) 1 , Z (b) 2 , · · · , Z (b) µ ) with Z(b)i = {z(i−1)(a+b)+a+1, · · · , z(i−1)(a+b)+a+b}. (6)
Based on this blocking technique, we require the following definition of Rademacher complexity.
Definition 2. Let {X̃i}µi=1 be an iid sample from the strict stationary distribution π. Let D be the unit-L2(π) ball centered at f∗. For every γ > 0, define
ωµ(F − F , γ) := inf { r > 0 : E [ sup
h∈(F−F) ∩ rD
∣∣∣∣∣ 1µ µ∑
i=1
ϵih(X̃i) ∣∣∣∣∣ ] ≤ γr } , (7)
where {ϵi}µi=1 are iid Rademacher variables taking values ±1 with probability 1/2.
The quantity ωµ(H, γ) provides a localized complexity measure for the function class F , and serves as a generalization of the standard Rademacher complexity in the non-iid setting. For the case of locally strongly-convex losses, we need the following related measures of complexity. Definition 3. For a function class H⊂L2(π), a sample of size N from a strictly stationary β-mixing sequence with stationary distribution π satisfying Assumption 2.1-(a), and ζ1, ζ2 > 0, we define:
ω1(H, N, ζ1) = inf { r > 0 : E [∥G∥H∩rD] ≤ ζ1rN η1 2(1+η1) } and ω2(H, µ, ζ2) = ωµ(H, ζ2),
where µ = N η1
(1+η1) , ∥G∥H = suph∈HGh, and {Gh : h ∈ H} is the canonical Gaussian process indexed by H with a covariance induced by L2(π). Moreover, we let
ωQ(F − F , N, ζ1, ζ2) := max(ω1(H, N, ζ1), ω2(H, µ, ζ2)).
Note that in Definitions 2 and 3, the scaling is in terms of number of blocks µ instead of N . The number of blocks µ can be thought of as the effective sample size under dependency, and as η1 → ∞, one has µ→ N . The term E [∥G∥H∩rD] appearing in Definition 3 is termed as the localized Gaussian width and is also a widely used complexity measure in the literature. Note that while in the case of quadratic loss function, the bound is in terms of local Rademacher-based complexity measure, whereas in the convex case, we require both Gaussian and Rademacher-based complexity measures to establish the bound, mostly due to technical reasons in the proof. An intuitive explanation for this has eluded us thus far; see also [Men18, Lemma 4.2].
Concentration Inequalities for Heavy-tails: We now restate [MPR11, Theorem 1], in a form adapted to our setting below. This result is required to handle interactions of noise and input that satisfy condition (c)-(i) of Assumption 2.1. It is straightforward to check that the conditions required by [MPR11] are immediately satisfied under our Assumption 2.1. Indeed, while the results in [MPR11] are stated for τ -mixing sequences, condition (a) in Assumption 2.1 implies that the process {Zi}∞i=−∞ is exponentially τ -mixing [CG14], i.e., for a constant c′ > 0, τ(k) ≤ e−c
′kη1 . Lemma 2.1 ([MPR11]). Let {Wj}j≥1 be a sequence of zero-mean real-valued random variables satisfying conditions (a), (c)-(i), and (d) of Assumption 2.1. Define κM (x) = (x ∧M) ∨ (−M), for some M , where (x ∧ y) = min(x, y), and (x ∨ y) = max(x, y), and set,
V := sup M≥1 sup i>0
( VAR(κM (Wi)) + 2 ∑ j>i |COV (κM (Wi),κM (Wj))| ) . (8)
Note that V is finite. Then, for any N ≥ 4, there exist positive constants C1, C2, C3, and C4 depending only on c, η1, η2 such that, for any t > 0, we have
P ( sup j≤N ∣∣∣∣∣ j∑
i=1
Wi ∣∣∣∣∣ ≥ t ) ≤ Ne− tη C1 + e− t2 C2NV + e − t2C3N exp ( tη(1−η) C4(log t) η ) . (9)
To deal with polynomially tailed interactions, i.e., under condition (c)-(ii) of Assumption 2.1, we prove a concentration inequality for the sum of exponentially β-mixing random variables with polynomially heavy-tails, which may be of independent interest. Lemma 2.2 (Concentration for heavy-tailed β-mixing sum). Let {Wj}j≥1 be a sequence of zeromean real valued random variables satisfying conditions (a) and (c)-(ii) of Assumption 2.1, for some η2 > 2. Then for any positive integer N , 0 ≤ d1 ≤ 1, and d2 ≥ 0, and for any t > 1, we have,
P ( sup j≤N ∣∣∣∣∣ j∑
i=1
Wi ∣∣∣∣∣ ≥ t ) ≤ 2 η2+3
(d2 log t) 1−η2 η1
N
t(1+d1(η2−1)) + 8
N
t(1+c′d2) + 2e−
t2−2d1 (d2 log t) 1/η1
9N ,
(10)
where c′ > 0 is a constant.
Note that we do not need condition (d) of Assumption 2.1 for Lemma 2.2. Since the tail probabilities decay polynomially and the mixing coefficients decay exponentially fast, the effect of heavy-tail dominates and hence, there is no trade-off between η1 and η2. Lemma 2.2 extends the results of [BMdlP20] (on iid heavy-tailed random variables) to the exponentially β-mixing setting. The results of [BMdlP20] show that even for the iid case, the tail probability of the sum decays polynomially with t. We show that similar polynomial tail bounds can be obtained (up to log factors) even in the dependent setting. Furthermore, when β(k) = 0, k > 0, the sequence is iid, in which case we recover the result of [BMdlP20]. We also remark that the above two results are crucial to derive our convergence rates for ERM. As a preview, when the DGP has only (2 + δ)-moments, δ > 0, Lemma 2.2 leads to risk bounds that hold with polynomial probability, whereas for sub-Weibull DGP (see Definition 4), Lemma 2.1 would lead to risk bounds that hold with exponential probability.

heading: 3 Main Results
text: In this section, we state our main results on the rates of convergence of ERM for both squared and convex loss functions. First, we consider the squared loss. Theorem 3.1 (Rates of ERM with squared loss). Consider the ERM procedure with the squared error loss. For τ0 < τ2QH(2τ)/8, setting µ = NrQH(2τ)c 1 η1 /4, for some constants c, c′ > 0, and
0 < r < 1, we have, for sufficiently large N , and some positive constants C̃1, C̃2, the following:
1. Under conditions (a), (b), (c)-(i), and (d) of Assumption 2.1, for 0 < ι < 1/4, ∥f̂ − f∗∥L2 := (∫ (f̂ − f∗)2dπ ) 1 2 ≤ max { N− 1 4+ι, ωµ ( F − F , τQF−F (2τ)
16
)} , (11)
with probability at least (for V as defined in (8))
1− C̃1NrQH(2τ)c 1 η1 exp(−N (1−r)η1)− C̃2Nexp ( −(N 12+2ιτ0)η ) . (12)
2. Under conditions (a), (b), and (c)-(ii) of Assumption 2.1, for 0 < ι < (1− 1/η2)/4, ∥f̂ − f∗∥L2 ≤ max { N − 14 ( 1− 1η2 ) +ι , ωµ ( F − F , τQF−F (2τ)
16
)} . (13)
with probability at least
1− C̃1NrQH(2τ)c 1 η1 exp(−N (1−r)η1)− C̃2τ − 2η21+η2 0 N − 4ιη21+η2 . (14)
The detailed expression of the probabilities are provided in the Appendix (Theorem B.1). Remark 1. To the best of our knowledge, the above result is the first result on understanding rates of convergence of ERM with squared error loss functions with unbounded noise (as well as the loss) for heavy-tailed dependent data. For a wide range of function classes F used in practice (see Section 4), the dominant term in the rate of convergence is N− 1 4+ι, for 0 < ι < 1/4. Furthermore, under the stronger condition (c)-(i) of Assumption 2.1, the risk bound holds with exponential probability, whereas under the weaker condition (c)-(ii), it holds only with polynomial probability.
We now show that when the small-ball condition in Assumption 2.1 is replaced with the stronger normequivalence assumption, also considered in [MZ20], one could obtain improved rates. Examples of random vectors that satisfy the norm-equivalence conditions include multivariate student t-distribution and sub-exponential random variables. We refer to Section 4 for illustrative examples. Assumption 3.1 (Lp − L2 norm-equivalence). Let F ⊂ Lq(π) be a class of functions for some q ≥ 3. The function class F −F = {f − h : f, h ∈ F} is Lp −L2 norm-equivalent for some p > 2, if there exists an M1 > 0 such that, ∥h∥Lp := ( ∫ |h|p dπ)1/p ≤M1∥h∥L2 , ∀h ∈ F − F . Corollary 3.1. For the ERM procedure with squared error loss, under Assumptions 2.1, with condition (b) replaced by Assumption 3.1 with p = 8, for some 0 < ι < 12 and r, µ and τ0 same as in Theorem 3.1, for sufficiently large N , we have
∥f̂ − f∗∥L2 ≤ max { N− 1 2+ι, ωµ(F − F , τQF−F (2τ)/16) } , (15)
with probability at least (for some constants C̃1 and C̃2)
1− C̃1NrQH(2τ)c 1 η1 exp(−N (1−r)η1)− C̃2Nexp ( −(N2ιτ0)η/M1 ) .
Remark 2. In the model-based nonparametric regression setting (as discussed in Related Works) with Xi allowed to be dependent on ξi for all i = 1, . . . , n, but (Xi, ξi) being iid, in [HW19, Proposition 3] authors show a lower bound of N− 1 2+ϵ for some ϵ > 0, for sufficiently heavy-tailed input Xi. The above result provides an upper bound of similar order, for a more general setting in comparison to [HW19]. We also remark that for Corollary 3.1, in the model-based setting, if we assume that ξi is independent of Xi, we have the same conclusion with just p = 4 instead of p = 8. Furthermore, note that in part 2 of Theorem 3.1, L8 norm does not exist for η2 ≤ 8. We have elaborated more on how this result compares Theorem 3.1 in [Men15] later in Section F.
We now present our results for the class of convex loss functions that are locally strongly-convex. Assumption 3.2 (Convex loss). The loss function ℓ : R → R+∪{0} is a convex loss function which is strongly convex in the neighborhood of 0, i.e., there exists a t2 > 0 such that for any x, y ∈ [−t2, t2], ℓ(y) ≥ ℓ(x) + ℓ(x)′(y − x) + µc(y − x)2/2 for some constant µc > 0. Theorem 3.2 (Rates of ERM with convex loss). Consider ERM with loss functions that satisfy Assumption 3.2. For τ0 < c2QF−F (2τ)ρ(0, t2)τ2, t2 = O((κ0 + 1/ √ QH(2τ))∥ξ∥L2), setting µ = Nη1/(1+η1), for some constants c, c′ > 0, we have, for any N ≥ 4, the following:
1. Under conditions (a), (b), (c)-(i), and (d) of Assumption 2.1, for 0 < ι < 14 , ∥f̂ − f∗∥L2 ≤ max { N− 1 4+ι, 2ωQ(F − F , N, ζ1, ζ2) } , (16)
with probability at least (for V is defined in (8) and some positive c9, c10, C̃3)
1− c9QH(2τ)1− 1 η1Nη1/(1+η1)e−c10QH(2τ)
1+ 1 η1 N η1 1+η1 − C̃3Nexp ( −(N 12+2ιτ0)η/C1 ) .
2. Under conditions (a), (b), and (c)-(ii) of Assumption 2.1, for 0 < ι < (1− 1/η2)/4, ∥f̂ − f∗∥L2 ≤ max { N− (1−1/η2) 4 +ι, 2ωQ(F − F , N, ζ1, ζ2) } , (17)
with probability at least (for constants c9, c10, C̃4 > 0)
1− c9QH(2τ)1− 1 η1Nη1/(1+η1)e−c10QH(2τ) 1+ 1 η1 Nη1/(1+η1) − C̃4τ − 2η21+η2 0 N − 4ιη21+η2 . (18)
To the best of our knowledge, the above result is the first result on understanding rates of convergence of ERM with convex loss functions with unbounded noise (as well as the loss) for heavy-tailed dependent data. The above result highlights the advantage of using a robust loss function, e.g. Huber loss, over a quadratic loss function. For example, if (f − f∗)(X) has a sub-Weibull tail and the noise ξ has polynomial tail, one can still obtain risk bounds with exponential probability. This is because in this case ℓ′(ξ)(f − f∗)(X) can still be sub-Weibull for a suitable chosen ℓ′(ξ). Such a situation arises, for example, when there are outliers even if the data is light-tailed. With a squared error loss, one won’t be able to obtain a risk bound with exponential probability in this scenario. We will illustrate this in Section 4.3 through Huber loss, a popular choice of robust loss function in robust statistics. Similar to the quadratic case, we also have an improved result when the small-ball condition is replaced with the norm-equivalence condition. Due to space constraints, we state and prove it in the Appendix B.2.

heading: 4 Illustrative Examples
text: We illustrate the results of Section 3 with three examples, based on sub-Weibull random variables and Pareto random variables, that are canonical models of heavy-tailed data in the literature.

heading: 4.1 Example 1:β-mixing Sub-Weibull DGP with Squared Error Loss
text: Here, we consider sub-Weibull random variables to model the heavy-tailed behavior in the DGP. Definition 4 (Sub-Weibull random vectors). A real-valued random variable X is said to be subWeibull with parameter η > 0, if there are constants K1,K2 > 0, such that we have P (|X| > t) ≤ 2exp (−(t/K1)η) , or equivalently ∥X∥p = E [|X|p]1/p ≤ K2p 1 η . Based on this, a random vector X ∈ Rd is said to be marginally sub-Weibull with parameter η > 0 if each coordinate of X is sub-Weibull with η. We use X ∼ SW(η) to represent this fact.
The above family of distributions define a rich class of random variables, allowing for heavier tails than sub-Gaussian tails (η = 2) or sub-exponential tails (η = 1). Let {δi}i∈Z+ be an iid sequence of d-dimensional random vectors with independent coordinates with δi ∼ SW(ηδ). Assume that the dependent input vectors are generated according to the model
Xi = AXi−1 + δi, (19)
where A ∈ Rd×d with spectral radius less than 1. For simplicity, let A = σ20Id where σ20 < 1, and {Yi}i∈Z+ be a univariate response sequence given by Yi = θ∗⊤Xi + ξi, where θ∗ ∈ Bd1 (R) belongs to the ℓ1-norm ball in Rd with the radius R, and {ξi}ni=1 is an i.i.d sequence independent of Xi ∀i, and ξi ∼ SW(ηξ) for 0 < ηξ < 1 has independent coordinates. To proceed with learning framework, we consider ERM with squared loss and the function class F := FR = { ⟨θ, ·⟩ : θ ∈ Bd1 (R) } . We denote the difference function class FR −FR by HR. We verify conditions (a), (b), (c), and (d) of Assumption 2.1, and Assumption 3.1, provide the rate of ERM in this setting, and compare our result to that of [WLT20] in a similar setting in Appendix C.

heading: 4.2 Example 2: β-mixing Pareto DGP with Squared Error Loss
text: Let X̃t,i denote the i-th coordinate of the vector X̃t ∈ Rd. We consider the process given in [Pil91]: For i = 1, 2, · · · , d, η3 > 2 + 2ι, where ι > 0 is a small number, and t = 0, 1, · · · , define
X̃t,i =
{ 2 1 η3 X̃t−1,i with probability 1/2
min ( 2 1 η3 X̃t−1,i, δt,i ) with probability 1/2
(20)
where {δt,i}i=1,2,··· ,d,t=1,2,··· is a sequence of iid Pareto random variables with the distribution L+(δ; η3, di) = η3(diδ) η3−1/ (1 + (diδ) η3)
2 for δ > 0, η3 > 2 + 2ι, and we write X ∼ L+(η, σ) to denote that X is a Pareto random variable with parameters η and σ. The survival function of δ ∼ L+(η3, d) is given by, P(δ > t) = (1 + (dt)η3)−1, for t > 0. Let X̃1,t and X̃2,t be two independent trails of the process in (20). Let {Ut}t=0,1,··· be a sequence of iid U [0, 1] random variables. Now consider the process Xt,i = X̃1,t,i1(Ut ≤ 1/2)− X̃2,t,i1(Ut > 1/2). The marginal distribution of Xt,i is given by symmetric Pareto distribution, i.e.,
L(x; η3, di) = η3(di|δ|)η3−1
2 (1 + (di|δ|)η3)2 −∞ < δ <∞, η3 > 2 + 2ι. (21)
Now, let {Yi}i∈Z+ is given by Yi = θ∗⊤Xi + υi, where θ∗ ∈ Bd1 (R), as before. Let {υi}ni=1 be an iid sequence of 0 mean random variables independent of Xi for all i, with heavy tails such that for all t > 0, P (|υi| ≥ t) ≤ 1/(1 + tη4), η4 > 2 + 2ι. Like in Example 4.1, we consider ERM with squared loss and the function class F := FR = { ⟨θ, ·⟩ : θ ∈ Bd1 (R) } , where Bd1 (R) denotes the d-dimensional ℓ1-ball with radius R. We denote the difference function class FR −FR by HR. We show that conditions (a), (b), and (3)-(ii) of Assumption 2.1 hold here in the following section.

heading: 4.2.1 Verification of Assumption 2.1 for Example 4.2
text: [Pil91] shows that the AR(1) process given by (20) is strictly stationary if X̃0,i ∼ L+(η3, di), and X̃0,i are independent of each other for i = 1, 2, · · · , d. The stationary distribution is given by L+(η3, di). Let π(Xt1 , Xt2 , · · · , Xtn) be the joint distribution of Xt1 , Xt2 , · · · , Xtn for a set of time points t1, t2, · · · , tn. Now for any positive integer k,
π({Xti+k}ni=1) = π({X̃1,ti+k1(Uti+k ≤ 1/2)− X̃2,ti+k1(Uti+k > 1/2)}ni=1) =π({X̃1,ti1(Uti ≤ 1/2)− X̃2,ti1(Uti > 1/2)}ni=1) = π({Xti}ni=1). (22)
The second equality above follows from the fact X̃1,t and X̃2,t are strictly stationary process and {Ut} is iid. So Xt is a strictly stationary process with marginal distribution of Xt,i given by symmetric Pareto distribution
L(x; η3, di) = η3(di|δ|)η3−1
2 (1 + (di|δ|)η3)2 −∞ < δ <∞, η3 > 2 + 2ι. (23)
Without loss of generality we will assume that d1 ≤ d2 ≤ · · · ≤ dd, ∑d
i=1 1/di = K0, and d1 ≥ C ′6 for some constants K0, C ′6 > 0. It is shown in Lemma 1 of [Ris08], that the AR(1) process in (20) is ϕ-mixing with ϕ(k) = k log 2/(2k − 1), k = 1, 2, · · · . where ϕ-mixing coefficients are defined as in [Bra05]. We also have β(k) ≤ ϕ(k) ≤ e−k/3. Since Xt depends only on X1,t, X2,t, and Ut, X1,t, and X2,t are independent and exponentially β-mixing, and Ut is iid, Xt is also exponentially β-mixing.
• Since Yi depends only on Xi, and υi are iid, {(Xi, Yi)} is a strictly stationary β-mixing sequence with β(k) ≤ e−k/3. So condition (a) of Assumption 2.1 is true here.
• Now we will verify condition (b) of Assumption 2.1. Let σX,p denote ∥X∥p for p > 0. Then
E [( θ⊤X )2] = d∑ i=1 θ2i σ 2 Xi,2 d2i ≥ σ20 d∑ i=1 θ2i d2i , (24)
where σ0 = mini σXi,2, i = 1, 2, · · · , d. Since Xi ∼ L(η3), for any θi ∈ R we have ∥θiXi∥η3−0.5ι ≤ K1|θi|σX,η3−0.5ι/di. So,
∥θ⊤X∥η3−0.5ι ≤ K1σX,η3−0.5ι d∑
i=1
|θi| di ≤ K1σX,η3−0.5ι
∑d i=1 |θi| di
σ0 √∑d i=1 θ2i d2i
∥θ⊤X∥2 (25)
≤K1σX,η3−0.5ι √ d
σ0 ∥θ⊤X∥2. (26)
Then from Lemma 4.1 of [Men15] we have that the condition (b) of Assumption 2.1 is true here.
• Since υi and Xi are independent, and E [υi] = 0, we have E [ υit ⊤Xi ] = 0. Let η2 =
min(η3, η4). Then using Markov’s inequality, for any t ∈ HR and ∀i, we have, P (∣∣υit⊤Xi − E [υit⊤Xi]∣∣ ≥ τ) ≤ ∥υi∥η2η2 (∑d j=1 ∥tjXi,j∥η2 )η2 τη2 ≤ (RσυσX,η2dH) η2 τη2 ,
where dH = ∑d i=1 d −1 j , and for all i, j, ∥υi∥η2 = συ, and ∥Xi,j∥η2 = σX,η2/dj . This
implies that condition (c)-(ii) of Assumption 2.1 is true in this setting. Proposition 4.1. Consider the learning problem described in Section 4.2. Then with probability at least (14), we have,
∥f̂ − f∗∥L2 ≤ max { N − 14 ( 1− 1η2 ) +ι ,
C9R
τQHR(2τ) 3/2
d1/(η2−0.5ι)+ι/8N−1/2+ι } . (27)
The proof of Proposition 4.1 could be found in the Appendix D. Remark 3. Note that the heaviness of the tail dominates the exponential β-mixing rate η1 in determining the rates of convergence. Observe that as η2 → ∞ all the moments exist, and N−(1−1/η2)/4+ι → N−1/4+ι. This rate is the same as the one that we obtain under condition (c)-(i), although with weaker polynomial probability. Furthermore, the dimension dependency is polynomial here. On a related note, [ZZ18] analyzes ℓ1-regression with a truncated loss in the iid setting with the assumption η2 > 2 and gets √ d/N rate with exponential probability. We point out that in the iid setting Medians-of-mean method achieves the optimal rate with exponential probability under the stronger assumption that log d moments exist [LM19].

heading: 4.3 Example 3: β-mixing sub-Gaussian data and Pareto noise with Huber Loss
text: In this example, we consider Huber loss which satisfies Assumption 3.2:
ℓTh(t) = { t2/2 if |t| ≤ Th Th |t| − T 2h/2 if |t| ≥ Th .
We now formally establish the benefits of using Huber loss when the noise ξ has a polynomial tail but (f − f∗)(X) has a sub-Weibull tail. Let {δi}i∈Z+ be an iid sequence of d-dimensional standard
Gaussian random vectors δi ∼ N(0, Id). To compare the performance of ERM under Huber loss to that of squared loss, for simplicity, we allow X to be Gaussian; but a similar result will hold for sub-Weibull δi. Assume that the input vectors are generated according to (19) where A is a d× d matrix with spectral radius less than 1. For simplicity, let A = σ20Id where σ 2 0 < 1. [WZLL20] showed that this time series is stable, strict sense stationary, with Xi ∼ N(0, 1/(1 − σ40)Id). Let {Yi}i∈Z+ , be the response sequence given by Yi = θ∗⊤Xi + ξi, where θ∗ ∈ Bd1 (R), ℓ1-norm ball in Rd. Let {ξi}ni=1 be independent of Xi for all i, and be an iid sequence of 0 mean random variables with heavy tails such that for all t > 0, P (|ξi| ≥ t) ≤ 1/(1 + tη4), η4 > 2 + 2ι and VAR(ξ) ≤ σξ. Set Th = 3σξ. We verify the required assumptions for this example in the following section.

heading: 4.3.1 Verification of Assumptions for Example 4.3
text: Condition (a) of Assumption 2.1 and condition (b) of Assumption 2.2 are true here by the same argument as in Example 4.1. This also implies that Assumption 3.1 is true in this case for p = 8. Since we consider Huber loss, a lipschitz continuous loss, we have l′(ξ) ≤ min(|ξ/2|, Th). Set Th = c11(κ0 + 1/ √ ϵ)(σξ + 2R). As an immediate consequence of [VGNA20, Proposition 2.3], we have that condition (c)-(i) in Assumption 2.2 is valid here with η2 = 2. Note that if a sequence is exponentially β-mixing, i.e., satisfies condition (a) of Assumption 2.1 with coefficient η1 > 0, then the same condition is true for all η < η1. So choosing η′1 < η2/(η2 − 1) we get Condition (d) of Assumption 2.1 is true. Proposition 4.2. Consider the learning problem described in Section 4.3. Then, for some constant c12 > 0, with probability at least
1− c9ϵ1− 1 η1Nη1/(1+η1)e−c10ϵ 1+ 1 η1 Nη1/(1+η1) − C̃2Nexp ( −(N2ιτ0)η/M1 ) ,
we have ∥f̂ − f∗∥L2 ≤ max { N− 1 2+ι, c12R √ log(ed/N)N− 1 2 } .
The proof of Proposition 4.2 could be found in the Appendix E. Remark 4. In this setting, using squared loss would mean that condition (c)-(ii) of Assumption 2.1 is true. So, by part 2 of Theorem 3.1, the obtained rate would be of order N−1/8 with polynomial probability given by (14), which is significantly worse than that of Proposition 4.2 with Huber loss.

heading: 5 Proof Sketch of Theorem 3.1 and 3.2
text: Recall the decomposition (1). The first and the last terms in the RHS of (1) are handled respectively by condition (b) and (c) of Assumption 2.1. The idea is to show that if for some f ∈ F , ∥f − f∗∥L2 is large, then with high probability T1 := N−1 ∑N i=1(f − f∗)2(Xi) ≥ B (Lemma B.1) and
T2 := 2N −1∑N
i=1 ξi(f − f∗)(Xi) ≥ B̄ (see (67)) where B + B̄ > 0. But since f̂ minimizes PNLf and f∗ ∈ F , PNLf ≤ 0. So with high probability ∥f̂ − f∗∥L2 is small. In contrast to [Men15], we face two major challenges: 1- For the lower bound on T1, the symmetrization argument used in the iid case (e.g. [Men15, Men18]) is not applicable under our dependency structure. 2- For the term T2, [Men15, Men18] use the complexity measure α∗N (γ, δ) (see (108)) in the iid case to control the noise-input interactions, which is not possible to do in our setting; our analysis to control T2 is different, through which we can show how different tail conditions on the data and noise affect the high-probability statements on the learning rate. The proof of Theorem 3.2 for the locally strongly-convex loss functions, follows a similar strategy although the technical details are more involved.

heading: 6 Conclusion
text: In this work, we analyzed the performance of ERM with squared and convex loss functions, when the DGP is both dependent (specifically, exponentially β-mixing) and heavy-tailed. We derived explicit rates using a combination of small-ball method and concentration inequalities. We demonstrated the applicability of our results on a high-dimensional linear regression problem, and showed that our assumptions are easily verified for a certain classes of sub-Weibull and Pareto DGP. Our results clearly show the benefits of using Huber loss over the squared error loss for ERM with heavy-tailed data in our setting. For future work, we plan to study median-of-means based techniques and examine establishing similar rates of convergence for dependent heavy-tailed data.

Review summaries:
This submission is a technical paper, which is concerned with obtaining bounds for the Empirical Risk Minimization (ERM) procedure for regression, when the data is both heavy-tailed and dependent (non-iid). The performance of ERM in the heavy-tailed context is by now fairly well-understood, in part thanks to the observation that two-sided concentration for the empirical risk is not needed to establish excess risk bounds for ERM, and that one-sided guarantee suffice (functions with large population risk have large empirical risk). The authors' contribution is to extend previous results in the iid setting to the dependent (specifically, mixing or exponentially-mixing) setting. While both the heavy-tailed and dependent settings have been considered separately, this works combines these two cases.  As in the iid case, the analysis involves both obtaining lower isomorphic (small ball) estimates for the function class, and controlling the interaction between the noise and the functions in the class (the so-called multiplier inequality).  The main applications of the generic results obtained for square and convex loss (Section 3) are described in Section 4. For regression over L^1 balls, the authors obtain bounds under square loss for both exponential and polynomial tails of the noise, as well as bounds for Huber loss under polynomially-tailed noise (but sub-Gaussian data).

The paper considers the performance of ERM on heavy-tailed and dependent data. It proves general convergence rates with a few illustrative examples.  Update  ---- I have read the author(s)' response to my review and maintain my rating. 

The paper analyses empirical risk minimization over a locally strongly convex function class in the context where the data generating process returned non-IID data, and can be potentially heavy-tailed. The authors give specific realizations of their analysis to sparse linear function class.  This work aims to bridge the gap between analyses that handle IID data and heavy-tailed data separately.

The authors derive risk bounds for empirical risk minimization under heavy-tailed $\beta$-mixing observations. They provide general bounds for convex, locally strongly convex, loss functions and tighter bounds for the squared loss. They give three illustrative models for which they provide explicit bounds, the most important one being the Huber loss for which they find a justification of its usage as a robust lost in presence of sub-Weibull noise distributions.