Proceedings of the 2023 Winter Simulation Conference
C. G. Corlu, S. R. Hunter, H. Lam, B. S. Onggo, J. Shortle, and B. Biller, eds.

RESAMPLING STOCHASTIC GRADIENT DESCENT CHEAPLYHenry Lam
Zitong Wang
Industrial Engineering andOperations Research
Columbia University
500 West 120th Street
New York, NY 10027, USA

ABSTRACT
Stochastic gradient descent (SGD) or stochastic approximation has been widely used in model training andstochastic optimization. While there is a huge literature on analyzing its convergence, inference on theobtained solutions from SGD has only been recently studied, yet is important due to the growing need for
uncertainty quantification. We investigate two easily implementable resampling-based methods to constructconfidence intervals for SGD solutions. One uses multiple, but few, SGDs in parallel via resampling withreplacement from the data, and another operates this in an online fashion. Our methods can be regarded
as enhancements of established bootstrap schemes to substantially reduce the computation effort in termsof resampling requirements, while at the same time bypasses the intricate mixing conditions in existingbatching methods. We achieve these via a recent cheap bootstrap idea and Berry-Esseen-type bound for
SGD.
1

INTRODUCTION

Stochastic optimization commonly arises in many applications across machine learning, operations research,and scientific analysis. The problem can be formulated as:
min H(x) ≜ Eζ ∼P [h(x, ζ )]

x∈Rd

(1)

where P is some underlying data distribution governing the randomness ζ ∈ Ω, and h is some known
real-valued function. Stochastic gradient descent (SGD) or stochastic approximation is a popular numericalapproach to solve (1). With an initial guess x0 ∈ Rd , SGD iteratively updates the solution using
xt+1 = xt − ηt ∇h(xt , ζt+1 ),t = 0, . . . , n − 1(2)where ζt is a sample drawn using a Monte Carlo model generator or real data. While the Robbins-Monro(Robbins and Monro 1951) procedure outputs xn after a large number of iterations, one could also take thenaverage x̄n ≜ n1 ∑t=1
xt as the output, which is known as the Polyak-Ruppert-Juditsky averaging (Polyak andJuditsky 1992) or for convenience in this paper we call it averaged stochastic gradient descent (ASGD).
Both approaches are common, with ASGD known to be more robust with respect to the step size ηt (Rakhlin
et al. 2012).
We are interested in conducting inference, or quantifying the uncertainty, in SGD. That is, we aim
to construct, using the iterates (2), a 1 − γ confidence interval for (each component of) the true optimal
solution x∗ of problem (1). Despite the popularity of SGD, to our best knowledge, this problem has been

Lam and Wang
systematically studied only recently, driven by applications in exploration (Lattimore and Szepesvári 2020)
and as stopping criteria (Su and Zhu 2018; Fang et al. 2018; Chen et al. 2020).
1.1 Existing Methods and Challenges
A main challenge in SGD inference is the serial dependence incurred by the sequence {xt }, which makes
the construction of a consistent standard error estimator intricate. We discuss the ideas and challenges
of several recent works that address this issue. Chen et al. (2020) proposed two methods, one based on
the delta method that directly approximates the asymptotic covariance of gradient ∇h(x, ζ ) and Hessian
at optimal ∇2 H(x∗ ). However, this requires the computation and storage of the Hessian matrix, which
can be computationally demanding. For example, backpropagation can only provide first-order gradient
information (Rumelhart et al. 1986), and arguably, a major advantage of SGD lies in its Hessian-freenature. In addition, storing a Hessian matrix requires an expensive O(d 2 ) space. These put aside the subtle
regularity assumptions needed for consistency as noted by Chen et al. (2020) themselves.
Motivated by these, Chen et al. (2020)’s second method borrows the batch-means idea in stochastic
simulation output analysis (Glynn and Iglehart 1990; Schmeiser 1982; Schruben 1983; Glynn and Lam2018) and Markov Chain Monte Carlo (Geyer 1992; Flegal and Jones 2010; Jones et al. 2006). Thisapproach divides the iterations of SGD into M batches of increasing sizes and aggregates the means of thesebatches to construct confidence intervals. However, the batch-mean method introduces hyperparameter M,
the number of batches, to tune. Additionally, experiments show that this method is more sensitive to thequality of converges of SGD and could underperform other methods. Relatedly, Li et al. (2018) presented
a batch-means method for inference in M-estimation by using SGD trajectory with a constant step size.
Instead of using batches with increasing lengths, they use batches with a fixed length but separated bygaps to overcome the dependence between iterations of SGD. Zhu and Dong (2021) studied a batch-mean
algorithm to construct a d-dimensional confidence region for the optimal solution to problem (1). Theirmethod works by canceling out the asymptotic covariance matrix of the rescaled residue of SGD using an
F-type statistic. Additionally, they explore the influence of varying the number and sizes of batches on theefficacy of their algorithm.
Another approach is to use the bootstrap which, advantageously, does not succumb to the computationload of variance estimation nor the tuning and sensitivity of batch sizes. Fang et al. (2018) developed
an online bootstrap method that persistently maintains B perturbed version of SGD estimates, updated
upon each data arrival. However, B is required to be large for their method to work properly. For linear
regression problems of dimensions 10 or 20, they set B = 200, which means 200 times more computational
cost compared to running the SGD itself or using batch means.
Yet another approach called HiGrad is proposed by Su and Zhu (2018) based on “splitting” an SGD
trajectory. More specifically, it first runs SGD for several steps. Then, with the outcome of this thread as astarting point, perform multiple SGD threads using different data, and continue this branching process foreach thread’s outcome until data are used up. A confidence interval is constructed using all the obtained
split outcomes. HiGrad requires a substantial modification to the original SGD runs; in fact, there is no
more “original" run of SGD in HiGrad.
Finally, we briefly mention a line of work on quantifying algorithmic randomness, including Lundeet al. (2021) that applied the bootstrap on streaming principal component analysis (Oja 1982), and Chen
and Lopes (2020) on randomized Newton methods. Moreover, Nesterov and Vial (2008) gave a complexity
bound on the number of iterations of their method in relation to the confidence level on reaching the optimalvalue via SGD. However, all these works focus on assessing the uncertainty from algorithmic randomness
and treat the data as fixed, and hence are less relevant to our focus in this paper.

