Statistical Inference for Stochastic Gradient Descent
arXiv:2605.26000v1 [stat.ML] 25 May 2026
Beyond Finite Variance Jose Blanchet, Peter Glynn, and Wenhao Yang Management Science and Engineering, Stanford University
May 26, 2026
Abstract Stochastic gradient descent (SGD) is a foundational algorithm for large-scale statistical learning and stochastic optimization. However, statistical inference based on SGD iterates remains challenging when stochastic gradients have infinite variance, as the relevant limiting distributions depend on unknown nuisance parameters. In this paper, we develop an efficient, model-agnostic methodology for constructing confidence regions from SGD trajectories that applies in both finite- and infinitevariance regimes. The procedure is based on a joint weak convergence result for the Polyak–Ruppert averaged estimator and an empirical second-moment normalizer constructed from stochastic gradients along the SGD trajectory. This joint limit yields a self-normalized statistic in which the leading tail-dependent scaling terms cancel. We then use a subsampling calibration scheme to estimate the relevant critical values, avoiding explicit estimation of tail indices, slowly varying functions, or stable-law parameters. The resulting confidence regions are straightforward to implement and are asymptotically valid under both the finite- and infinite-secondmoment regimes. Simulation studies show reliable coverage in various settings, supporting the proposed method as a practical tool for uncertainty quantification in stochastic optimization.
1
1
Introduction
Many parameter estimation problems in statistics can be formulated as stochastic optimization problems. Given a convex population loss function ℓ(·), we consider θ∗ ∈ argmin ℓ(θ). θ∈Rd
A classical example is M-estimation, where the loss is an empirical average of a sample-wise P criterion, ℓ(θ) := N1 N i=1 f (θ, Xi ). Building on this optimization viewpoint, Robbins and Monro [29] introduced stochastic approximation (SA), in which the parameter is updated iteratively according to
θn+1 = θn − ηn g(θn , ξn+1 ), where E[g(θ, ξ)] = ∇ℓ(θ), ηn = c · n−ρ with ρ ∈ (1/2, 1], and ξn denotes the randomness in the stochastic gradient query. This recursion is commonly referred to as stochastic gradient descent (SGD), and it is now widely used in large-scale statistical learning and stochastic optimization. A large body of work has studied its convergence properties in various settings [4, 11, 22, 23, 25, 26, 30]. However, convergence guarantees do not quantify the statistical uncertainty of the estimator produced by the SGD trajectory. This motivates the problem of constructing asymptotically valid confidence regions for θ∗ from the iterates generated by SGD. Formally, we aim to construct a confidence region CIn for the true parameter θ∗ with nominal level 1 − δ such that lim P (θ∗ ∈ CIn ) = 1 − δ,
n→+∞
where n is the number of SGD iterations. A standard route to asymptotically valid confidence regions is to first establish a central limit theorem and then estimate the relevant quantiles of the limiting distribution. In early work, Sacks [33] established 2
asymptotic normality for the final SGD iterate under suitable conditions. However, this limit theorem is difficult to use directly for inference: the asymptotic variance is complicated to estimate, and the constant c in the learning rate is hard to tune. Subsequently, Polyak and Juditsky [27] established asymptotic normality for Polyak–Ruppert averaged SGD. A key advantage of Polyak–Ruppert averaging is that its asymptotic distribution does not depend on the choice of the learning-rate constant. Building on this result, recent work has developed several inference methodologies for SGD, including plug-in estimators [6], batch means [6, 42], random scaling [16], replication [41], and bootstrap methods [40]. These methods, however, typically rely on the assumption that the stochastic gradients have finite second moments. While analytically convenient, this assumption can be inconsistent with empirical findings in deep learning, where stochastic gradient noise often exhibits heavy-tailed1 behavior [35, 36]. This gap between theory and practice raises important questions about the reliability of classical inference procedures under more realistic noise regimes. Heavy-tailed stochastic gradients may arise from several sources, including the structure of neural networks, the choice of hyperparameters such as learning rate and batch size, and rare but influential data points [35, 36]. When the classical finite-second-moment assumption is violated, standard uncertainty quantification procedures can fail in several ways: • The classical central limit theorem [27] may fail to hold. • Plug-in variance estimators [6] become inconsistent. • Bootstrap methods [1] exhibit erratic behavior and fail to converge. These failures make existing approaches for constructing confidence regions unreliable in the presence of heavy-tailed noise. Recently, Blanchet et al. [5] showed that, under suitable 1
In this manuscript, a distribution is referred to as heavy-tailed if it does not have a finite second moment, and light-tailed if it has a finite second moment.
3
conditions, the final iterate of SGD converges weakly to a non-Gaussian distribution: ∗
d
Z +∞
h(ηn ) (θn − θ ) → 0
where h(·) is a
1(ρ = 1) exp − H − t dLt , 2c
1 − 1 -regularly varying function and Lt is a Lévy process with index α. α
This limit theorem provides a theoretical foundation for understanding infinite-variance SGD, but it does not immediately yield implementable confidence regions because of several difficulties: • The choice of c in the learning rate ηn = c · n−ρ depends on the Hessian H. • regularly varying function h(·) and tail index α depend on the detailed tail behavior of the stochastic gradient g(θ, ξ). • The quantiles of the limiting distribution depend on several nuisance parameters. In certain settings where heavy-tailed noise is artificially injected into the iterates, the nuisance parameters and limiting distribution become explicitly known, making statistical inference more direct. We discuss this approach in Appendix B. However, such approaches require careful calibration of the injected perturbations and can obscure the intrinsic stochastic behavior of the learning process. It is therefore desirable to develop inference procedures that do not rely on artificially injected randomness. In this work, we address these challenges by developing a model-agnostic methodology for constructing confidence regions from SGD trajectories. The proposed procedure does not require artificial perturbations, tail-index estimation, or a parametric specification of the stochastic gradient noise. Our main contributions are as follows: • We establish the joint weak convergence of the Polyak–Ruppert averaged SGD estimator and an empirical second-moment normalizer computed from the same SGD trajectory. This result enables the construction of a self-normalized statistic in which the unknown scaling function h(·), the tail index α, and other leading nuisance quantities cancel out. As a result, the statistic can be computed from 4
the observed SGD trajectory without estimating the detailed tail behavior of the stochastic gradients. • To avoid direct evaluation of the non-Gaussian limiting distribution, we introduce a subsampling calibration scheme based on short auxiliary SGD trajectories. This calibration does not require repeated independent runs of the full n-step SGD procedure and avoids explicit estimation of the tail index, slowly varying function, or stable-law parameters. And this procedure applies without knowing whether the stochastic gradients have finite or infinite variance. We support the theory with experimental studies that examine coverage and interval length across a range of tail behaviors. Together, these contributions bridge the gap between heavy-tailed limit theorems and statistical inference for SGD, yielding a practical tool for uncertainty quantification in stochastic optimization. Related Literature For weak convergence of heavy-tailed SGD, Krasulina [14] first established a limit theorem in the one-dimensional case with Pareto-type noise. More recently, Blanchet et al. [5] extended this result to high dimensions and relaxed the assumptions on the stochastic gradients to cover a broader class of heavy-tailed distributions. In related work, Goodsell and Hanson [9], Li [17] investigated conditions for almost sure convergence, while Wang et al. [38] analyzed the weak convergence of Polyak–Ruppert averaged SGD under Pareto-type noise, without giving a detailed characterization of the limiting distribution or developing inference methodology. Beyond limit theorems for SGD, heavy-tailed phenomena in stochastic optimization have also been studied extensively. Both theoretical analyses and empirical evidence show that hyperparameters such as learning rate, batch size, and loss structure can induce heavy-tailed behavior [8, 10, 12, 13, 34]. On the algorithmic side, Cutkosky and Mehta [7], Liu and Zhou [18], Wang et al. [38] investigated non-asymptotic convergence rates under various noise assumptions. Moreover, several works demonstrate that heavy-tailed 5
behavior is positively correlated with generalization properties across different learning settings [3, 20, 37]. In particular, for SGD, Simsekli et al. [35, 36] showed that heavy-tailed dynamics favor wider local minima over sharper ones, and Wang et al. [39] proposed truncated heavy-tailed SGD to mitigate convergence to sharp minima. In the broader literature on heavy-tailed statistics, Logan et al. [19] introduced selfnormalization techniques that cancel dominant nuisance parameters in heavy-tailed mean estimation problems. Building on this idea, Romano and Wolf [31] proposed a subsampling approach to estimate quantiles of the limiting distribution of self-normalized statistics in i.i.d. heavy-tailed mean estimation. More recently, Bai et al. [2] extended this framework to accommodate long-range dependent time series. d
Notation. Let R̄ := R ∪ {±∞} and R̄d := (R̄)d . We write → for convergence in v
distribution and → for vague convergence. The (d − 1)-dimensional unit sphere is denoted by Sd−1 := {x ∈ Rd : ∥x∥ = 1}. We denote by DJ1 ([0, 1], Rd ) the Skorokhod space endowed with the J1 topology. Sd++ denotes the set of all d × d symmetric positive definite matrices. We write a ≲ b if there exists C > 0 such that a ≤ Cb.
2
SGD Setup and Limit Theory
In this section, we introduce the stochastic gradient descent setup and the assumptions used throughout the paper. We then summarize the relevant limit theory for heavy-tailed SGD and explain why these results alone do not immediately yield a practical inference procedure. Foundational results from heavy-tailed statistics are deferred to Appendix A.
2.1
Assumptions
Given an optimization problem
min ℓ(θ),
θ∈Rd
6
where ℓ(θ) is a strongly convex function, stochastic gradient descent iteratively approximates the minimizer θ∗ ∈ argminθ ℓ(θ) via
θn+1 = θn − ηn g(θn , ξn+1 ), where Eξ∼P [g(θ, ξ)] = ∇ℓ(θ) and {ξi }i≥1 are i.i.d. random variables. Motivated by empirical observations that stochastic gradients may exhibit infinite variance [36], we assume that the stochastic gradients satisfy a multivariate regular variation condition. Assumption 2.1 (Blanchet et al. [5]). For any θ ∈ Rd , g(θ, ξ) satisfies: P (∥g(θ, ξ)∥ > r) = r−α b0 (r, θ), g(θ, ξ) v ∈ · ∥g(θ, ξ)∥ > r → µ(θ, ·), when r → +∞. P ∥g(θ, ξ)∥ We also impose regularity assumptions on the loss ℓ(θ) and the stochastic gradient g(θ, ξ). Assumption 2.2 (Wang et al. [38]). The set {∇2 ℓ(θ) : θ ∈ Rd } is bounded and uniformly p-positive definite for any p ∈ (1, α). Equivalently, for any θ ∈ Rd and ∥u∥p = 1, ⊤
∇2 ℓ(θ) is bounded and satisfies u⊤ ∇2 ℓ(θ) (sgn(u1 )|u1 |p−1 , · · · , sgn(ud )|ud |p−1 )
> 0 for
any p ∈ (1, α). Assumption 2.3 (Blanchet et al. [5]). There exists CLip > 0 such that for any θ1 , θ2 :
∥g(θ1 , ξ) − g(θ2 , ξ)∥ ≤ CLip ∥θ1 − θ2 ∥, a.s. Assumption 2.4 (Blanchet et al. [5]). There exists q ∈ (1, α), K > 0, such that for any θ,
∥∇ℓ(θ) − ∇2 ℓ(θ∗ )(θ − θ∗ )∥ ≤ K∥θ − θ∗ ∥q .
We denote H := ∇2 ℓ(θ∗ ). 7
Assumption 2.2 strengthens the standard strong convexity condition and is commonly used in the analysis of heavy-tailed stochastic optimization [38]. Assumption 2.3 imposes a smoothness requirement on the stochastic gradients, while Assumption 2.4 bounds the approximation error of the nonlinear loss function ℓ(·). Similar assumptions are also adopted in Polyak and Juditsky [27], Wang et al. [38].
2.2
Limit Theory for SGD Iterates
In prior work [5], a central limit theorem for the final iterate of SGD has already been established, as summarized below. Theorem 2.1 (Blanchet et al. [5]). Under Assumptions 2.1, 2.2, 2.3, and 2.4, if the learning rate satisfies ηn = c · n−ρ (ρ ∈ (α−1 , 1]), then we have:
n
1 ρ(1− α )
ρ
∗
d
d
1 1− α
Z ∞
b1 (n ) (θn − θ ) → Zfinal,ρ = c
exp (−Hρ t) dLt , 0
where Lt is a Lévy process with characteristics 0, ν(θ , ·), − ∥x∥>1 xν(θ , dx) , and Hρ = ∗
R
∗
H − 1−αc 1(ρ = 1). −1
In practice, one naturally prefers algorithms that achieve the fastest possible convergence rate. Theorem 2.1 shows that, for SGD, the optimal rate is attained when the step size is chosen as ηn = c · n−1 . However, this choice requires the constant c to satisfy
c>
1 − α−1 , σmin (H)
which is difficult to verify or control in practice. To address this issue, we draw on the averaging technique of Polyak and Juditsky [27], which removes such parameter dependence and yields convergence guarantees without the delicate tuning of c. Theorem 2.2 (Polyak-Ruppert Averaging). Under Assumptions 2.1, 2.2, 2.3, and 2.4, if the learning rate satisfies ηn = c · n−ρ (ρ ∈ (α−1 , 1)) and denote θ̄n =
8
Pn
k=1 θk
n
, then we
have: d 1 n1− α b1 (n) θ̄n − θ∗ → ZPolyak , d
where ZPolyak = H −1 L1 and Lt is defined in Theorem 2.1. In the finite-variance setting, the Polyak–Ruppert averaged estimator for SGD is known to achieve the Cramér–Rao lower bound under mild regularity conditions [27], implying that it is the most efficient estimator attainable. In contrast, in the infinite-variance regime, such minimax efficiency criteria are no longer applicable, since the second moment is unbounded. Instead, what can be established is that the Polyak–Ruppert averaged estimator exhibits strictly better asymptotic performance than the final-iterate estimator within certain subclasses of heavy-tailed distributions. Theorem 2.3. If Lt is elliptically contoured with scale matrix Σ, i.e. log EX∼Lt [eiuX ] = α
−σα (L1 ; u)t with σα (L1 ; u) := (u⊤ Σu) 2 , and ΣH = HΣ, then the scales of the limit distributions Zfinal,1 in Theorem 2.1 and ZPolyak Theorem 2.2 satisfy, for any ∥u∥ = 1: (1)
σα (Zfinal,1 ; u) ≥ σα (ZPolyak ; u) .
We defer the detail definitions of elliptically contoured distribution and scale of a distribution to Appendix A. In Theorem 2.3, this comparison does not generally extend beyond the elliptically contoured stable subclass. If the stochastic gradient noise is completely skewed, the comparison in Eqn (1) may no longer hold.
2.3
The Inference Challenge
The limit theorems above identify the asymptotic behavior of SGD iterates in the heavytailed regime, but they do not provide an implementable confidence region. To turn 1
Theorem 2.2 into inference, one would need to know the normalization rate n1− α b1 (n) and the relevant quantiles of the stable limiting distribution. Both quantities depend on the 9
tail behavior of the stochastic gradients and are generally unavailable from the observed SGD trajectory. One situation in which these quantities may become accessible is when the heavytailed noise is artificially injected. To distinguish this case from the usual data-driven randomness, we decompose the stochastic gradient into two sources: intrinsic randomness and artificially injected randomness. For simplification, we represent the stochastic gradient as
g(θ, ξ) := g(θ, ξint ) + ξart ,
(2)
where ξint denotes the intrinsic randomness and ξart denotes the artificially injected randomness. For example, in Perturbed Gradient Descent (PGD), one may inject heaviertailed random noise ξart than g(θ, ξint ) into each gradient descent iterate. In this case, the parameters α, b1 (·), and the distribution of L1 in Theorem 2.2 are all determined by the injected noise and are therefore accessible. Consequently, Theorem 2.2 can be used directly to establish confidence regions for θ∗ , provided that a consistent estimator of H is available. We discuss this approach in Appendix B. However, artificially injecting noise is not always desirable. It requires prior knowledge of the intrinsic stochastic gradient noise in order to ensure that the injected component dominates the tail behavior, and it can obscure the intrinsic dynamics of the original SGD trajectory. We therefore seek an inference procedure that does not rely on artificial perturbations. In this setting, the nuisance parameters in Theorem 2.2 are determined by the intrinsic stochastic gradients, leading to two major obstacles: • The tail index α and the slowly varying function b1 (·) are unknown. • The quantiles of ZPolyak are unknown. In practice, it is often infeasible to verify whether the underlying data-generating mechanism follows an α-regularly varying distribution, and accurate estimation of b1 (·) is especially difficult. These obstacles motivate the self-normalized inference procedure 10
developed in Section 3.
3
Self-normalized Subsampling Inference
In this section, we develop a self-normalized inference procedure that addresses the difficulties identified in Section 2. The key idea is to normalize the Polyak–Ruppert averaged estimator by an empirical second-moment estimator computed from the same SGD trajectory. This normalization cancels the leading tail-dependent scaling factors, while subsampling provides a data-driven way to estimate the remaining critical values.
3.1
Self-normalized Statistic
Denote the empirical second-moment normalizer as n
1X Σn = g(θk , ξk+1 )g(θk , ξk+1 )⊤ . n k=1 We define the self-normalized statistic √ Tn⋆ (θ; h) :=
n · φ θ̄n − θ p , Tr(Σn )
(3)
where φ : Rd → R is a continuous homogeneous function with degree 1. The role of Σn is to provide a data-dependent normalization that cancels the unknown tail scaling in the heavy-tailed regime. The following theorem gives the formal justification. Theorem 3.1. Under Assumptions 2.1, 2.2, 2.3, and 2.4, the joint weak convergence holds:
n
1 1− α
b1 (n) θ̄n − θ
∗
2 1− α
,n
d b1 (n) Σn → H −1 L1 , W , 2
where W is infinitely divisible. For any continuous homogeneous function φ : Rd → R
11
with degree 1, by continuous mapping theorem, −1 L1 )| d |φ (H := Z ⋆ (φ). |Tn⋆ (θ∗ ; φ)| → p Tr(W )
In Theorem 3.1, the distribution of W need not be full rank. For instance, W fails to be full rank when some coordinates of g(θ∗ , ξ) have tail indices strictly larger than α. The benefit of Tn⋆ (θ∗ ; φ) is that it cancels out the dependence on the tail index α and the slowly varying function b1 (·). Consequently, if the (1 − δ) quantile value of Z ⋆ (φ) was known, a valid confidence region would be Rn (δ; φ) = θ ∈ Rd |Tn⋆ (θ; φ)| ≤ q(δ; φ) , where q(δ; φ) = inf {x |P (Z ⋆ (φ) > x) ≤ δ }. The remaining task is to estimate q(δ; φ) without evaluating the limiting distribution. Remark 3.1. In Eqn (3), there are many choices for the homogeneous function φ(·). Typical choices of φ(·) include linear projections and norms such as ∥ · ∥∞ and ∥ · ∥p . Similarly, the trace normalizer may also be replaced by other matrix functionals, such as the Frobenius norm or the largest eigenvalue, provided the corresponding limiting statistic is well defined.
3.2
Subsampling Calibration
To estimate q(δ; φ), we construct short auxiliary SGD trajectories to achieve it, which is motivated by the subsampling idea of [31]. Let tn = ⌊nr ⌋ for some r ∈ (0, 1) and (b)
(b)
let Bn = ⌊n/tn ⌋. For each sub-procedure b ∈ {1, . . . , Bn }, let θ̄tn and Σtn denote the averaged estimator and empirical second-moment estimator computed from a length-tn
12
SGD trajectory. We estimate the distribution of the self-normalized statistic by √ (b) Bn tn · φ θ̄tn − θ̄n ⋆,(b) 1 X b r ≤ x Fn (x; φ) = 1 Tn,φ := . Bn b=1 (b) Tr Σtn
(4)
n o b Let qb(δ; φ) = inf x : Fn (x; φ) ≥ 1 − δ . The resulting data-driven confidence region is ( b n (δ; φ) = R
θ ∈ Rd :
√
n · φ θ̄n − θ p Tr(Σn )
) ≤ qb(δ) .
(5)
A detailed procedure is described in Algorithm 1. It is worth noting the two stochastic gradients queried in each iteration may use different randomness. It is also valid to (b) (b) set ξk = ξk′ b , in which case the statistics θ̄tn , Σtn across sub-procedures still remain independent. Algorithm 1 records the full matrix Σn for notational clarity. In practice, when the confidence region is based only on Tr(Σn ), it is enough to update the trace directly. Theorem 3.2. Assuming continuous homogeneous functional |φ(·)| is sublinear, for Fbn (x; φ) in (4), if Bn → ∞, Bnn → 0 and x is a continuous point of the distribution function of Z ⋆ (φ), we have p Fbn (x; φ) → P (Z ⋆ (φ) ≤ x) .
If (1−δ)-quantile of Z ⋆ (φ) is a continuity point, the output confidence region in Algorithm 1 has asymptotic 1 − δ coverage when n → ∞, which is b n (δ; φ) → 1 − δ. P θ∗ ∈ R The continuity of the limiting distribution is an essential requirement for the asymptotic validity of the confidence region. If (1 − δ)-quantile of Z ⋆ (φ) is not a continuity point,
13
Algorithm 1 Data-driven Confidence Region for SGD Require: Fixed number of iterations n, confidence level δ, subsampling size tn = ⌊nr ⌋ (r ∈ (0, 1)), test function φ(·). (b) (b) (b) 1: Initialize b = 1, θ0 = θ0 = θ̄0 = θ̄0 = θ0 , and Σ0 = Σ0 = 0. 2: for k = 1, 2, · · · , n do 3: kb = k − (b − 1)tn . (b) 4: Query stochastic gradients g(θk−1 , ξk ) and g(θkb −1 , ξk′ b ). 5: Update the parameters: θk = θk−1 − ηk g(θk−1 , ξk ), (b)
(b)
(b)
θkb = θkb −1 − ηkb g(θkb −1 , ξk′ b ). 6:
Update the statistics: k θ̄k−1 + θk , k+1 (b) (b) kb θ̄kb −1 + θkb (b) θ̄kb = , kb + 1 (k − 1)Σk−1 + g(θk−1 , ξk )g(θk−1 , ξk )⊤ Σk = , k (b) (b) (b) (kb − 1)Σkb −1 + g(θkb −1 , ξk′ b )g(θkb −1 , ξk′ b )⊤ (b) . Σkb = kb θ̄k =
7: if kb = tn then 8: b = b + 1. (b) (b) (b) 9: θ0 = θ̄0 = θ0 and Σ0 = 0. 10: end if 11: end for 12: for i = 1, 2, · · · , b do ⋆,(i)
13: Calculate Tn,φ in Eqn (4). 14: end for bn in Eqn (4) and the quantile qb(δ; φ). 15: Calculate F
b n (δ; φ) in Eqn (5). 16: Output: Confidence region R one can still prove that the confidence region is conservative: b n (δ; φ) ≥ 1 − δ, lim P θ∗ ∈ R n→+∞
which guarantees that the resulting confidence region is conservative rather than overly aggressive. 14
Remark 3.2. It is worth noting that the joint convergence result in Theorem 3.1 continues to hold even when the stochastic gradient noise has finite variance [27]. In that case, the √ 1 coupling rate changes from the n1− α scaling to the classical n scaling, and the limiting distribution becomes Gaussian rather than stable. Nevertheless, the self-normalized statistic Tn⋆ (θ; φ) remains well-defined, and the inference procedure described in Algorithm 1 can be applied without modification. Importantly, the algorithm does not require prior knowledge of whether the stochastic gradient noise has finite or infinite variance, and therefore remains applicable across both regimes.
3.3
Discussion on Ellipsoidal Confidence Region
In constructing confidence regions for multi-dimensional parameters, several geometric forms can be considered. A popular choice in classical statistical inference is ellipsoidal confidence regions as they naturally account for correlations among coordinates through the covariance matrix. This idea underlies the well-known Hotelling t2 statistic and has also been adopted in recent work on statistical inference for stochastic gradient descent [16]. In the infinite variance regime considered in this paper, ellipsoidal confidence regions remain possible in principle. In particular, Corollary 3.1 shows that a self-normalized quadratic statistic can be constructed, which admits a non-degenerate limiting distribution provided that the limiting variance matrix W is invertible almost surely. Corollary 3.1. If W in Theorem 3.1 is invertible almost surely, the following limit theorem holds:
Tn† := n θ̄n − θ∗
⊤
d ⊤ −1 −1 −1 ∗ Σ−1 θ̄ − θ → L1 H W H L1 . n n
For Tn† in Corollary 3.1, the nuisance parameters α and b1 (·) are also canceled out on the left via self-normalization. However, the condition that W is invertible almost surely can only be satisfied in limited cases. In particular, one sufficient condition is that each coordinate of g(θ, ξ) has the same regularly varying tail index and the limit distribution 15
L1 does not lie in any subspace of Rd [21]. If some coordinates of g(θ, ξ) have lighter tails, the statistic Tn† will be invalid. To illustrate this issue, we consider the following simple quadratic optimization problem with dimension d = 2: min θ⊤ Hθ,
H :=
θ
1 1 −1 . 1 2
The stochastic gradient satisfies g(θ, ξ) = Hθ + ξ, where ξ = (ξ1 , ξ2 )⊤ and ξ1 ⊥ ⊥ ξ2 . Here ξ1 is an α1 -regularly varying random variable and ξ2 is an α2 -regularly varying random variable with 1 < α1 < α2 < 2. In this case, the joint central limit theorem holds in the following:
n
1− α1
1
θ̄n − θ
∗
1− α2
,n
1
d
Σn →
W1 0 , 0 0 , L L1 1
where the variance scaling term is singular and Tn† cannot be applied to establish a confidence region for θ∗ . In contrast, the confidence region based on the trace of secondmoment matrix does not require inversion and remains well-defined even when different coordinates exhibit heterogeneous tail indices.
4
Simulation Experiments
In this section, we evaluate the empirical performance of the proposed inference procedure through simulation studies. Our experiments are designed to examine the finite-sample behavior of the confidence regions constructed by Algorithm 1 under different stochastic gradient noise regimes, with particular focus on the empirical coverage probability and the size of the resulting confidence intervals. We consider two commonly used statistical models, linear regression and logistic regression discussed in [5]. For each setting, we run stochastic gradient descent with Polyak–Ruppert averaging and construct confidence regions using the proposed method. The nominal coverage level is set to 95% throughout all experiments. Moreover, we evaluate our methodology in a challenging setting where 16
the tail index varies across gradient queries, thereby assessing its robustness. The results are deferred to the Appendix E.1.
4.1
Linear Regression
For the linear regression model with heavy-tailed noise in [5]: ∗ yi = x⊤ i θ + εi ,
where {(xi , εi )}ni=1 are i.i.d. random variables satisfying xi ∈ Rd follows N (0, Σ) and εi follows a symmetric Pareto distribution with mean zero and index α. Moreover, xi ⊥ ⊥ εi . Given one pair sample (xi , yi ), the stochastic gradient is: g(θ, xi , yi ) = xi (x⊤ i θ − yi ). For the choice of Σ, we consider two cases: • Identity: Σ = I. • Toeplitz: Σij = q |i−j| with q = 0.3. For this model, we consider the sample size n = 106 and dimension d = 5, 20. The true parameter θ∗ is a d-dimensional vector and generated by the standard Gaussian distribution N (0, Id ). We consider three different choices of subsample size n0.6 , n0.7 , n0.8 . The nominal coverage probability is set to be 95% and we conduct 500 independent runs to calculate the empirical coverage rate. The performance of Algorithm 1 is measured by the average coverage rate of the confidence regions and the average length of the intervals. The results for linear regression are in Table 1. We summarize the Table 1 in the following. (a) Across both identity and Toeplitz covariance structures, the coverage rates are generally close to the nominal 95% target. (b) For dimension d = 5, the coverage rates are relatively insensitive to the choice 17
Figure 1: Linear Regression: The average coverage rate and length of confidence intervals for the cases α = 1.5 and d = 5 with different choices of subsample size. = 1.5, d = 5, Linear Regression.
1.0
= 1.5, d = 5, Linear Regression. 0.8
Confidence Interval Size
Coverage Rate
0.9
0.8
0.7
0.6
0.4
0.2
0.6
0.5
r=0.6 r=0.7 r=0.8
r=0.6 r=0.7 r=0.8 Baseline: 0.95 0.0
0.2
0.4 0.6 Number of SGD Iterations
0.8
1.0 1e6
0.0 0.0
0.2
0.4 0.6 Number of SGD Iterations
0.8
1.0 1e6
of subsample size when α = 1.5. In contrast, under Gaussian noise, coverage deteriorates when the subsample size is large (e.g., n0.8 ). (c) For dimension d = 20, the coverage rates become sensitive to the choice of subsample size. Among the three options considered, n0.8 yields the best performance, while the other two choices lead to conservative coverage. A plausible explanation is the limitation imposed by the total sample budget. (d) In general, coverage rates tend to decrease as the subsample size increases. This pattern is consistent with the usual bias–variance trade-off in subsampling: smaller blocks may introduce approximation bias, whereas larger blocks reduce the effective number of blocks available for quantile estimation. Hence, there should exist an intermediate “critical” subsample size that balances these errors under a fixed sample budget. q n) (e) For the average confidence region length qb Tr(Σ , we observe a decreasing trend as n the stochastic gradients become lighter-tailed. Theoretically, the confidence region 1 −1 e min{α,2} length is expected to scale as O n , which is consistent with our empirical findings.
18
Table 1: Linear Regression: The average coverage rate and length of confidence intervals, std 500 for the nominal coverage probability 95%. Standard errors √500 are reported in the brackets. Coverage Rate (%)
d n0.6
n0.7
Average Length (×10−2 ) n0.8
n0.6
n0.7
n0.8
94.4(1.0E-2)
11.33(1.7E-2)
12.57(2.5E-2)
15.47(4.0E-2)
98.2(5.9E-3)
93.8(1.1E-2)
94.26(8.0E-1)
46.84(3.1E-1)
17.00(2.7E-2)
95.8(9.0E-3) 95.6(9.2E-3)
91.0(1.3E-2)
1.91(7.1E-4)
1.91(8.8E-4)
1.92(1.3E-3)
97.8(6.6E-3)
91.8(1.2E-2)
3.48(7.7E-3)
2.98(4.5E-3)
2.35(1.0E-3)
94.6(1.0E-2) 94.4(1.0E-2)
86.4(1.5E-2)
0.26(4.5E-6)
0.25(8.1E-6)
0.23(1.2E-5)
98.0(6.3E-3)
95.6(9.2E-3)
0.50(1.4E-4)
0.48(1.5E-4)
0.47(2.4E-4)
97.2(7.4E-3) 96.6(8.1E-3)
94.0(1.1E-2)
12.53(1.9E-2)
13.70(2.8E-2)
16.82(4.7E-2)
99.0(4.4E-3)
95.8(9.0E-3)
106.4(8.9E-1)
60.94(4.4E-1)
18.65(2.9E-2)
96.0(8.8E-3) 95.2(9.6E-3)
91.4(1.3E-2)
2.10(8.3E-4)
2.07(9.9E-4)
2.09(1.5E-3)
98.2(5.9E-3)
93.6(1.1E-2)
3.99(8.0E-3)
3.54(6.0E-3)
2.68(1.0E-3)
94.6(1.0E-2) 92.6(1.2E-2)
87.0(1.5E-2)
0.29(4.8E-6)
0.27(8.9E-6)
0.25(1.4E-5)
95.8(9.0E-3)
0.64(2.2E-4)
0.62(2.4E-4)
0.62(4.6E-4)
Identity Σ α = 1.5
α = 1.8
Gaussian
5
97.0(7.6E-3) 96.8(7.9E-3)
20 99.8(2.0E-3) 5
20 99.6(2.8E-3) 5
20 98.6(5.3E-3)
Toeplitz Σ α = 1.5
α = 1.8
Gaussian
4.2
5
20 99.8(2.0E-3) 5
20 99.6(2.8E-3) 5
20 99.0(4.4E-3)
99.4(3.5E-3)
Logistic Regression
For the binary logistic regression model with heavy-tailed covariates in [5]: 1 , ∗ 1 + exp −x⊤ i β ∗ exp −x⊤ i β . P (yi = −1|xi ) = ∗ 1 + exp −x⊤ i β P (yi = 1|xi ) =
where {(xi , yi )}ni=1 are i.i.d. random variables satisfying xi ∈ Rd and yi follows a binary distribution. For the choice of distribution of xi , we consider two cases: (1)
(d)
• Homogeneous tail indices: xi , · · · , xi with index α.
19
are i.i.d symmetric Pareto distributions
Figure 2: Logistic Regression: The average coverage rate and length of confidence intervals for the cases α = 1.5 and d = 5 with different choices of subsample size. = 1.5, d = 5, Logistic Regression.
1.0
= 1.5, d = 5, Logistic Regression.
4
Confidence Interval Size
Coverage Rate
0.9
0.8
0.7
0.6
0.5
r=0.8 r=0.7 r=0.6
0.2
0.4 0.6 Number of SGD Iterations
0.8
(1)
2
1
r=0.8 r=0.7 r=0.6 Baseline: 0.95 0.0
3
1.0 1e6
0.0
(d)
• Heterogeneous tail indices: xi , · · · , xi
0.2
0.4 0.6 Number of SGD Iterations
0.8
1.0 1e6
are independently distributed symmetric
Pareto distributions with different indices α1 , · · · , αd with α1 = α, αd = 2.5 and α2 , · · · , αd−1 ∼ Uniform(α, 2). Given a pair of sample (xi , yi ), the stochastic gradient at θ is:
g(θ, xi , yi ) =
−yi exp(−yi x⊤ i θ) xi . ⊤ 1 + exp(−yi xi θ)
For this model, we consider the sample size n = 106 and dimension d = 5. The true parameter θ∗ is a d-dimensional vector generated by the standard Gaussian distribution N (0, Id ). We consider three different choices of subsample size n0.6 , n0.7 , n0.8 . The nominal coverage probability is set to be 95% and we conduct 500 independent runs to calculate the empirical coverage rate. The performance of Algorithm 1 is measured by the average coverage rate of the confidence regions and the average length of the intervals. The results for logistic regression are in Table 2. We summarize the Table 2 in the following. (a) Across both homogeneous and heterogeneous covariates, the coverage rates are generally close to the nominal 95% target. (b) Several observations are consistent with the linear regression results: (1) the coverage 20
Table 2: Logistic Regression: The average coverage rate and length of confidence intervals, for the nominal coverage probability 95%. Standard errors are reported in the brackets. Coverage Rate (%)
d n0.6
n0.7
Average Length (×10−2 ) n0.8
n0.6
n0.7
n0.8
Homogeneous xi α = 1.5
5
91.0(1.3E-2)
93.0(1.1E-2)
93.2(1.1E-2) 28.99(3.0E-2)
55.42(6.0E-2)
103.3(1.3E-1)
α = 1.8
5
98.2(5.9E-3)
98.0(6.3E-3)
96.4(8.3E-3)
8.73(1.2E-3)
9.17(3.8E-3)
15.05(1.9E-2)
Gaussian
5
99.8(2.0E-3)
97.4(7.1E-3) 90.4(1.3E-2)
0.92(2.0E-5)
0.74(2.9E-5)
0.63(3.7E-5)
Heterogeneous xi , 4 heaviest coordinates α = 1.5
5
92.4(1.2E-2)
94.2(1.0E-2)
93.4(1.1E-2) 25.10(2.9E-2)
46.14(5.9E-2)
87.89(1.2E-1)
α = 1.8
5
98.4(5.6E-3)
98.0(6.3E-3)
96.6(8.1E-3)
6.91(1.0E-3)
6.98(3.3E-3)
12.62(1.8E-2)
Heterogeneous xi , 3 heaviest coordinates α = 1.5
5
92.8(1.2E-2)
94.0(1.1E-2)
93.6(1.1E-2) 23.12(2.7E-2)
41.77(5.5E-2)
79.29(1.2E-1)
α = 1.8
5
98.4(5.6E-3)
98.2(5.9E-3)
97.0(7.6E-3)
6.72(9.8E-4)
6.74(3.2E-3)
12.16(1.7E-2)
Heterogeneous xi , 2 heaviest coordinates α = 1.5
5
94.0(1.1E-2)
94.2(1.0E-2)
93.6(1.1E-2) 19.10(2.6E-2)
31.53(5.3E-2)
59.42(1.1E-1)
α = 1.8
5
98.6(5.3E-3)
98.0(6.3E-3)
97.0(7.6E-3)
5.78(8.3E-4)
5.45(2.6E-3)
10.01(1.6E-2)
Heterogeneous xi , 1 heaviest coordinate α = 1.5
5
95.2(9.6E-3)
95.2(9.6E-3)
94.6(1.0E-2) 24.74(5.3E-2)
15.85(2.6E-2)
48.48(1.1E-1)
α = 1.8
5
98.6(5.3E-3)
98.0(6.3E-3)
97.0(7.6E-3)
4.88(2.3E-3)
9.13(1.5E-2)
5.37(7.7E-4)
rates are still relatively insensitive to the choice of subsample size when α = 1.5; (2) the coverage rates tend to decrease as the subsample size increases; (3) the average q n) confidence region length qb Tr(Σ decreases as the stochastic gradients become n lighter-tailed; (4) standard errors of the average confidence region length generally increase with the subsample size. (c) In the heterogeneous case, the coverage rates remain consistent with those observed in the homogeneous setting.
21
Figure 3: Tail diagnostics for the ℓ2 -norm of stochastic gradients in the AlexNet CIFAR-10 experiment. The stochastic gradients are evaluated near the reference parameter using mini-batches of size 64. Left: empirical histogram of the gradient norms. Right: Hill estimates of the tail index α. Hill Plot for Gradient Norms
Histogram of Gradient Norms 1600
4.0
1400
3.5
1000
Hill estimator
Count
1200
800 600
200
5
2.5 2.0
400
0
3.0
1.5 0.0
0.2
0.4
0.6 Gradient norm: l2
0.8
1.0
1.2
0
250
500
750 1000 1250 Number of upper order statistics k
1500
1750
2000
Application to CIFAR-10 dataset
In this section, we illustrate the proposed inference procedure on the CIFAR-10 dataset. This real-data experiment complements the simulation results in the previous section by considering an image-classification problem based on learned feature representations. The CIFAR-10 dataset is a standard benchmark for image classification. It consists of 60,000 color images of size 32 × 32 from 10 object categories, with 50,000 training images and 10,000 test images. We first train an AlexNet convolutional neural network on the CIFAR-10 training set using stochastic gradient descent. This network is designed for 32 × 32 inputs and uses a smaller fully connected representation, which is consistent with [36]. After training, we freeze the network and use the output of its penultimate layer, i.e., the layer before the final linear classifier, as the covariate vector. In our implementation, this produces a learned feature representation of dimension 448. We then construct a binary classification task by selecting one CIFAR-10 class as the positive class and grouping the remaining classes as the negative class, and fit a logistic regression model to the extracted features. In Fig. 3, we show tail diagnostics for the ℓ2 -norm of the stochastic gradient evaluated near the optimal parameter of AlexNet on CIFAR-10. The mini-batch size is set to 64, consistent with the experimental setup used in the subsequent inference results. The 22
Figure 4: Empirical reference coverage, coverage MSE, and fraction within 95% band in the CIFAR-10 experiment. Coverage Rate vs Subsampling Size r
MSE vs Subsampling Size r
target = 0.95
0.8
0.002 0.001
0
8
0.8
6
0.7
4
Subsampling size r
0.7
2
0.7
0
0.7
8
0.7
6
0.6
0
8
0.8
6
0.7
4
0.7
2
0.7
0
0.7
8
0.7
6
0.6
4
0.6
2
0.6
0
0.6
0.6
0
8
0.8
6
0.7
4
0.7
2
0.7
0
0.7
8
0.7
6
0.6
4
0.6
2
0.6
0
0.6
0.6
Subsampling size r
4
0.0
0.000 Subsampling size r
0.2
0.6
0.88
0.4
2
0.90
0.003
0.6
0.6
0.92
0.004
0
0.94
Fraction within 95% band
Mean Squared Error
0.005
0.6
0.96 Coverage rate
Fraction within 95% band vs Subsampling Size r
0.006
0.6
0.98
histogram on the left shows a highly right-skewed distribution with extreme values. The Hill plot is consistent with heavy-tailed behavior in this finite-sample experiment. In the Appendix E.3, we also consider a case in which the stochastic gradients exhibit lightertailed behavior under a different choice of feature mapping. The result is qualitatively similar to the results in this section. To evaluate the coverage performance, we conduct 500 independent runs. Since the true optimizer θ∗ is unknown in this real-data setting, we compute a highly accurate reference solution using the L-BFGS algorithm and use it as a proxy for θ∗ . For each independent run, we construct coordinate-wise 95% confidence intervals for all 448 parameters. For each coordinate, the empirical coverage rate is then computed as the fraction of the 500 confidence intervals that contain the corresponding coordinate of the L-BFGS reference solution. We provide additional details on the experimental setup and hyperparameter selection in Appendix E.2. Here, we summarize the results by reporting the average coverage rate across all coordinates. In addition, we report the mean squared error (MSE) of the coordinate-wise coverage rates relative to the nominal level 95%, defined as 448
2 1 Xc CRj − 0.95 , 448 j=1 c j denotes the empirical coverage rate for the jth coordinate. Meanwhile, where CR we report a diagnostic measure denoted by “Fraction within 95% band.” If the true c j based on coverage probability is exactly 0.95, then the Monte Carlo fluctuation of CR 23
Table 3: Results with reference coverage rate and mean squared error for 95% confidence intervals in the CIFAR-10 experiment. Average Coverage (%)
Mean Squared Error (×10−4 )
Fraction within 95% band (%)
Random Scaling [16]
92.05
10.43
20.54
Algorithm 1 (r = 0.7)
95.45
1.20
92.63
500 independent trials is approximately that of a binomial proportion with standard q . We therefore compute the proportion of parameters whose empirical error 0.95(1−0.95) 500 coverage rates fall within
r 0.95 ± 1.96
0.95(1 − 0.95) . 500
We first examine the sensitivity of the proposed method to the subsampling size parameter. Specifically, we vary r ∈ {0.60, 0.62, . . . , 0.78, 0.80}. Fig. 4 reports the average coverage rate and the coverage MSE under different choices of r. The average coverage rate decreases as r increases: smaller values of r tend to produce conservative confidence intervals, whereas larger values of r lead to under-coverage. The coverage MSE exhibits a U-shaped pattern. It decreases as r increases from 0.60 to around 0.70, and then increases for larger values of r. Similarly, the fraction within 95% also reaches the maximum at r = 0.7. This pattern reflects the usual bias-variance trade-off in subsampling: smaller subsamples may yield conservative inference, while overly large subsamples leave fewer effective blocks for stable quantile estimation. Based on Fig. 4, we select r = 0.7 for the proposed method in the following comparison. We compare Algorithm 1 with the random scaling method [16], which is effective for SGD in finite-variance settings. As shown in Table 3, random scaling leads to noticeable under-coverage for nominal 95% confidence intervals. In contrast, the proposed self-normalized subsampling method with r = 0.7 achieves an average coverage rate closer to the nominal level and yields a substantially smaller coverage MSE. The proposed method also attains a substantially larger fraction within this band, indicating that its coordinate-wise coverage rates are more consistently aligned with the nominal level. Overall, these results suggest that the self-normalized
24
subsampling procedure provides more reliable uncertainty quantification than the random scaling approach in this CIFAR-10 image-classification experiment.
6
Conclusion
In this paper, we develop a statistical inference methodology for SGD with infinitevariance stochastic gradients. Our approach is based on a joint central limit theorem involving Polyak–Ruppert averaged estimator and empirical second-moment normalizer of stochastic gradients. To address nuisance parameters in the limiting distribution, we construct a self-normalized statistic that cancels their influence. Building on the subsampling approach of Romano and Wolf [31], we further design an online procedure for estimating the quantiles of the limiting distribution and constructing confidence regions for the true parameter. We evaluate the method in synthetic experiments and in CIFAR-10 classification experiments. The results suggest that self-normalized subsampling can provide stable uncertainty quantification in settings where standard finite-variance methods may under-cover. Beyond these results, several important questions remain open. Progress on these questions would further broaden the scope of statistically principled uncertainty quantification for modern learning algorithms. Choice of the Subsampling Level. A key hyperparameter in our inference procedure is the subsampling level, which determines the block size used to approximate the distribution of the self-normalized statistic. This choice reflects a trade-off between CLT approximation bias and quantile estimation accuracy. The block size n/Bn should be large enough for the block-level statistic to be close to its limiting distribution, but not so large that the effective number of blocks becomes too small to estimate quantiles accurately. An optimal choice would therefore require explicit rates for both error terms. In particular, controlling the approximation bias would likely require Berry–Esseen-type bounds, or related non-asymptotic bounds, for the self-normalized statistic constructed from SGD iterates. We leave the development of finite-sample results as an important direction for 25
future research. Extensions to Other Optimizers. Although our analysis focuses on stochastic gradient descent, the main ideas underlying our inference framework are potentially applicable to a broader class of stochastic optimization algorithms, such as momentum and Adam. Extending our methodology to such algorithms would require understanding the appropriate central limit theorems and the corresponding normalization statistics. Once such results are established, the subsampling principle developed in this paper could potentially be used in a fully data-driven manner. Extensions to Non-convex Optimization. Another important direction is to extend the current theory beyond strongly convex objectives. Many optimization problems arising in modern machine learning are non-convex, where the algorithm may converge to different local minima. One possible route is to pursue a local analysis around an isolated local minimizer. If the iterate sequence can be shown to enter and remain in a neighborhood where the objective is approximately strongly convex, then the asymptotic framework developed here may still apply after localization [26]. Under such a regime, one might construct confidence regions for the local minimizer, conditional on convergence to that basin of attraction. Moreover, rare but large gradients in the heavy-tailed regime further complicate the non-convex setting, potentially causing repeated escapes from shallow basins or transitions between metastable regions [39]. Understanding the interplay among heavy tails, non-convex geometry, and uncertainty quantification remains largely open.
References [1] Krishna B Athreya. Bootstrap of the mean in the infinite variance case. The annals of statistics, pages 724–731, 1987. [2] Shuyang Bai, Murad S Taqqu, and Ting Zhang. A unified approach to self-normalized block sampling. Stochastic Processes and their Applications, 126(8):2465–2493, 2016. 26
[3] Melih Barsbey, Milad Sefidgaran, Murat A Erdogdu, Gael Richard, and Umut Simsekli. Heavy tails in sgd and compressibility of overparametrized neural networks. Advances in neural information processing systems, 34:29364–29378, 2021. [4] Dimitri P Bertsekas and John N Tsitsiklis. Gradient convergence in gradient methods with errors. SIAM Journal on Optimization, 10(3):627–642, 2000. [5] Jose Blanchet, Aleksandar Mijatović, and Wenhao Yang. Limit theorems for stochastic gradient descent with infinite variance. arXiv preprint arXiv:2410.16340, 2024. [6] Xi Chen, Jason D Lee, Xin T Tong, and Yichen Zhang. Statistical inference for model parameters in stochastic gradient descent. 2020. [7] Ashok Cutkosky and Harsh Mehta. High-probability bounds for non-convex stochastic optimization with heavy tails. Advances in Neural Information Processing Systems, 34:4883–4895, 2021. [8] Ewa Damek and Sebastian Mentemeier. Analysing heavy-tail properties of stochastic gradient descent by means of stochastic recurrence equations. arXiv preprint arXiv:2403.13868, 2024. [9] CA Goodsell and DL Hanson. Almost sure convergence for the robbins-monro process. The Annals of Probability, pages 890–901, 1976. [10] Mert Gurbuzbalaban, Umut Simsekli, and Lingjiong Zhu. The heavy-tail phenomenon in sgd. In International Conference on Machine Learning, pages 3964–3975. PMLR, 2021. [11] J Harold, G Kushner, and George Yin. Stochastic approximation and recursive algorithm and applications. Application of Mathematics, 35(10), 1997. [12] Liam Hodgkinson and Michael Mahoney. Multiplicative noise and heavy tails in stochastic optimization. In International Conference on Machine Learning, pages 4262–4274. PMLR, 2021. 27
[13] Zhezhe Jiao and Martin Keller-Ressel. Emergence of heavy tails in homogenized stochastic gradient descent. Advances in Neural Information Processing Systems, 37: 14066–14092, 2024. [14] Tatiana Pavlovna Krasulina. On stochastic approximation processes with infinite variance. Theory of Probability & Its Applications, 14(3):522–526, 1969. [15] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. Imagenet classification with deep convolutional neural networks. Advances in neural information processing systems, 25, 2012. [16] Sokbae Lee, Yuan Liao, Myung Hwan Seo, and Youngki Shin. Fast and robust online inference with stochastic gradient descent via random scaling. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, pages 7381–7389, 2022. [17] Gang Li. Almost sure convergence of stochastic approximation procedures. Statistica Sinica, pages 361–372, 1994. [18] Zijian Liu and Zhengyuan Zhou. Stochastic nonsmooth convex optimization with heavy-tailed noises: High-probability bound, in-expectation rate and initial distance adaptation. arXiv preprint arXiv:2303.12277, 2023. [19] Benjamin F Logan, CL Mallows, SO Rice, and Larry A Shepp. Limit distributions of self-normalized sums. The Annals of Probability, 1(5):788–809, 1973. [20] Michael Mahoney and Charles Martin. Traditional and heavy tailed self regularization in neural network models. In International Conference on Machine Learning, pages 4284–4293. PMLR, 2019. [21] Mark M Meerschaert and Hans-Peter Scheffler. Sample covariance matrix for random vectors with heavy tails. Journal of Theoretical Probability, 12(3):821–838, 1999. [22] Wenlong Mou, Koulik Khamaru, Martin J Wainwright, Peter L Bartlett, and Michael I
28
Jordan. Optimal variance-reduced stochastic approximation in banach spaces. arXiv preprint arXiv:2201.08518, 2022. [23] Eric Moulines and Francis Bach. Non-asymptotic analysis of stochastic approximation algorithms for machine learning. Advances in neural information processing systems, 24, 2011. [24] John P Nolan. Multivariate elliptically contoured stable distributions: theory and estimation. Computational statistics, 28(5):2067–2089, 2013. [25] Mariane Pelletier. On the almost sure asymptotic behaviour of stochastic algorithms. Stochastic processes and their applications, 78(2):217–244, 1998. [26] Mariane Pelletier. Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing. Annals of Applied Probability, pages 10–44, 1998. [27] Boris T Polyak and Anatoli B Juditsky. Acceleration of stochastic approximation by averaging. SIAM journal on control and optimization, 30(4):838–855, 1992. [28] Sidney I Resnick. Heavy-tail phenomena: probabilistic and statistical modeling. Springer Science & Business Media, 2007. [29] Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathematical statistics, pages 400–407, 1951. [30] Herbert Robbins and David Siegmund. A convergence theorem for non negative almost supermartingales and some applications. In Optimizing methods in statistics, pages 233–257. Elsevier, 1971. [31] Joseph P Romano and Michael Wolf. Subsampling inference for the mean in the heavy-tailed case. Metrika, 50(1):55–69, 1999. [32] Matt Rosenzweig. Lp spaces for 0< p< 1. 29
[33] Jerome Sacks. Asymptotic distribution of stochastic approximation procedures. The Annals of Mathematical Statistics, 29(2):373–405, 1958. [34] Adrien Schertzer and Loucas Pillaud-Vivien. Stochastic differential equations models for least-squares stochastic gradient descent. arXiv preprint arXiv:2407.02322, 2024. [35] Umut Simsekli, Mert Gürbüzbalaban, Thanh Huy Nguyen, Gaël Richard, and Levent Sagun. On the heavy-tailed theory of stochastic gradient descent for deep neural networks. arXiv preprint arXiv:1912.00018, 2019. [36] Umut Simsekli, Levent Sagun, and Mert Gurbuzbalaban. A tail-index analysis of stochastic gradient noise in deep neural networks. In International Conference on Machine Learning, pages 5827–5837. PMLR, 2019. [37] Umut Simsekli, Ozan Sener, George Deligiannidis, and Murat A Erdogdu. Hausdorff dimension, heavy tails, and generalization in neural networks. Advances in Neural Information Processing Systems, 33:5138–5151, 2020. [38] Hongjian Wang, Mert Gurbuzbalaban, Lingjiong Zhu, Umut Simsekli, and Murat A Erdogdu. Convergence rates of stochastic gradient descent under infinite noise variance. Advances in Neural Information Processing Systems, 34:18866–18877, 2021. [39] Xingyu Wang, Sewoong Oh, and Chang-Han Rhee. Eliminating sharp minima from sgd with truncated heavy-tailed noise. arXiv preprint arXiv:2102.04297, 2021. [40] Yanjie Zhong, Todd Kuffner, and Soumendra Lahiri. Online bootstrap inference with nonconvex stochastic gradient descent estimator. arXiv preprint arXiv:2306.02205, 2023. [41] Wanrong Zhu, Zhipeng Lou, Ziyang Wei, and Wei Biao Wu. High confidence level inference is almost free using parallel stochastic optimization. arXiv preprint arXiv:2401.09346, 2024.
30
[42] Yi Zhu and Jing Dong. On constructing confidence region for model parameters in stochastic gradient descent via batch means. In 2021 Winter Simulation Conference (WSC), pages 1–12. IEEE, 2021.
31
A
Generalized Central Limit Theorem (GCLT)
For i.i.d. random vectors X1 , · · · , Xn with E[∥X1 ∥] < +∞ and finite covariance E[X1 X1⊤ ] = Σ, classic CLT states the following weak convergence: n
1 X d √ (Xi − E[Xi ]) → N (0, Σ) . n i=1
(6)
In the setting considered in this paper, some components of X1 may not have finite variance. The weak convergence (6) is no longer correct. Instead, to establish limit theorem for infinite variance case, the following assumption should be made for X. Assumption A.1 (Multivariate regularly varying [28]). There exists α ∈ (1, 2), a slowly varying function b0 : R → R≥0 and a probability measure µ(·) on Sd−1 , such that, as r → +∞ ,
P (∥X∥ > r) = r−α b0 (r), X v ∈ · ∥X∥ > r → µ(·). P ∥X∥ And equivalently, there exists a Radon measure ν(·) on Rd \ {0} and a slowly varying function b1 (·), such that, as n → +∞: 1 v nP n− α b1 (n)X ∈ · → ν(·).
(7)
Assumption A.1 states that the distribution of X can be decomposed into a radial component and an angular component. The radial part ∥X∥ has an α-regularly varying tail. While the angular part X/∥X∥ converges vaguely to µ, which describes the limiting 1
direction of extremes. The scaling rate n− α b1 (n) in (7) determines the scaling rate in the generalized CLT. Theorem A.1 (Generalized CLT [28]). For i.i.d. random vectors X1 satisfy Assump-
32
tion A.1, the following weak convergence holds on DJ1 ([0, 1], Rd ):
n
1 −α
b1 (n)
⌊nt⌋ X
(Xi − E[Xi ]) ⇒ Lt ,
i=1
R where Lt is a Lévy process with characteristics (0, ν(·), γ) where γ = − ∥x∥>1 xν(dx). That is, for any ∥u∥ = 1, the characteristic function for Lt is: Z ⊤ iu⊤ x ⊤ E[exp (−iuLt )] = exp −itu γ + t e − 1 − iu x1(∥x∥ ≤ 1) ν(dx) . For the infinite-variance case, Theorem A.1 establishes a generalized central limit theorem that characterizes the weak convergence of normalized sums of i.i.d. random vectors. This generalized CLT depends on several structural parameters: the tail index α; the slowly varying function b1 (·); and the Radon measure ν(·). Together, these parameters determine the precise form of the limiting law. In contrast to the classical result (6), the scaling factor changes from n−1/2 to n−1/α b1 (n), and the limiting distribution is no longer Gaussian but rather a heavy-tailed distribution. In particular, when t = 1, the random vector L1 has heavy-tailed marginals, and its radial component ∥L1 ∥ follows an α-regularly varying distribution. One of the examples for L1 is the elliptically contoured stable distribution. Definition A.1 (Elliptically contoured stable distributions [24]). We denote a sub-class of stable distributions: α Sellipse (α) = F (·) EX∼F [exp(iu⊤ X)] = exp(−(u⊤ Σu) 2 + iu⊤ x), Σ ∈ Sd++ and x ∈ Rd .
α And we also define the scale function σα (Z; u) = u⊤ Σu 2 for distribution function of Z belonging to Sellipse (α). The elliptically contoured multivariate stable distribution is a symmetric subclass of the general multivariate α-stable family. When α = 2, it reduces to the multivariate 33
Gaussian distribution. The parameters x and Σ specify the location and the scale–shape structure of the distribution, respectively. It is important to note, however, that when α < 2, the matrix Σ should not be interpreted as a covariance or correlation matrix; instead, it is referred to as the scale matrix.
B
Confidence Region for Artificially Injected Noise
In the regime of PGD algorithm, the random noise ξart to the stochastic gradient is designed by practitioners, which implies the characterization of the random noise is aware. Practically, the choice of the injected noise is chosen to be i.i.d. symmetric Pareto distribution among the coordinates. Meanwhile, the intrinsic randomness may still exist, and we require the tail distribution of the intrinsic stochastic gradient is lighter than the injected noise. Assumption B.1. For the stochastic gradient g(θ, ξ) in Eqn (2), each coordinate of ξart ∈ Rd follows i.i.d. symmetric Pareto distribution with index α ∈ (1, 2) and scale λ > 0, i.e., for i ∈ [d], αλα d (i) P ξart ≤ t = . dt 2(λ + |t|)α+1 Moreover, lim tα P (∥∇g(θ, ξint )∥ > t) = 0 is required to dominate the intrinsic randomt→+∞
ness. In this case, the implicit slowly varying function in Theorem 2.2 is eliminated and the limit distribution is known except for the Hessian matrix H. Thus, if there is a consistent p bn → estimator for Hessian H H, we can simply plug it in the construction of confidence
region. Traditionally, we consider the following estimator of Hessian: n 1X b Hn = ∇g(θk , ξk+1 ). n k=1
And we also assume a smoothness assumption for the stochastic Hessian matrix. 34
(8)
Assumption B.2. We denote the stochastic Hessian matrix as H(θ, ξ) = ∇g(θ, ξ) and assume for any θ1 , θ2 ,
∥H(θ1 , ξ) − H(θ2 , ξ)∥ ≤ CH,Lip ∥θ1 − θ2 ∥ a.s. b n in Eqn (8) Theorem B.1. Under Assumptions 2.2, 2.3, 2.4, B.1, and B.2, estimator H p bn → is consistent, i.e. H H. Additionally, if there exists αH > 1 such that:
E ∥H(θ∗ , ξ) − H∥αH ≤ CH ,
then for any p ∈ (1, α), it holds:
bn − H ≲ E H
1 1− α1
n
H
1
+ n
ρ(1− p1 )
.
b n in Eqn (8) is a Under the smoothness assumption B.2, the Hessian estimator H consistent estimator for H. By Theorem 2.2 and Assumption B.1, we have: Pn
∗ i=1 (θi − θ )
n
1 α
d
→ H −1 L1 ,
(9)
where the characteristic function of H −1 L1 is given in the following proposition. Proposition B.1. Under Assumption B.1, the characteristic function of H −1 L1 in Eqn (9) is:
E exp −iu⊤ H −1 L1 = exp −σ∥H −1 u∥α . Thus, given consistent Hn−1 , the quantile of Hn−1 L can be calculated or simulated. The confidence region can be constructed: n o b † (δ; h) = θ ∈ Rd n1− α1 · h θ̄n − θ ≤ q † (δ; h) , R n n
35
(10)
where qn† (δ; h) = inf t P h(Hn−1 L) > t|Hn−1 ≤ δ . b † (δ; h) is asymptotically valid. It can be verified the confidence region R n b †n (δ; h) in Eqn (10), if h(H −1 L) has no atom, Theorem B.2. For the confidence region R we have b †n (δ; h) = 1 − δ. lim P θ∗ ∈ R
n→+∞
While this artificially injected noise approach provides a conceptually simple route for inference, it also has several limitations. First, it requires the injected noise to dominate the intrinsic stochastic gradient noise so that the asymptotic distribution is determined solely by the artificial perturbation. In practice, however, the tail index of the intrinsic stochastic gradient noise is typically unknown, making it difficult to guarantee such dominance. Second, when the injected noise dominates the learning dynamics, the resulting confidence regions primarily reflect the variability of the artificial perturbation rather than the intrinsic stochastic gradient noise. Consequently, the resulting confidence regions may become overly conservative and fail to capture useful information about the curvature of the loss landscape. These considerations motivate the development of an inference procedure that does not rely on artificially injected noise.
B.1
Proof of Theorem B.1
b n − H can be decomposed into: We notice the error H n n 1X 1X ∗ b Hn − H = (H(θ , ξk+1 ) − H) + (H(θk , ξk+1 ) − H(θ∗ , ξk+1 )) n k=1 n k=1
:= ∆1 + ∆2 .
36
By strong law of large numbers, we have: a.s.
∆1 → 0.
(11)
1X E∥∆2 ∥ ≤ E ∥H(θk , ξk+1 ) − H(θ∗ , ξk+1 )∥ n k=1
(12)
Besides, by Assumption B.2, we have: n
n
CH,Lip X E∥θk − θ∗ ∥ ≤ n k=1 p−1 = O n−ρ p , where p ∈ (1, α) and the last inequality is due to Lemma D.3. Therefore, we have: p
∆2 → 0.
(13)
Combining Eqn (11) and (13), we conclude: p bn → H H.
Additionally, if there exists αH ∈ (1, 2] such that: E ∥H(θ∗ , ξ) − H∥αH ≤ CH , the error ∆1 can be further controlled via Lemma D.1. That is, for any ∥u∥ = 1: α 1 E u⊤ ∆1 ≤ E u⊤ ∆1 H αH α1 α 1− H H 22−αH d 2 CH ≤ . 1− 1 n αH
37
(14)
Thus, combining Eqn (12) and (14), we conclude:
bn − H ≲ E H
B.2
1 1− α1 H
n
+
1 1
nρ(1− p )
.
Proof of Theorem B.2
We only need to prove qn† (δ; h) is consistent, i.e. p qn† (δ; h) → inf t P h H −1 L > t ≤ δ . In fact, for any matrix A, we denote
g(A, t; h) = P (|h (AL)| > t) . As h(H −1 L) has no atom, it implies g(·, t; h) is continuous at H −1 for any t > 0. By continuous mapping theorem, for any t > 0, we have: p
b −1 , t; h) → g(H −1 , t; h). g(H n Thus, consistency of qn† (δ; h) is obtained by inverse functional mapping.
C
Proof in Section 2 and Section 3
C.1
Proof of Theorem 2.2
We first revisit the iteration error: θn+1 − θ∗ = θn − θ∗ − ηn g(θn , ξn+1 ) = θn − θ∗ − ηn ∇ℓ(θn ) − ηn (g(θn , ξn+1 ) − ∇ℓ(θn )) = (I − ηn H)(θn − θ∗ ) − ηn (∇ℓ(θn ) − H(θn − θ∗ )) − ηn (g(θn , ξn+1 ) − ∇ℓ(θn )) = (I − ηn H)(θn − θ∗ ) − ηn Rn,1 − ηn Rn,2 − ηn δn , 38
where Rn,1 = ∇ℓ(θn ) − H(θn − θ∗ ), Rn,2 = g(θn , ξn+1 ) − ∇ℓ(θn ) − g(θ∗ , ξn+1 ), and δn = g(θ∗ , ξn+1 ). Thus, we have: θn+1 − θ∗ =
n Y
n n X Y
Ak (θ0 − θ∗ ) −
k=0
Ai ηk (Rk,1 + Rk,2 + δk ),
k=0 i=k+1
where Ak = I − ηk H. Then, the Polyak-Ruppert averaging estimator for θ satisfies: n n k−1 n k−1 k−1 1X 1 XX Y 1 XY ∗ ∗ Aj · (θ0 − θ ) − Ai · ηj (Rj,1 + Rj,2 + δj ) (θk − θ ) = n k=1 n k=1 j=0 n k=1 j=0 i=j+1
:= ∆1 (n) − ∆2 (n) − ∆3 (n) − ∆4 (n), where: n
k−1 k−1
n
k−1 k−1
n
k−1 k−1
1 XX Y ∆2 (n) = Ai · ηj δj , n k=1 j=0 i=j+1 ∆3 (n) =
1 XX Y Ai · ηj Rj,1 , n k=1 j=0 i=j+1
1 XX Y Ai · ηj Rj,2 . ∆4 (n) = n k=1 j=0 i=j+1
For the term ∆1 (n), as
j=0 (I − ηk H) exponentially decay, we have:
Qn−1
lim
ℓ1 (n) 1
n→+∞ n α −1
∥∆1 (n)∥ = 0.
Arranging terms in ∆2 (n), we have:
∆2 (n) =
n−1 n k−1 1X X Y Ai · ηj δj n j=0 k=j+1 i=j+1 n−1
1X n = A δj n j=0 j
39
n−1 n−1 1 X −1 1X n = H δj + (A − H −1 )δj , n j=0 n j=0 j
where Anj = ηj
Pn
k=j+1
n i=j+1 Ai . We notice that St :=
Qk−1
n −1 )δj is a martingale j=0 (Aj −H
Pt−1
w.r.t. t = 1, 2, · · · , n. By Lemma D.1 and β ∈ [1, α), we have: β
n−1
n−1
ℓ1 (n) X n E (Aj − H −1 )δj 1 n α j=0
≲
ℓ1 (n)β X n
β α
∥Anj − H −1 ∥β
j=0 n−1
≲
ℓ1 (n)β X n
β α
∥Anj − H −1 ∥.
j=0
Thus, as long as max{αρ, 1} < β < α, by Lemma D.2, we have: β
n−1
E
ℓ1 (n) X 1
nα
(Anj − H −1 )δj
≲
ℓ1 (n)β
−1 n ∥ n→+∞ j=0 ∥Aj − H
β α −ρ 2
n
j=0
Pn−1
n
β α +ρ 2
→ 0.
Thus, we have: ℓ1 (n) n
1 −1 α
∆2 (n) ⇒ H −1 L1 ,
R where Lt is a Lévy process on t ∈ [0, 1] with characteristics 0, ν(θ∗ , ·), − ∥x∥>1 xν(θ∗ , dx) . For the term ∆3 (n), we have: n−1
E∥∆3 (n)∥ = E
1X n A Rj,1 n j=0 j n−1
1X ≲ E∥Rj,1 ∥ n j=0 n−1
1X E∥θj − θ∗ ∥q n j=0 ε−ρq α−1 α ≲o n ,
≲
40
where the last step is due to Lemma D.3. Thus, we have: ℓ1 (n)
ε−(ρq−1) α−1 α E∥∆ (n)∥ = o ℓ (n)n . 3 1 1 −1
nα
p
As ρ > q −1 , we have above term tending to zero as n → +∞. Thus, ℓ11(n) ∆3 (n) → 0. For −1 nα
the term ∆4 (n), we write it as: n−1
∆4 (n) =
We also denote Stn =
n j=0 Aj Rj,2 ,
Pt−1
1X n A Rj,2 . n j=0 j
n which satisfies E[St+1 |σ({ξ1 , · · · , ξt })] = Stn . By
Lemma D.1, for β ∈ (1, α) we have: n−1
1 X E∥∆4 (n)∥ ≲ β E∥Rj,2 ∥β n j=0 β
n−1
1 X E∥θj − θ∗ ∥β ≲ β n j=0 n−1
1 X ε−ρβ α−1 α o j ≲ β n j=0 α−1 ≲ o nε−ρβ α −β+1 . Thus, we have: ℓ1 (n)β β
n α −β
β
β β ε−ρβ α−1 −α +1 α
E∥∆4 (n)∥ ≲ o ℓ1 (n) n p
α As long as 1+ρ(α−1) < β < α, we have ℓ11(n) ∆4 (n) → 0. −1 nα
C.2
Proof of Theorem 2.3
First, as Lt ∈ Sellipse (α), there exists Σ ∈ Sd++ such that: α σα (Lt ; u) = t · u⊤ Σu 2 . 41
.
d
Thus, the scale of Z =
R +∞ 0
e exp(−Ht)dL t satisfies:
σα (Z; u) =
Z +∞
⊤
e e u exp(−Ht)Σ exp(−Ht)u
α2
dt.
0
e and Σ commute, we assume H e = diag{e As H h1 , · · · , e hd } and Σ = diag{σ1 , · · · , σd }, where −1 e hi = hi − 1−αc > 0 and σi > 0. As α ∈ (1, 2), by Lemma D.5, we have
c
α−1
Z +∞
e e u⊤ exp(−Ht)Σ exp(−Ht)u
α2
dt
0
=cα−1
Z +∞ 0
≥cα−1
! α2 u2i σi exp(−2e hi t)
dt
i=1
d Z +∞ X i=1
=cα−1
d X
α2 ! α2 |ui |α σi exp(−αe hi t)dt α 2
0
d X u2i σi 2 (αe hi ) α
! α2 (15)
i=1
As f (t) = tα − αt + α − 1 ≥ 0 for t ≥ 0, we have: cα−1 1 ≥ α. hi αe hi
(16)
Apply (16) to (15), we have:
cα−1
Z +∞
e e u⊤ exp(−Ht)Σ exp(−Ht)u
α2
dt ≥
0
d X u2 σi i
i=1
Thus, our result is concluded.
42
h2i
! α2
α = u⊤ H −1 ΣH −1 u 2 .
C.3
Proof of Theorem 3.1
Lemma C.1. The empirical second-moment normalizer Σn satisfies the following decomposition: n
n
2 1− α
b1 (n) X
2
b1 (n) Σn =
n
2 α
g(θ∗ , ξi )g(θ∗ , ξi )⊤ + oP (1).
i=1
Proof. We also decompose g(θi , ξi ) into: g(θi , ξi ) = g(θ∗ , ξi ) + (g(θi , ξi ) − g(θ∗ , ξi )) := Xi + Yi . Then, we have: n n n n 1X 1X 1X 1X ⊤ ⊤ ⊤ Sn = X i Xi + Xi Yi + Yi Xi + Yi Yi⊤ n i=1 n i=1 n i=1 n i=1 n
1X Xi Xi⊤ + ∆n . := n i=1 2
For any u ∈ Sd−1 , multiplying n1− α ℓ1 (n−1 )2 both sides, we have: n
n
b1 (n)2 ⊤ 2b1 (n)2 X ⊤ b1 (n)2 X ⊤ 2 ⊤ u ∆ u = (u X )(u Y ) + (u Yi ) . n i i 2 2 2 n α −1 n α i=1 n α i=1 We notice that ∥Yi ∥ ≤ ∥θi − θ∗ ∥. Thus, for 1 < p < α we have v ! p1 u n n X X b1 (n) u b (n) 1 Et (u⊤ Yi )2 ≤ E ∥Yi ∥p 1 1 α nα n i=1 i=1 ! p1 n X b1 (n) ≤ E∥Yi ∥p 1 nα i=1 ! p1 n X b1 (n) ≤ ηip−1 1 nα i=1 ≤ b1 (n)n
43
1−(p−1)ρ 1 −α p
.
Thus, choosing p ∈
1+ρ α−1 +ρ
, α , we have: n
b1 (n)2 X n
2 α
p
(u⊤ Yi )2 → 0.
i=1
Similarly, we have: n
2b1 (n)2 X n
2 α
n
⊤
2b1 (n)2 X
⊤
E|u Xi ||u Yi | ≤
n
i=1
2 α
2b1 (n)
=
n
i=1 n 2 X
2 α
E∥Xi ∥∥θi − θ∗ ∥ E∥Xi ∥E∥θi − θ∗ ∥
i=1 n 2 X
2Cb1 (n)
≤
2
nα
1− p1
ηi
i=1 1
2
= b1 (n)2 n1−(1− p )ρ− α → 0. Thus, we have: n
2b1 (n)2 X n
2 α
p
(u⊤ Xi )(u⊤ Yi ) → 0.
i=1
2
2
Thus, for any u ∈ Sd−1 , we have b12(n) u⊤ ∆n u = oP (1), which implies b12(n) ∥∆n ∥ = −1 −1 nα
nα
oP (1). By Theorem 2.2, we notice that the dominate randomness for Polyak-averaging estimator is: 1 1− α
n
Pn b1 (n)
k=1 θk
n
−θ
∗
n
=
b1 (n) X n
1 α
H −1 g(θ∗ , ξk ) + oP (1).
k=1
By Lemma C.1, the dominate randomness for Σn is: n
2
n1− α b1 (n)2 Σn =
b1 (n)2 X n
2 α
g(θ∗ , ξk )g(θ∗ , ξk )⊤ + oP (1).
k=1
The, the joint convergence holds by Lemma D.4.
44
C.4
Proof of Theorem 3.2
Proof. We first introduce the infeasible version centered at the true parameter θ∗ : √ (b) ∗ tn h θ̄tn − θ (b) r Ten,h := , (b) Tr Σtn and define
Bn n o 1 X (b) e Fn (x; h) := 1 Ten,h ≤ x . Bn b=1
We first show that p Fen (x; h) → P (Z ⋆ (h) ≤ x)
for every continuity point x of Z ⋆ (h). Indeed, by Theorem 3.1 applied to each length-tn sub-procedure, (b) d Ten,h → Z ⋆ (h),
as n, tn → ∞. Therefore, for every continuity point x, (b) P Ten,h ≤ x → P (Z ⋆ (h) ≤ x) . Moreover, since the sub-procedures are independent, n o 1 1 (1) Var Fen (x; h) = Var 1 Ten,h ≤ x ≤ → 0. Bn Bn Hence p Fen (x; h) − EFen (x; h) → 0,
which proves p Fen (x; h) → P (Z ⋆ (h) ≤ x) .
45
It remains to show that replacing the infeasible centering θ∗ by θ̄n is asymptotically negligible. Define
√ (b) Dn,h :=
tn h θ̄n − θ∗ r . (b) Tr Σtn
By Theorem 3.1, 1
n1− α b1 (n)(θ̄n − θ∗ ) = Op (1), and 1− 2
(b) d
tn α b1 (tn )2 Σtn → W. Since P(Tr(W ) > 0) = 1, and since h is homogeneous and continuous, we obtain 1− 1
tn α b1 (tn )
(b) Dn,h = Op
1
n1− α b1 (n)
! = op (1),
1
where the last step follows from tn /n → 0 and the regular variation of n1− α b1 (n). Consequently, for every ε > 0, P
(b) Dn,h > ε
→ 0.
As |h(·)| is sublinear, then (b) (b) (b) Tbn,h − Ten,h ≤ Dn,h .
Thus, for every ε > 0, n o n o n o (b) (b) (b) 1 Ten,h ≤ x − ε − 1 Dn,h > ε ≤ 1 Tbn,h ≤ x and o n o n o n (b) (b) (b) 1 Tbn,h ≤ x ≤ 1 Ten,h ≤ x + ε + 1 Dn,h > ε . Averaging over b = 1, . . . , Bn , we obtain
Fen (x − ε; h) − An (ε) ≤ Fbn (x; h) ≤ Fen (x + ε; h) + An (ε), 46
where An (ε) :=
Bn n o 1 X (b) 1 Dn,h > ε . Bn b=1
Since (1) EAn (ε) = P Dn,h > ε → 0, Markov’s inequality implies p
An (ε) → 0. Therefore, Fen (x − ε; h) + op (1) ≤ Fbn (x; h) ≤ Fen (x + ε; h) + op (1). Letting n → ∞ and then ε ↓ 0, and using the continuity of the distribution of Z ⋆ (h) at x, we conclude that p Fbn (x; h) → P (Z ⋆ (h) ≤ x) .
Finally, let q(δ; h) := inf {x : P (Z ⋆ (h) ≤ x) ≥ 1 − δ} , and n o qb(δ; h) := inf x : Fbn (x; h) ≥ 1 − δ . If the distribution function of Z ⋆ (h) is continuous at q(δ; h), the preceding convergence implies p
qb(δ; h) → q(δ; h). Moreover, by Theorem 3.1, d
|Tn⋆ (θ∗ ; h)| → Z ⋆ (h). Hence, by Slutsky’s theorem, bn (δ; h) = P (|Tn⋆ (θ∗ ; h)| ≤ qb(δ; h)) P θ∗ ∈ R → P (Z ⋆ (h) ≤ q(δ; h)) = 1 − δ. 47
This proves the desired asymptotic coverage.
D
Auxiliary Lemmas
Lemma D.1 (Lemma 7 in [38]). Suppose p ∈ [1, 2] and {St } is a d-dimensional martingale with finite p moment for every t, and S0 = 0. Then: t X
2−p 1− p2
p
E∥St ∥ ≤ 2
d
E∥Si − Si−1 ∥p .
i=1
Lemma D.2. Suppose H is symmetric positive definite matrix and ηn = c · n−ρ , where P Q ρ ∈ (0, 1). We denote Ak = I − ηk H and Anj = ηj nk=j+1 k−1 i=j+1 Ai . Then we have: 1. (Lemma 1 in [27]) sup0≤j≤n−1 ∥Anj ∥ ≤ CA . 2. (Lemma 17 in [38]) for any ρ < κ ≤ 1, Pn−1
n −1 ∥ j=0 ∥Aj − H = 0. nκ
lim
n→+∞
Lemma D.3 (Corollary 4 in [38]). For any ε > 0 and q ∈ (1, α), then: α−1 E∥θn − θ∗ ∥q = o nε−ρq α . Lemma D.4 (Lemma 5 in [21]). For an i.i.d. sequence α-regularly varying random vector {Xi }ni=1 , the joint weak convergence holds: n
ℓ1 (n) X n
1 α
i=1
!
n
Xi ,
ℓ1 (n)2 X n
2 α
Xi Xi⊤
d
→ (Z, W ).
i=1
Lemma D.5 (Proposition 4 in [32]). For p ∈ (0, 1),
R X
|f |p dµ < +∞ and
R
we have: Z X
p1 Z p1 Z p1 p p ||f | + |g|| dµ ≥ |f | dµ + |g| dµ . p
X
48
X
X
|g|p dµ < +∞,
E
More Experimental Results
E.1
Simulation: Robustness to Varying Tail Indices
The theoretical analysis in this paper assumes that the stochastic gradient noise follows a heavy-tailed distribution with a fixed tail index α ∈ (1, 2). This assumption enables us to characterize the asymptotic behavior of the Polyak–Ruppert averaged SGD estimator via stable limit theorems. However, in practical scenarios, the heaviness of stochastic gradient noise may vary over time. For instance, the distribution of stochastic gradients may change across iterations due to different sampled mini-batches or evolving model parameters. As a result, the tail indices of the stochastic gradient noises may not remain constant throughout the optimization process. To evaluate the robustness of our inference procedure under such heterogeneous noise regimes, we consider a setting where the tail index of the stochastic gradient noise varies randomly across iterations. Specifically, at iteration k, the stochastic gradient noise follows a Pareto distribution with tail index αk , where αk is independently drawn from a uniform distribution. In the experiment, we consider three different ranges for αk : [1.5, 1.8), [1.8, 2.2), [2.2, 4), corresponding respectively to a heavy-tailed regime, a mixed heavy- and light-tailed regime, and a light-tailed regime. This setting departs from the theoretical framework developed in earlier sections, where the tail index is assumed to be fixed. We deploy this varying index setting to both linear regression and logistic regression. For logistic regression, we also consider both homogeneous and heterogeneous settings. In the homogeneous setting, a single αk is sampled at each iterate and applied to all coordinates, whereas in the heterogeneous setting, different αk values are sampled for different coordinates. Nevertheless, it is plausible that the asymptotic behavior of the averaged SGD estimator remains governed by the heaviest tail in the system. A rigorous theoretical analysis of such varying tail regimes is beyond the scope of the present work and is left for future research. Empirically, however, we observe that the proposed inference procedure continues to 49
Table 4: Varying Indices: The average coverage rate and length of confidence intervals, for the nominal coverage probability 95%. Standard errors are reported in the brackets. Coverage Rate (%) n0.6
Average Length (×10−2 )
n0.7
n0.8
n0.6
n0.7
n0.8
Linear Regression (Identity Σ) α ∈ (1.5, 1.8)
96.2(8.6E-3)
94.6(1.0E-2)
93.0(1.1E-2)
4.828(3.6E-3)
5.161(4.9E-3)
5.657(7.3E-3)
α ∈ (1.8, 2.2)
93.8(1.1E-2)
92.4(1.2E-2)
87.6(1.5E-2)
1.136(2.3E-4)
1.119(2.6E-4)
1.065(2.9E-4)
α ∈ (2.2, 4)
94.6(1.0E-2)
94.4(1.0E-2)
86.4(1.5E-2) 2.644(4.5E-6)
2.529(8.1E-6)
2.327(1.3E-5)
Linear Regression (Toeplitz Σ) α ∈ (1.5, 1.8)
96.6(8.1E-3)
94.6(1.0E-2)
91.8(1.2E-2)
3.940(2.4E-3)
3.892(1.8E-3)
4.083(3.2E-3)
α ∈ (1.8, 2.2)
94.2(1.0E-2)
90.6(1.3E-2)
87.8(1.5E-2)
1.248(2.5E-4)
1.217(2.9E-4)
1.156(3.4E-4)
α ∈ (2.2, 4)
94.6(1.0E-2)
92.6(1.2E-2)
87.0(1.5E-2) 2.877(4.8E-6)
2.722(8.9E-6)
2.494(1.4E-5)
Logistic Regression (Homogeneous xi ) α ∈ (1.5, 1.8)
94.0(1.1E-2)
94.6(1.0E-2)
94.2(1.0E-2)
19.03(2.0E-2)
31.25(4.0E-2)
65.88(8.6E-2)
α ∈ (1.8, 2.2)
99.6(2.8E-3)
99.2(4.0E-3)
97.8(6.6E-3)
5.012(1.1E-3)
4.335(2.5E-3)
8.688(1.2E-2)
α ∈ (2.2, 4)
99.8(2.0E-3)
99.2(4.0E-3)
94.2(1.0E-2) 2.016(5.8E-5)
1.446(9.9E-5)
1.211(1.0E-3)
Logistic Regression (Heterogeneous xi ) α ∈ (1.5, 1.8)
93.6(1.1E-2)
95.2(9.6E-3)
94.6(1.0E-2)
18.25(1.3E-2)
31.65(2.7E-2)
55.61(6.1E-2)
α ∈ (1.8, 2.2)
99.2(4.0E-3)
98.6(5.3E-3)
97.0(7.6E-3)
5.000(8.1E-4)
4.312(1.7E-3)
8.271(1.5E-2)
α ∈ (2.2, 4)
100.0(0E-0)
99.8(2.0E-3)
94.4(1.0E-2)
2.025(6.5E-5)
1.465(1.0E-4)
1.498(2.8E-3)
achieve reliable coverage probabilities in both linear regression and logistic regression models in Table 4, suggesting that the method remains robust even when the tail behavior of stochastic gradients varies across iterations. We summarize the Table 4 in the following. (a) Across all models considered, including linear regression with identity and Toeplitz covariance structures and logistic regression with homogeneous and heterogeneous covariates, the empirical coverage rates are in general reasonably close to the nominal 95% level. This suggests that the proposed inference procedure remains robust even when the tail index varies randomly across iterations. (b) In the linear regression experiments, the coverage rates are fairly stable when the 50
tail indices are relatively heavy, namely α ∈ (1.5, 1.8). As the noise becomes lightertailed, however, the coverage rates tend to deteriorate for larger subsample sizes, especially when the subsample size is chosen as n0.8 . (c) In the logistic regression experiments, the procedure often yields conservative coverage when the tail indices are in the lighter-tailed regimes α ∈ (1.8, 2.2) and α ∈ (2.2, 4), with empirical coverage rates frequently exceeding 99%. This phenomenon is more pronounced for smaller subsample sizes, indicating that the method may overestimate uncertainty in these settings. (d) The effect of the subsample size is model dependent. For linear regression, the coverage rates generally decrease as the subsample size increases in the range r ∈ (0.6, 0.8), especially under lighter-tailed noises. In contrast, for logistic regression with the heaviest-tailed regime α ∈ (1.5, 1.8), the coverage rates remain relatively stable across different subsample sizes. Overall, these results again reflect a biasvariance tradeoff in the choice of subsample size. (e) The average confidence interval length decreases substantially as the stochastic gradient noise becomes lighter-tailed. This trend is consistent across both linear and logistic regression models, and agrees with the theoretical intuition that heavier tails lead to slower concentration and hence wider confidence intervals.
E.2
CIFAR-10: Experimental Setup
Pretraining. We use a CIFAR-10 adapted AlexNet architecture following Simsekli et al. [36]. The network is trained on the CIFAR-10 training set using stochastic gradient descent with mini-batch size 128, learning rate 0.1, and 100 epochs. After pretraining, we freeze the network and use the output of the penultimate layer as the feature vector. In our implementation, this gives a 448-dimensional representation for each image. We then normalize the extracted features coordinate-wise using the empirical mean and standard deviation computed from the training data. In Section E.3, we further report an additional 51
Figure 5: Tail diagnostics for the ℓ2 -norm of stochastic gradients in the full AlexNet CIFAR-10 experiment. The stochastic gradients are evaluated near the reference parameter using mini-batches of size 64. Left: empirical histogram of the gradient norms. Right: Hill estimates of the tail index α. Histogram of Gradient Norms 300
16
250
14 12 Hill estimator
Count
200 150
10 8
100
6
50
4
0
Hill Plot for Gradient Norms
18
2
2.0
2.5
3.0
3.5
4.0 Gradient Norm: l2
4.5
5.0
0
5.5
200
400 600 Number of upper order statistics k
800
1000
experiment based on the full AlexNet architecture of Krizhevsky et al. [15], which produces a 4096-dimensional penultimate-layer feature representation. Subsampling. Using the normalized AlexNet features, we construct a binary classification task by selecting one CIFAR-10 class as the positive class and grouping the remaining classes as the negative class. We then fit a regularized logistic regression model on top of the extracted features. No intercept term is included, so the parameter dimension is exactly 448. Specifically, we consider N 1 X λ ⊤ min log 1 + exp(z θ) − yi zi⊤ θ + ∥θ∥22 , i 448 θ∈R N i=1 2
where zi denotes the normalized feature vector and yi ∈ {0, 1} denotes the binary label. We set the ℓ2 regularization coefficient to λ = 0.1. Since the true optimizer θ∗ is unknown in this real-data experiment, we compute a high-accuracy reference solution using L-BFGS and use it as a proxy for θ∗ . The L-BFGS optimization tolerance is set to 10−10 . For the SGD trajectories used in the inference procedure, we use mini-batches of size 64 and the step-size schedule ηt = 0.5 t−0.6 . Each trajectory is run for 106 iterations, and Polyak–Ruppert averaging is applied to the iterates.
52
Figure 6: Empirical reference coverage, coverage MSE, and fraction within 95% band in the CIFAR-10 experiment with full AlexNet. Coverage Rate vs Subsampling Size r
MSE vs Subsampling Size r
Fraction within 95% band vs Subsampling Size r
0.001
0
8
0.8
6
0.7
4
Subsampling size r
0.7
2
0.7
0
0.7
8
0.7
6
0.6
0
8
0.8
6
0.7
4
0.7
2
0.7
0
0.7
8
0.7
6
0.6
4
0.6
2
0.6
0
0.6
0.6
0
8
0.8
6
0.7
4
0.7
2
0.7
0
0.7
8
0.7
6
0.6
4
0.6
2
0.6
0
0.6
0.6
Subsampling size r
4
0.0
0.000 Subsampling size r
0.2
0.6
0.88
0.4
2
0.90
0.002
0.6
0.6
0.92
0.003
0
0.004
0.94
0.6
0.96
0.8 Fraction within 95% band
0.005
0.6
0.98
Mean Squared Error
Coverage rate
target = 0.95
Random Scaling. As a benchmark, we compare the proposed method with the random scaling method of Lee et al. [16]. The same regularized logistic regression model is used, again without an intercept term. The training hyperparameters are the same as those used for the subsampling approach; only the critical value differs. In Lee et al. [16], for an SGD trajectory {θt }nt=1 , the random-scaling variance estimator is n 1 X 2 s (θ̄s − θ̄n )(θ̄s − θ̄n )⊤ . VbnRS = 2 n s=1
For coordinate-wise nominal 95% confidence intervals, we use s RS θ̄n,j ± q0.975
RS Vbn,jj , n
RS where q0.975 = 6.747 is the tabulated random-scaling critical value used in [16].
E.3
CIFAR-10: Full AlexNet
We further conduct an additional CIFAR-10 experiment using the original AlexNet architecture of Krizhevsky et al. [15]. This experiment is intended as a robustness check for the small AlexNet experiment reported in Section 5. The experimental protocol is the same as that described in the preceding subsection, except that the feature extractor is replaced by the original AlexNet model. In Fig. 5, we observe that, with the larger feature dimension, the gradient-norm diagnostics appear more consistent with a lighter-tailed, 53
Table 5: Coverage rate, mean squared error, and fraction within 95% band for 95%confidence intervals in CIFAR-10 dataset with full AlexNet. Average Coverage (%)
mean squared error (×10−4 )
Fraction within 95% band (%)
Random Scaling [16]
93.07
5.90
46.97
Algorithm 1 (r = 0.72)
94.81
1.76
88.53
possibly finite-variance, regime. Nevertheless, our methodology remains directly applicable. The corresponding inference results are reported below. Fig. 6 is the performance of the proposed method for different subsampling levels r ∈ {0.60, 0.62, . . . , 0.78, 0.80}. The average coverage rate decreases as r increases: smaller values of r tend to produce conservative intervals, whereas larger values of r lead to undercoverage. The coverage MSE exhibits a U-shaped pattern, decreasing as r increases from 0.60 to around 0.72 and then increasing for larger values of r. This pattern is consistent with the bias–variance trade-off in subsampling: smaller blocks may yield conservative inference, whereas larger blocks reduce the effective number of blocks available for estimating the critical value. Based on the coverage MSE, we select r = 0.72 for the comparison below. Table 5 compares the proposed method with random scaling. Random scaling leads to noticeable under-coverage, with an average coverage rate of 93.07% for nominal 95% confidence intervals. In contrast, Algorithm 1 with r = 0.72 achieves an average coverage rate of 94.81%, much closer to the nominal level. The coverage MSE is also reduced from 5.90 × 10−4 under random scaling to 1.76 × 10−4 under the proposed method. Moreover, the fraction within the nominal 95% Monte Carlo band increases from 46.97% to 88.53%.
54