Lam and Wang
1.2 Our Contributions
Our discussion above reveals that existing approaches in SGD inference encounter either intricate algorithmic
tuning that relates to mixing conditions (batching), substantial modification on the SGD itself (HiGrad),or computation and storage challenges (delta method and online bootstrap). In this paper, we investigate amethodology that resolves these three challenges simultaneously. More precisely, we adopt the bootstrap
approach, which does not require mixing-related tuning nor substantial modification to the original SGD.
At the same time, we enhance the bootstrap to make it substantially lighter in terms of resampling cost.
The latter is made possible by using a recent “cheap bootstrap" idea (Lam 2022b; Lam and Liu 2023; Lam
2022a).
Our methodology can be implemented in both offline and online fashions. The offline version, which
we call the Cheap Offline Bootstrap (COfB), reruns the SGD using resampling with replacement from thedata B times and constructs confidence intervals from these resampled iterates via an approach similar tothe standard error bootstrap. However, while this approach may appear to require heavy resampling effort,
our key assertion is that the B in our implementation can be very small (such as 3). In this way, ourapproach is computationally less demanding than the delta method (Chen et al. 2020) and online bootstrap(Fang et al. 2018), does not require hyperparameter tuning in batch-means (Chen et al. 2020; Zhu andDong 2021), and also does not substantially modify the SGD trajectory in HiGrad (Su and Zhu 2018).
A caveat of COfB is that we can only rerun SGD after all the data become available. Thus, it cannot
be used in a single-pass streaming fashion. To address this, our online version, Cheap Online Bootstrap(COnB), runs multiple (B + 1) SGDs in parallel on the fly as new data comes in. COnB borrows the idea
of Fang et al. (2018) in perturbing the gradient estimate in the SGD iteration. However, like COfB, it iscomputationally much cheaper than Fang et al. (2018) as it only needs to maintain a very small numberof SGD runs. In both our theory and experimentation, we illustrate that using B = 3 already producesconsistently better coverage than the existing approaches.
Our methodology synthesizes two recent ideas. One, as mentioned earlier, is the recent cheap bootstrap
idea. Roughly speaking, instead of using the resemblance of the resample distribution to the sampling
distribution as in classical bootstraps, the cheap bootstrap exploits the approximate independence between
the resample and original estimates. Coupled with asymptotic normality, this allows to construct a pivotal
statistic with an extremely small number (as low as 1) of resample runs B. While the cheap bootstrap
gives rise to asymptotically exact-coverage intervals, it also comes with the cost of (arguably fair) larger
interval lengths when B is small. Nonetheless, as discussed in Lam (2022b), Lam and Liu (2023), theinterval length advantageously shrinks quickly as B increases away from 1. Our second methodologicalelement, which constitutes our main technical development, is to show the asymptotic independence, more
precisely a joint central limit theorem, for the original and the resampled SGD runs under resampling withreplacement, and how to suitably aggregate the outputs guided by this theorem using the cheap bootstrap
logic. To attain this independence, we generalize the recent non-asymptotic bounds for ASGD studied by
Shao and Zhang (2022), Anastasiou et al. (2019) to hold uniformly for both the original and resampled
runs, under SGD and ASGD settings.
2

METHODOLOGY

Denote the underlying data distribution by P. Let xout be the output of (A)SGD, using step sizes ηt = ηt −αnnxt , and in SGD xout = xn ,drew from P. More precisely, in ASGD xout = 1n ∑t=1
and i.i.d. data {ζt }t=1where xt is the solution obtained in the t-th iteration of (2). Let P̂n denote the empirical distribution fromn , i.e., P̂ (·) = 1 n I(ζ ∈ ·), where I(·) denotes the indicator function. We also use (·) todata {ζt }t=1ntin ∑t=1denote the i-th entry of a vector and (·)i, j to denote the (i, j)-th entry of a matrix.
n , we repeatedly
Our first method, COfB, works as follows. After obtaining xout with data {ζt }t=1resample with replacement from the data (i.e., draw n observations from P̂n ) and run (A)SGD on the∗b , b = 1, . . . , B. As we will see in laterresampled data for B times. Denote the resample outputs by xCOfB

Lam and Wang
discussions, the number of reruns B is not necessarily large as long as it is greater than 2. Then, the 1 − γ
confidence interval for the i-th entry of x∗ is given by
h
i
COfB
COfB
γ
Ii,n
= (xout )i − tB−1,1− γ sCOfB
,
(x
)
+
t
s
(3)
out i
i
B−1,1− i
2

2

q

2 ∗
γ
1 B
1
∗b
∗b ) − (x∗γ≜ B−1where sCOfB∑Bb=1 (xCOfB
ii
COfB )i , (xCOfB )i ≜ B ∑b=1 (xCOfB )i , and tB−1,1− 2 denotes the 1− 2
quantile of the student-t distribution with degree of freedom B − 1. A pseudo-code for COfB can be foundin Algorithm 1.
Algorithm 1 Cheap Offline Bootstrap (COfB)n , number of reruns B ≥ 2, step size sequence {η }, initial guess x , nominal
Input: i.i.d. data {ζt }t=1t0coverage level 1 − γ.
COfB , i = 1, . . . , d
Output: Ii,nRun (A)SGD (2) to obtain xout .
for b ← [1, 2, . . . , B] donResample with replacement from {ζt }t=1
to obtain {ζ1∗b , . . . , ζn∗b }.
∗b
n
∗b .
Run (A)SGD for n steps on {ζt }t=1 with initialization x0 to obtain xCOfB
end forfor i ← [1, 2, . . . , d] do∗b )(x∗COfB )i ← B1 ∑Bb=1 (xCOfBiq21∗b ) − (x∗
sCOfB
)i
← B−1
∑Bb=1 (xCOfB
i
i
COfB
i
h
COfB ← (x ) − t
COfB , (x ) + t
COfB
Ii,n
out i
out i
B−1,1− 2γ si
B−1,1− 2γ si
end for
Note that COfB is an offline algorithm since resampling from {ζi }ni=1 can only be accomplished when
all these data points have been obtained. In contrast, our second method, COnB, works by maintaining
B + 1 parallel runs of ASGD starting from the same initialization. One of these trajectories is the original
run following exactly (2). The other B trajectories update similarly, except that the gradient estimate∇h(xt , ζt+1 ) is perturbed by a factor Wt,b following exponential distribution with rate 1. The confidence
COnB are constructed similarly as COfB with x∗bnintervals Ii,nout and {xCOnB }b=1 , except that the standardq2∗berror term sCOnB)i − (xout )i has xout instead of x∗COnB as the center of the squares. A≜ B1 ∑Bb=1 (xCOnB
i
pseudo-code of COnB is in Algorithm 2. Note that, when a new data ζt arrives, COnB uses only B + 1
gradient calculations to update the original and resampled outputs.
3

MAIN THEORETICAL GUARANTEESOur main theoretical guarantees on COnB and COfB is on the asymptotic coverage exactness, for B as low aseither one or two. To explain and state this result more precisely, Let Hn (·) = n1 ∑ni=1 h(x, ζi ) denote the sample
1

average approximation (SAA) of (1) and x̂n the minimizer of Hn (·). ∥x∥ p denotes (E[∥x∥ p ]) p for a randomvariable x and ∥ · ∥ denotes the standard Euclidean 2-norm for vectors. Let X1 be a bounded subset of Rd ,containing x∗ in its interior, and let X = {x| supy∈X1 ∥x − y∥ ≤ ε1 } for some ε1 > 0. For each i, j, define thefunction classes Fi, j = {∂i,2j h(x, ζ )|x ∈ X } and F˜i = {(∂i h(x1 , ζ ) − ∑ j ∂i,2j h(x2 , ζ )(x1 − x2 ))/∥x1 − x2 ∥|x1 ∈X , x2 ∈ X1 , x1 ̸= x2 }. These function classes represent the scopes of the higher-order terms of the Taylor
expansion of H at x∗ , which are crucial in developing the required asymptotic properties. Let G(x) = ∇2 H(x)and S(x) = E[∇h(x, ζ )(∇h(x, ζ ))t op] be the Hessian of H and covariance matrix of ∇h(x, ζ ) respectively.
Given n data points, define Gn (x) = n1 ∑ni=1 ∇2 h(x, ζi ) and Sn (x) = 1n ∑ni=1 ∇h(x, ζi )(∇h(x, ζi ))t op.

Lam and Wang
Algorithm 2 Cheap Online Bootstrap (COnB)n , number of bootstrap runs B ≥ 1, step size sequence {η }, initial guess x ,Input: i.i.d. data {ζt }t=1t0nominal coverage level 1 − γ.
COnB , i = 1, . . . , d
Output: Ii,nfor t ← [1, 2, . . . , n] do
xt ← xt−1 − ηt ∇h(xt−1 , ζt )for b ← [1, 2, . . . , B] do
Randomly generate Wb,t from exponential distribution with rate 1.
∗b − η W ∇h(x∗b , ζ )xt∗b ← xt−1
t b,t
t−1 t
end forend fornxout ← n1 ∑t=1
xt
for b ← [1, 2, . . . , B] do
n
∗b
xt∗b
xCOnB
← n1 ∑t=1
end forfor i ← [1, 2,q. . . , d] do
2
∗b
sCOnB
)i − (xout )i← B1 ∑Bb=1 (xCOnBihi
COnB ← (x ) − t
COnB , (x ) − t
COnB
Ii,n
out i
out i
B,1− 2γ si
B,1− 2γ si
end forAssumption 1 h and H are twice continuously differentiable in x. And the eigenvalues of ∇2 h(x, ζ ) lies
in [l, L] for some positive real numbers 0 < l < L for all x, ζ .
nAssumption 2 The noise of estimated gradient {∇h(xt−1 , ζt ) − ∇H(xt−1 )}t=1
is i.i.d. with mean 0.
Assumption 1 specifies that the objective function h exhibits strong convexity along with a bounded
Hessian, which implies the same property holds for H, in particular its strong convexity. Thus, it guarantees
the existence and uniqueness of x∗ that satisfies the first-order optimality condition ∇H(x∗ ) = 0. Assumption 2
stipulates that the evaluation noise in the first-order gradient oracle is unbiased, which is a standard assumption
to ensure the convergence of (A)SGD. To establish asymptotic normality, an additional assumption on thevariability of ∇h(x, ζ ) is required:
Assumption 3 There are τ0 , τ > 0 such that ∥x0 − x∗ ∥ ≤ τ0 and ∥∇h(x∗ , ζ )∥4 ≤ τ. The eigenvalues ofS(x∗ ) = E[∇h(x∗ , ζ )(∇h(x∗ , ζ ))t op] lie in the interval [λ1 , λ2 ] for some positive constants λ1 < λ2 .
We also need the SAA solution, namely, x̂n , to be consistent in the sense that the difference between
Hn (x̂n ) and Hn (x∗ ) converges to 0 in probability. The following assumption is sufficient for this requirement.
P

Assumption 4 supx∈Rd |Hn (x) − H(x)| −→ 0.
A further sufficient condition for Assumption 4 is that the function class {h(x, ζ )|x ∈ Rd } is GlivenkoCantelli, which can be implied by Assumption 1 if the space of x is a bounded subset of Rd (Van der Vaart
2000), though we do not assume the latter here. Essentially, a function class F is Glivenko-Cantelli if thelaw of large numbers holds uniformly in functions over F .
The following two assumptions are specialized for ASGD and SGD considered in this work respectively.
The specific choice of step size guarantees the convergence of (A)SGD in distribution. The GlivenkoCantelli assumptions help us analyze the vanishing property of some terms in our analysis of the residual
xout − x∗ .
Assumption 5 The step size satisfies ηt = ηt −α for some α ∈ ( 12 , 1]. For each i, j, function class Fi, j is
P-Glivenko-Cantelli.

Lam and Wang
Assumption 6 The step size is ηt = ηt −1 , and the initial step size η satisfies ηl > 12 . For each i, j, function
classes Fi, j and F˜i are P-Glivenko-Cantelli and P(xt = x∗ ) = 0, ∀t.
With the above assumptions, we have the following theorem:
Theorem 7 Under Assumptions 1, 2, 3, 4 and 5 for COfB running ASGD, or Assumptions 1, 2, 3, 4 and
6 for COfB running SGD, we have, for any fixed B ≥ 2, i = 1, . . . , d, the COfB 1 − γ confidence interval
for i-th entry is asymptotically exact in the sense
COfB
lim P(xi∗ ∈ Ii,n
) = 1−γn→∞(4)Moreover, under Assumptions 1, 2, 3, 4 and 5 for COnB, we have, for any fixed B ≥ 1, i = 1, . . . , d, theCOnB 1 − γ confidence interval for i-th entry is asymptotically exact in the sense
COnB
lim P(xi∗ ∈ Ii,n
) = 1−γn→∞(5)Theorem 7 states that COfB and COnB attain asymptotically exact coverage as the sample size n → ∞,regardless of any fixed choice of B ≥ 2 for COfB and B ≥ 1 for COnB. This light computation hinges
on our interval construction step at the end that differs from standard bootstraps. Note the subtlety that
COfB requires B ≥ 2, but COnB is valid even for B as small as 1. This discrepancy comes from the slightdifference in the joint asymptotic limits among the original and resample (A)SGD runs of COfB and COnB
respectively, which will be discussed in Theorem 8 in the following section. One may also notice that
COnB works only for ASGD. Whether it will work for SGD is still open to us, as the asymptotic behavior
for SGD is actually more delicate.
4

IDEA BEHIND THE MAIN GUARANTEES

We present the development for Theorem 7 in three layers. First is the cheap bootstrap idea that relies
on asymptotic independence among original and resample (A)SGD runs. Second is the conversion fromconditional convergence of resample estimates given the data, a condition that is widely utilized in classicalbootstraps, into asymptotic independence. The third and most challenging step is the development of
uniform non-asymptotic bounds to argue this conditional convergence.
4.1 Bootstrap via Asymptotic Independence
We start with the following result on asymptotic independence among original and resample runs.
Theorem 8 Under the same assumptions as Theorem 7, we have




Z0
xout − x∗
∗1



√ 
 xCOfB − x̂n  d  Z1 n .
−−−→  .  ..
 n→∞  .. 
∗B
xCOfB
− x̂n

and

(6)

ZB




Z0
xout − x∗
∗1



√ 
 xCOnB − xout  d  Z1 n .
−−−→  .  ..
 n→∞  .. 


∗B
xCOnB
− xout(7)

ZB

where Zb , b = 0, . . . , B are i.i.d. d-dimensional Gaussian random variables with mean 0. When xout standsnfor the ASGD output n1 ∑t=1
xt , the covariance matrix of Zb is G(x∗ )−1 S(x∗ )G(x∗ )−1 , b = 0, . . . , B, whereG(x∗ ) = ∇2 H(x∗ ) and S(x∗ ) = E[∇h(x∗ , ζ )(h(x∗ , ζ ))⊤ ]. When xout stands for the SGD output xn , consider

Lam and Wang
the singular value decomposition G(x∗ ) = QDQ⊤ with D = diag(d1 , . . . , dd ), where d1 , . . . , dd are eigenvaluesof G(x∗ ) in decreasing order and Q the matrix consisting of eigenvectors. For i, j = 1, . . . , d and b = 0, . . . , B,
the covariance between (Zb )i and (Zb ) j is given by σi,2 j = η 2 (ηdi + ηd j − 1)−1 (Q⊤ S(x∗ )Q)i, j .
Theorem 8 involves two aspects. First, the asymptotic distribution of the error of the resample run∗bcompared with the SAA solution in COfB, namely, xCOfB− x̂n (or compared with the original ASGD run∗b
in COnB, xCOnB − xout ) is the same as that of the original run compared with the true minimizer, xout − x∗ .
Second, more importantly, is the asymptotic independence among all these errors. Derivation of Theorem
8 will be discussed in the later sections. Taking Theorem 8 for granted, the following is a proof of Theorem
7.
Proof of Theorem 7.
Consider COfB first. Observe that√√√n ((xout )i − xi∗ )n ((xout )i − xi∗ )n ((xout )i − xi∗ )(xout )i − xi∗
q
r
r
=
=
=
√
2
2
sCOfB
∗b ) −(x̂ ) −√n (x∗i∑B ((x∗b )i −(x∗ )i )∑Bb=1 ( n((xCOfB( COfB )i −(x̂n )i ))n × (sCOfB
)2
i
n i)
i
n × b=1 COfBB−1 COfB
B−1
As n goes to infinity, we have√n ((xout )i − xi∗ )N(Z0 )idddr= q 2 = tB−1
→q B
2
√
√
2
χ
((Z
)
−(
Z̄)
)
B
∗b ) −(x̂ ) − n (x∗∑b=1 b ii
B−1∑b=1 ( n((xCOfB( COfB )i −(x̂n )i ))
i
n i)
B−1
B−1

B−1

2
a χ 2 -variable with B − 1
where Z̄ = (1/B) ∑Bb=1 Zb , and N stands for a standard normal variable, χB−1ddegree of freedom, tB−1 a student t-variable with B − 1 degree of freedom, and “=” equality in distribution.
The convergence in distribution above comes from the continuous mapping theorem. The first equality in
distribution comes from the i.i.d. normality limit in Theorem 8 and the elementary relation between χ 2
and normal. The second equality in distribution comes from the elementary construction of a t variable.
Thus, by a pivotal argument, we obtain the confidence interval generated from COfB.
A similar argument works for COnB, except that we use xout directly in place of x∗COfB in the pivotal
construction and correspondingly, it would result in student t-distribution with degree of freedom B.
4.2 From Conditional Convergence to Asymptotic Independence
Our next step is to prove Theorem 8. As a sub-step, the conclusions in Theorem 8 can be implied by aconditional convergence:
Theorem 9 Suppose√dn(xout − x∗ ) −−−→ Z0(8)n→∞

In addition, if√ ∗bdn(xCOfB − x̂n∗ ) −−−→ Z0 conditional on ζ1 , ζ2 , . . .

(9)

√ ∗b
d
n(xCOnB − xout ) −−−→ Z0 conditional on ζ1 , ζ2 , . . .

(10)

n→∞

then (6) holds, and ifn→∞then (7) holds, where Z0 denotes a d-dimensional Gaussian random variable with mean 0 and covariancematrix as described in Theorem 8.
Therefore, if we can show (8), (9) and (10), then we obtain the conclusions in Theorem 8 andsubsequently the guarantees of COfB and COnB in Theorem 7. Note that (8) is the classical asymptotic
normality of (A)SGD guaranteed by our assumptions (Chung 1954; Sacks 1958). On the other hand, theLam and Wang
type of conditional convergence in (9) and (10) is the main driver of classical bootstrap methods that allow
the approximation of a sampling distribution using the resample counterpart. Here, converting the latter
into the conclusion in Theorem 8 presents a new bootstrap methodological route that substantially reduces
B. The proof for Theorem 9 generalizes the proof of Proposition 1 in Lam (2022b) and we refer readers
to the details therein.
4.3 From Uniform Non-Asymptotic Bound to Conditional Convergence
For COnB, the desired conditional convergence result (10) is well established in the proof of Theorem 1
in Fang et al. (2018). We focus on proving (9) for COfB under both SGD and ASGD settings, whichconstitutes our main technical development.
Let us abstractize our discussion to denote ψ(P) as the minimizer for (1) with data distribution P,
where ψ is viewed as a mapping from the data distribution to Rd . Correspondingly, define ψn as the
mapping from the data distribution to the outcome of (A)SGD. Then ψn (P) ∈ Rd is the (random) outcomenof (A)SGD after n iterations, as a function of data distribution P with h and {ηt }t=1
implicitly
chosen.
√Classical results (Polyak and Juditsky 1992; Chung 1954) state that the weak limit of n(ψn (P)− ψ̂(P))
exists and equals Z0 . This Z0 is the Gaussian variable described in Theorems 8 and 9 whose variancedepends on P. Let Ẑm denote the normal variable that replaces P in its variance with P̂m , conditional on
the collected data. With these new notations, (9) holds if for any Borel measurable set D ⊂ Rd , we havethe following√lim |P∗ ( n ψn (P̂n ) − ψ̂(P̂n ) ∈ D) − P(Z0 ∈ D)| = 0 w.p.1
(11)
n→∞

where P∗ denotes the probability conditional on the data (we will also use E∗ to denote the corresponding
conditional expectation). By the triangle inequality, one can obtain
√|P∗ ( n ψn (P̂n ) − ψ̂(P̂n ) ∈ D) − P(Z0 ∈ D)|√≤|P∗ ( n ψn (P̂n ) − ψ̂(P̂n ) ∈ D) − P∗ (Ẑn ∈ D)| + |P∗ (Ẑn ∈ D) − P(Z0 ∈ D)|
It can be proved that the second term above vanishes with probability 1. On the other hand, we have thefollowing theorem for the first term:
Theorem 10 Under the same assumptions as in theorem 7 and focusing on COfB, for any measurable set
D, we have√lim |P∗ n(ψn (P̂n ) − ψ̂(P̂n )) ∈ D − P∗ (Ẑn ∈ D)| = 0 w.p.1(12)n→∞
The proof invokes an expansive analysis on the behavior of the (A)SGD output. From the iterative
scheme (2), one obtains the followingnnxn+1 = B0n x1 − ∑ ηm Bmn δm − ∑ ηm Bmn Em
m=1

(13)

m=1

where δk ≜ δ (xk ) = ∇H(xk ) − G(x∗ )(xk − x∗ ) is the second-order residual of the Taylor expansion of ∇H
at x∗ , Ek−1 = ∇h(xk−1 , ζk ) − ∇H(xk−1 ), and Bmn = ∏nj=m+1 (I − η j G) ∈ Rd×d .
In the ASGD case, from (13) and using results from Shao and Zhang (2022), we show that there is a
4-tuple (τˆ0 , τ̂, Ĉ, N) such that√sup |P∗ ( n ψn (P̂n ) − ψ̂(P̂n ) ∈ D) − P∗ (Ẑn ∈ D)| ≤ Ĉ(d 3/2 + τ̂ 3 + τ̂03 )(d 1/2 n−1/2 + n−α+1/2 )
(14)
n>N

for any measurable D. Shao and Zhang (2022) gave an inequality similar to (14) but with a fixed distribution
instead of a varying distribution P̂n depending on n. Additional arguments are required to show that theLam and Wang
Table 1: Average time for different methods for logistic regression.
MethodAverage Time (s)

delta
10.53

BM
1.10

OB (B = 100)
104.49HiGrid(2,2)
0.95

COfB ASGD
3.15

COnB
4.09

established rate is uniform across all data distributions including the empirical distribution. Detailed prooffor (14) is omitted here due to space limit.
For the SGD case, the first two terms in (13) correspond to the interaction of the error of the initialsolution and second-order residual in the Taylor expansion of ∇H. One can show the following vanishingpropertyn

lim E[∥B0n x1 − ∑ ηm Bmn δm ∥] = 0

n→∞

m=1
n
On the other hand, the last term ∑m=1 ηm Bmn Em consists of the difference between sample gradient ∇hand true gradient ∇H, which converges to a normal distribution. We use the Berry-Esseen-type result
from Lemma 4 in Shao and Zhang (2022) to give a non-asymptotic convergence result for this term andthus establish a bound similar to (14) that is uniform across all data distributions including the empirical
distribution.

(a)

(b)

Figure 1: Performance of methods with respect to different choices of initial step size η, with sample sizen = 104 . We compare the delta method, batch-mean method with M = n0.25 , COfB method with B = 10,
COnB with B = 3, 10, and the online bootstrap method with B = 10, when n = 104 . The left figure shows
the sensitivity of the average coverage probability against η; The right figure shows the sensitivity of theaverage length of confidence interval against η. We report the result for identity Σ.

5EXPERIMENTS

In this section, we illustrate the numerical performances of our approaches and compare with the other⊤methods. We consider logistic regression, with loss function h(x, ζ ) = log (1 + e−b×a x ). The data ζ = (a, b)
coming from distribution P consists of the independent variable a ∈ Rd and dependent variable b ∈ {−1, 1}.
where a follows a multivariate normal distribution with mean 0 and covariance matrix Σ, and b = 1 with1aa⊤probability. In this case, we have ∇h(x, ζ ) = −b×aand ∇2 h(x, ζ ) =. The−a⊤ x∗
b×a⊤ x
a⊤ x
−a⊤ x
1+e

1+e
2
Hessian information ∇ h(x, ζ ) above will only be used in the delta method.

(1+e

)(1+e

)

Baselines.
The batch-mean method (BM) splits {xi }ni=0 into M batches, with ek and sk denoting
the ending index and starting index of the k-th batch respectively. nk denotes the number of iterates in
1 ek⊤
k-th batch, and the estimator is defined to be M1 ∑M
k=1 nk (x̄nk − x̄M )(x̄nk − x̄M ) , where x̄nk = nk ∑i=sk xi andLam and Wang
1−α

1

eM
1
n
x̄M = eM −e0, and ek to be the closest integer to ((k + 1)N) 1−α for each k = 0, . . . , M as
x . Let N = M+1∑i=s
1 isuggested in Chen et al. (2020). The confidence interval for each entry of x∗ is constructed using diagonal
entries of the batch-mean estimator and a normal quantile. The delta method (Chen et al. 2020) generate
1 n−12
confidence intervals using normal quantile and Σ̃2n = G̃−1n S̃n G̃n , where G̃n = n ∑i=1 ∇ h(xi−1 , ζi ), and
1 nS̃n = n ∑i=1 ∇h(xi−1 , ζi )(∇h(xi−1 , ζi ))⊤ are computed on the fly. The online bootstrap method (OB) (Fang(b)

et al. 2018) works by generating B sequences of i.i.d. exponential random variables, {Wi }ni=1 , b = 1, . . . , B.
(b)
And for each b, run ASGD using {ζi }ni=1 and step sizes {Wi ηi }ni=1 . The confidence interval is thenconstructed using the outcomes and normal quantile. The HiGrad method (Su and Zhu 2018) takes twotuples (B1 , B2 , . . . , BK ) and (n0 , n1 , . . . , nK ) as hyperparameters describing when to break the SGD thread
into how many branches. Bk represents the number of branches a single branch divides into at the k-th
break, and nk represents the number of data each thread uses between the k-th and (k + 1)-th breaks. Afterall the breaks, there will be T = ∏Kk=1 Bk threads and one obtain {x( j) }tj=1 op by averaging each thread.
The confidence interval for the i-th entry of x∗ is calculated by aggregating {(x( j) )i }tj=1 op.
Hyperparameters. The nominal coverage probability is set to be 95%. We report results for threechoices of dimensions, d = 5, 20, 200. Two choices of Σ are tested, namely identity (Σ = Id ) and Toeplitz(Σi, j = 0.5|i− j| ). For all experiments, the decay rate of step size α is set to be 0.501, so ηt = ηt −0.501 . We
12
set x∗ = [0, d−1, d−1, . . . , 1] and x0 = [0, . . . , 0] ∈ Rd .
For each set of hyperparameters, we run 500 independent trials and report the mean and standard
deviation of the coverage probabilities and the average length of the intervals across d dimensions. We tune
the initial step size η within the range [0.2, 0.7] and report the result with the most accurate average coverage
probability. For the batch-mean method, M is selected to be the nearest integer to n0.25 as suggested in
Chen et al. (2020). For HiGrad, the architecture we experiment on is ((2, 2), (n/7, n/7, n/7)). We report
the performance of COfB and COnB with B = 3, and set B = 200 for the online bootstrap method unlessotherwise specified, as suggested in Fang et al. (2018).
Results.
Results for this logistic regression experiment can be found in Table 2. We use bold
numbers to denote good results in the sense that the coverage probability is between 92% and 98%, anduse italic numbers to denote bad results with coverage probability less than 80%.
Generally speaking, for any d and Σ discussed in this experiment, the delta method and batch-mean
method are outperformed by other methods. For the HiGrad method, we observe a significant drop in
performance when d = 200, which might be caused by its shorter SGD trajectory compared with othermethods. When d = 200, the number of data/iterations is inadequate for SGD to converge properly. The
performance of the delta method drops significantly as d increases since it requires a Monte-Carlo estimationof a full Hessian matrix, the accuracy of which suffers when the dimension increases, while the benefitof reusing data in bootstrap-type methods is more apparent. For the online bootstrap method, although it
gives a comparable coverage probability, our methods are significantly faster. As a trade-off, the average
length of the confidence interval of our method is larger compared with the batch-mean method and thedelta method, which is due to the t-quantile with degree of freedom B − 1 (B) in our COfB (COnB). Recall
that entries of x∗ are linearly spaced in [0, 1], and the confidence intervals are of magnitude 10−2 . Thus
the benefits in computational efficiency appear to significantly outweigh the increase in interval length.
Table 1 presents the computation time for different methods. The experiment was performed on a single
processor at 3.4GHz. HiGrad and the batch-mean method are the fastest as they require no extra gradientsteps or computing Hessian. Roughly speaking, the computation time is proportional to the number ofgradient steps except for the delta method.
Sensitivity Analysis. In Figure 1, we compare the performance of our methods, the delta method,
the batch-mean method, and the online bootstrap method with number of samples n = 104 for the linear
regression problem. Observe that the coverage probability of COfB methods remains stable around 95%
regardless of changes in the initial step size. On the other hand, the batch-mean method requires a carefulchoice of the initial step size to give a comparable coverage rate, and the optimal choice is not the same

Lam and Wang
Table 2: Coverage probabilities and average interval lengths in the logistic regression experiment.
Cov (%)

d=5
Len (×10−2 )

delta
BM
OB
HiGrad(2,2)
COfBASGD
COfBSGD
COnB

95.00 (0.07)
89.33 (0.10)
94.83 (0.07)
94.33 (0.07)
94.12 (0.07)
95.40 (0.07)
94.33 (0.07)

delta
BM
OB
HiGrad(2,2)
COfBASGD
COfBSGD
COnB

94.83 (0.07)
84.00 (0.12)
95.00 (0.07)
95.33 (0.07)
94.12 (0.07)
95.36 (0.07)
94.00 (0.08)

d = 20Cov (%)Len (×10−2 )Identity Σ, n = 105
3.10 (0.00) 94.12 (0.07) 3.68 (0.00)
2.64 (0.00) 87.29 (0.11) 3.11 (0.01)
3.18 (0.00) 96.96 (0.05) 4.41 (0.01)
5.77 (0.02) 95.46 (0.07) 7.10 (0.02)
4.21 (0.00) 95.03 (0.07) 7.32 (0.01)
9.38 (0.01) 94.99 (0.07) 9.42 (0.01)
4.71 (0.02) 95.62 (0.06) 6.56 (0.03)
Toeplitz Σ, n = 105
4.05 (0.00) 93.29 (0.08) 5.59 (0.00)
3.16 (0.01) 75.25 (0.14) 3.75 (0.01)
4.24 (0.00) 94.67 (0.07) 6.70 (0.01)
7.18 (0.03) 93.38 (0.08) 8.92 (0.03)
5.70 (0.00) 94.77 (0.07) 11.49 (0.01)
9.48 (0.01) 94.80 (0.07) 9.65 (0.01)
6.22 (0.02) 94.71 (0.07) 10.27 (0.05)

d = 200
Cov (%)Len (×10−2 )61.92 (0.15)
57.47 (0.16)
99.91 (0.01)
80.61 (0.13)
92.34 (0.08)
94.72 (0.07)
99.42 (0.02)

5.85 (0.00)
5.56 (0.02)
50.71 (0.23)
10.26 (0.04)
19.00 (0.02)
25.61 (0.03)
75.25 (0.43)

53.69 (0.16)
34.93 (0.15)
99.78 (0.01)
57.02 (0.16)
93.79 (0.08)
94.41 (0.07)
97.82 (0.05)

9.56 (0.00)
7.30 (0.03)
69.28 (0.26)
10.27 (0.04)
42.27 (0.05)
30.53 (0.03)
99.61 (0.60)

across different problems. The delta method suffers from a huge under-coverage and fails to give a valid
confidence interval. The delta method and the batch-mean method have smaller average lengths, which
can be associated with their under-coverages. It can be observed that the coverage probability and theaverage length of the batch-mean estimator both increase as η increases. Our COnB method has a similarsensitivity as the online bootstrap method. Nonetheless, as mentioned earlier, COnB is substantially faster.
Additionally, the average length of our COnB becomes almost the same as that of the online bootstrap
method when increasing B to 10.
ACKNOWLEDGEMENTS
We gratefully acknowledge support from the InnoHK initiative, the Government of the HKSAR, andLaboratory for AI-Powered Financial Technologies.
REFERENCES
Anastasiou, A., K. Balasubramanian, and M. A. Erdogdu. 2019. “Normal Approximation for Stochastic Gradient Descent via
Non-asymptotic Rates of Martingale CLT”. In Conference on Learning Theory, 115–137. PMLR.
Chen, J. X., and M. Lopes. 2020, 13–18 Jul. “Estimating the Error of Randomized Newton Methods: A Bootstrap Approach”.
In Proceedings of the 37th International Conference on Machine Learning, edited by H. D. III and A. Singh, Volume 119
of Proceedings of Machine Learning Research, 1649–1659: PMLR.
Chen, X., J. D. Lee, X. T. Tong, and Y. Zhang. 2020. “Statistical Inference for Model Parameters in Stochastic Gradient
Descent”. The Annals of Statistics 48(1):251–273.
Chung, K. L. 1954. “On a Stochastic Approximation Method”. The Annals of Mathematical Statistics:463–483.
Fang, Y., J. Xu, and L. Yang. 2018. “Online Bootstrap Confidence Intervals for the Stochastic Gradient Descent Estimator”.
Journal of Machine Learning Research 19(78):1–21.
Flegal, J. M., and G. L. Jones. 2010. “Batch Means and Spectral Variance Estimators in Markov Chain Monte Carlo”. The
Annals of Statistics 38(2):1034–1070.
Geyer, C. J. 1992. “Practical Markov Chain Monte Carlo”. Statistical Science 7(4):473–483.
Glynn, P. W., and D. L. Iglehart. 1990. “Simulation Output Analysis using Standardized Time Series”. Mathematics of Operations
Research 15(1):1–16.

Lam and Wang
Glynn, P. W., and H. Lam. 2018. “Constructing Simulation Output Intervals under Input Uncertainty via Data Sectioning”. In
Proceedings of the 2018 Winter Simulation Conference, 1551–1562: Institute of Electrical and Electronics Engineers, Inc.
Jones, G. L., M. Haran, B. S. Caffo, and R. Neath. 2006. “Fixed-width Output Analysis for Markov Chain Monte Carlo”.
Journal of the American Statistical Association 101(476):1537–1547.
Lam, H. 2022a. “Cheap Bootstrap for Input Uncertainty Quantification”. In Proceedings of the 2022 Winter Simulation
Conference, 2318–2329: Institute of Electrical and Electronics Engineers, Inc.
Lam, H. 2022b. “A Cheap Bootstrap Method for Fast Inference”. arXiv preprint arXiv:2202.00090.
Lam, H., and Z. Liu. 2023. “Bootstrap in High Dimension with Low Computation”. International Conference on Machine
Learning (ICML).
Lattimore, T., and C. Szepesvári. 2020. Bandit Algorithms. Cambridge University Press.
Li, T., L. Liu, A. Kyrillidis, and C. Caramanis. 2018, Apr.. “Statistical Inference Using SGD”. Proceedings of the AAAI
Conference on Artificial Intelligence 32(1).
Lunde, R., P. Sarkar, and R. Ward. 2021. “Bootstrapping the Error of Oja’s Algorithm”. Advances in Neural Information
Processing Systems 34:6240–6252.
Nesterov, Y., and J. P. Vial. 2008, jun. “Confidence Level Solutions for Stochastic Programming”. Automatica 44(6):1559–1568.
Oja, E. 1982. “Simplified Neuron Model as a Principal Component Analyzer”. Journal of Mathematical Biology 15(3):267–273.
Polyak, B. T., and A. B. Juditsky. 1992. “Acceleration of Stochastic Approximation by Averaging”. SIAM Journal on Control
and Optimization 30(4):838–855.
Rakhlin, A., O. Shamir, and K. Sridharan. 2012. “Making Gradient Descent Optimal for Strongly Convex Stochastic Optimization”.
In Proceedings of the 29th International Conference on International Conference on Machine Learning, 1571–1578.
Robbins, H., and S. Monro. 1951. “A Stochastic Approximation Method”. The Annals of Mathematical Statistics 22(3):400–407.
Rumelhart, D. E., G. E. Hinton, and R. J. Williams. 1986. “Learning Representations by Back-propagating Errors”. Nature 323(6088):533–536.
Sacks, J. 1958. “Asymptotic Distribution of Stochastic Approximation Procedures”. The Annals of Mathematical Statistics 29(2):373–405.
Schmeiser, B. 1982. “Batch Size Effects in the Analysis of Simulation Output”. Operations Research 30(3):556–568.
Schruben, L. 1983, dec. “Confidence Interval Estimation Using Standardized Time Series”. Operations Research 31(6):1090–1108.
Shao, Q.-M., and Z.-S. Zhang. 2022. “Berry–Esseen Bounds for Multivariate Nonlinear Statistics with Applications to M-estimators
and Stochastic Gradient Descent Algorithms”. Bernoulli 28(3):1548–1576.
Su, W. J., and Y. Zhu. 2018. “Uncertainty Quantification for Online Learning and Stochastic Approximation via Hierarchical
Incremental Gradient Descent”. arXiv preprint arXiv:1802.04876.
Van der Vaart, A. W. 2000. Asymptotic Statistics, Volume 3. Cambridge university press.
Zhu, Y., and J. Dong. 2021. “On Constructing Confidence Region for Model Parameters in Stochastic Gradient Descent Via
Batch Means”. In Proceedings of the 2021 Winter Simulation Conference, 1–12: Institute of Electrical and Electronics
Engineers, Inc.

AUTHOR BIOGRAPHIESHENRY LAM is an Associate Professor in the Department of Industrial Engineering and Operations Research at Columbia
University. His research interests include Monte Carlo methods, uncertainty quantification, risk analysis and data-driven
optimization. He serves as the area editor for Stochastic Models and Data Science in Operations Research Letters, and on theeditorial boards of Operations Research, INFORMS Journal on Computing, Applied Probability Journals, Stochastic Models,
Manufacturing and Service Operations Management, and Queueing Systems. His email address is henry.lam@columbia.edu
and his website is http://www.columbia.edu/~khl2114/.
ZITONG WANG is a Ph.D. student of Industrial Engineering and Operations Research at Columbia University. His primaryresearch interests are stochastic optimization and uncertainty quantification. His email address is zw2690@columbia.edu.