\section{Introduction}
In this paper, we consider the following strongly convex optimization problem:
\begin{align}\label{Problem1} \min_{x\in \mathbb{R}^m} f(x) \triangleq   \mathbb{E} [\us{F}(x,\xi)],
\end{align}
where  $\xi:  \Omega \to \Real^d$  is a random variable \blue{defined on the  probability space} $({ \Omega}, {\cal F},\mathbb{P})$, $F: \mathbb{R}^m \times \mathbb{R}^d \to \mathbb{R}$, and   the expectation is taken over the distribution of the random vector. Stochastic optimization problems  have been  extensively studied
 given  wide applicability in almost all areas of
science and engineering, ranging from communication and queueing systems  to
finance (cf.~\cite{birge97introduction,shapiro2009lectures}). However, in most
situations,  this expectation  \blue{along with} its derivative is unavailable in closed
form requiring the development of sampling approaches. Sample  average approximation (SAA)~\cite{shapiro2009lectures} and   stochastic approximation (SA)  represent two  commonly used approaches  for contending with
stochastic programs.
Our focus is on SA  schemes,  first considered by Robbins and Monro
\cite{robbins1951stochastic} for seeking roots of a regression
function with noisy observations.  The   standard  SA algorithm
$x_{k+1}=x_k-\alpha_k \nabla F(x_k,\xi_{k}),$ also known as  the
stochastic gradient  descent (SGD) algorithm,   updates the estimate
$x_{k+1}$ based on  a sampled gradient  $\nabla F(x_k,\xi_{k})$.  The  convergence analysis usually requires
suitable properties on the gradient map (such as Lipschitzian
requirements) and the steplength sequence (such as non-summable but
squuare summable). The almost sure convergence of $x_k$ to
$x^*$, the unique optimal solution of \eqref{Problem1},  was
established in
\cite{blum1954multidimensional,bottou1998online,chen2006stochastic}
on the basis of the  Robbins-Siegmund theorem
\cite{robbins1971convergence}   while ODE techniques were
employed for claiming similar statements   in
\cite{hasminskii1973stochastic,kushner2003stochastic,borkar2009stochastic}.
Tight bounds on the  rate of convergence can be obtained by establishing
the asymptotic distribution for the iterates (cf.
\cite{chen2006stochastic,hasminskii1973stochastic,kushner2003stochastic,borkar2009stochastic,fabian1968asymptotic}).
An instructive review of results up to around 2010 is provided by
\cite{pasupathy2011stochastic},  while  \cite{fabian1968asymptotic,hasminskii1973stochastic,kushner2003stochastic}
proved  the asymptotic normality of the suitably  scaled  iterates for   SA
with decreasing step-sizes. To be specific,  the rescaled   error process
$\sqrt{ \alpha_k}(x_k-x^*)$ asymptotically converges  in distribution to a
normal distribution  with zero mean  \blue{and} with covariance depending on the
Hessian matrix, the covariance of the gradient noise, and the steplength.
The asymptotic normality was further investigated in \cite{chen2006stochastic}
for SA with expanding truncations, an \blue{avenue} that does not require
Lipschitz continuity   of the gradients.    The sequence of iterates
generated by  the  constant steplength SA scheme  is  shown to be a
homogeneous Markov chain with a unique stationary distribution;  see
\cite[Chapter 9]{kushner2003stochastic},  \cite[Chapter 17]{meyn2012markov},
and \cite{dieuleveut2017bridging}.\\

Central limit results for SA schemes are significant from
the standpoint of algorithm design as well as inference. Indeed, the optimal
selection of the steplength depends on the Hessian matrix at  $x^*$, see e.g.,
\cite[Chapter 3.4]{chen2006stochastic}. Since the  Hessian at $x^*$ is
unavailable, the optimal value of $\alpha_k$ can only be estimated, which has  led to the development of   adaptive SA methods (cf. \cite{venter1967extension,wei1987multivariate,spall2000adaptive}).
Motivated by heavy dependence of SA schemes on steplength choices,
\cite{yousefian10stochastic} developed a self-tuned rule that
adapts the steplength sequence to problem parameters. In a similar vein,
there has been an effort to develop optimal constant steplengths.
Specifically,    it was shown \blue{in \cite{nemirovski2009robust}}  that with   suitably selected
 constant  step-sizes, the  expected  function values   at the averaged
iterate  converge to the optimum with rate  $\mathcal{O}(1/\sqrt{k})$ in
merely convex case and  rate $\mathcal{O}(1/k)$  in  strongly convex case,
matching the lower bounds~\cite{nemirovskii1983problem}.
\blue{It is further  shown in   \cite{mandt2017stochastic} that the constant SGD  simulates a Markov chain
with a stationary distribution, which   might   be  used  to adjust the tuning parameters of constant SGD
 so as to improve the  convergence rate,}    e.g.,  SGD for fitting  the  generalized linear
models \cite{toulis2014statistical}   and approximating the  Bayesian posterior
inference \cite{mandt2017stochastic}.  Furthermore, the central limit theorems (CLTs) for SA
schemes might  allow for  the possibility  of constructing  confidence
regions for the optimal solution,   see
\cite{bickel2013asymptotic,hsieh2002recent}.   For instance, Hsieh and
Glynn~\cite{hsieh2002recent}   designed an approach to rigorously
characterize confidence regions without explicitly estimating the covariance
of the limiting normal distribution.\\

Unfortunately,  SA schemes with  diminishing steps cannot  recover  the  deterministic convergence rates
 seen in exact gradient methods  while constant
steplength SA schemes are only characterized by convergence
guarantees to a  neighborhood of the optimal solution.
 Variance-reduction schemes employing an increasing batch-size of
sampled gradients (instead of the unavailable true gradient) appear
to have been first alluded to
in~\cite{bertsekas1989parallel,hdmello03variable,deng09variable,pasupathy10choosing} and
analyzed in smooth and strongly convex~\cite{friedlander12hybrid,
byrd12sample, shanbhag2015budget, schmidt12convergence}, smooth
convex regimes~\cite{ghadimi2016accelerated}, nonsmooth (but
smoothable) convex~\cite{jalilzadeh2018optimal},
nonconvex~\cite{lei18asynchronous},  and
game-theoretic~\cite{lei18game} regimes. Notably,  the linear rates in mean-squared error were  derived
in for strongly convex smooth~\cite{shanbhag2015budget} and a
subclass of nonsmooth objectives~\cite{jalilzadeh2018optimal}, while
a   rate $\mathcal{O}(1/k^2)$ and $\mathcal{O}(1/k)$ was obtained
for expected sub-optimality  in convex
smooth~\cite{ghadimi2016accelerated, jofre19variance} and
nonsmooth~\cite{jalilzadeh2018optimal}, respectively.  In each instance, the
schemes achieve the corresponding optimal deterministic rates of
convergence under suitable growth in
sample-size sequences while in almost all cases, the optimal sample-complexity bounds were obtained.  An excellent
discussion relating sampling rates and the canonical rate (i.e. the Monte-Carlo rate) in stochastic
gradient-based schemes has been provided by Pasupathy et
al.~\cite{pasupathy18sampling}.\\

 {\bf Gaps and motivation.}  This paper is motivated by
the following gaps. (i) Despite a surge of interest in
variance-reduced schemes (with and without acceleration), no limit
theorems are available for claiming asymptotic normality of the
limiting scaled sequence in either \blue{unaccelerated or} Nesterov accelerated regimes coupled
with variance reduction. In particular, we have little understanding
regarding whether such avenues have detrimental impacts in terms of
the limiting behavior. Such statements are also unavailable for the
related heavy-ball method. (ii) The geometric growth in sample-size
required for achieving linear rates of convergence may prove onerous
in some settings.  Are CLTs available in settings where the
sample-size grows at a polynomial rate? (iii) Finally, given a CLT,
there is little by way of availability of rigorous confidence
statements, barring the work by Hsieh and
Glynn~\cite{hsieh2002recent}. Can such statements be developed for
the proposed class of variance-reduced first-order methods?\\

  {\bf Outline and contributions.}  To address these
gapes, we present CLTs and confidence statements for   first-order
stochastic variance-reduced algorithms for resolving
\eqref{Problem1}, including the classical
SGD~\cite{robbins1951stochastic}, the stochastic variants of  the
Nesterov's  accelerated  method~\cite{nesterov1983method}  and the
heavy ball method~\cite{polyak1964some}.  We provide statements
when batch-sizes increases at either a geometric or a polynomial
rate. Our main contributions are summarized next.

\noindent   {{\bf (I) CLTs for variable sample-size gradient
methods.}  In Section \ref{Sec:Alg:SGD}, we recall the
variance-reduced (VR) stochastic gradient algorithm
\blue{with  a constant steplength}  (see Algorithm
\ref{Alg_1}), where the  gradient is  estimated by  the average
of  an  increasing  batch of sampled gradients. In Section
\ref{sec:CLT-geo}, when the  batch-size increases at a
geometric rate, the mean-squared error
diminishes at a geometric   rate (see Proposition
\ref{prp1}) and   provide  a preliminary Lemma
\ref{lem-CLT2}  for  establishing  a  CLT for a
noised-corrupted     linear  recursion. Based on  this Lemma
and    the linear approximation  of the gradient function
at  $x^*$,  we  may derive a  CLT (see Theorem
\ref{thm-CLT1}) in this setting. We proceed to show}  that
the  covariance   of the limiting normal distribution
depends on the  Hessian matrix at the solution,  the
condition number,  the covariance of   gradient noise, etc
and the steplength.  Additionally, we show in    Section
\ref{sec:CLT-poly} that when the batch-size  is increased at
a polynomial rate, the sequence of iterates  converge
at  a corresponding   polynomial  rate (see Proposition
\ref{prp3}).  Then based on the CLT shown in Lemma
\ref{lem-CLT3}  for the time-varying linear recursion, the
CLT for Algorithm \ref{Alg_1}  with polynomially increasing
batch-size is established in  Theorem \ref{thm-CLT3}.

\noindent {\bf (II)  CLTs for VR-accelerated gradient method.} In Section \ref{Sec:Alg:ACC},  we consider a VR-accelerated
gradient algorithm with constant steplengths (see Algorithm \ref{Alg_2}) and  by
leveraging the  geometric rate of convergence (Prop.~\ref{prp2}),
we establish amongst the first CLTs in accelerated regimes
(Theorem \ref{thm-CLT1-Alg2})  when the batch-size increases at a
geometric rate.   It is well-known that the accelerated variant
has better constants in the iteration complexity in deterministic regimes and this is also seen in stochastic settings.
 Akin to earlier,  when the  batch-size increases at a polynomial rate, a polynomial convergence   rate and the corresponding CLT are   in  Proposition \ref{prp4}   and Theorem \ref{thm-CLT2-Alg2}.

\noindent {\bf III. CLTs for VR-heavy-ball schemes.} In
Section \ref{Sec:Alg:HB},   we design  a VR-heavy ball scheme with constant steplengths (see Algorithm
\ref{Alg_3}).  In contrast with  Algorithms \ref{Alg_1} and \ref{Alg_2},  the heavy ball method imposes  a
twice-differentiability assumption on $f$.  When the batch-size
increases geometrically (or polynomially), a  geometric (or
polynomial)  rate  of convergence is shown in Proposition \ref{prp-alg3}
(or Proposition \ref{prp2-alg3}),   while  a CLT is
established in Theorem \ref{thm-CLT1-Alg3}  (or Theorem \ref{thm-CLT2-Alg3}).

\noindent {\bf IV. Confidence statements.} In Section
\ref{sec:CI}, inspired by  \cite{hsieh2002recent},  we
 provide rigorous  confidence regions  for the optimal solution
and  optimal function value.  Then in Section \ref {sec:simu},
we  implement  some simulations on a  parameter estimation  problem
in the stochastic  environment to validate the theoretical findings.\\

  \blue{{\bf Notations.} Let $\mathbf{I}_m$ denote the identity matrix of dimension $m $ and $\mathbf{0}_m \in \mathbb{R}^{m \times m}$  denote   the   matrix with all entries equal zero.
  Let $\{X_k\}$ be a sequence of   random variables.   $X_k  \xrightarrow [k \rightarrow \infty]{d}   N(0,\mathbf{S}) $ denotes that $ X_k $ converges  in distribution  to a  normal distribution  $N(0, \mathbf{S})$   with mean $0$ and covariance $\mathbf{S}, $ and $X_k  \xrightarrow [k \rightarrow \infty]{P} X $ denotes that  $ X_k $ converges  in probability  to $X.$}

\section{First-Order Variable Sample-size  Stochastic   Algorithms}\label{Sec:Alg}
Since the exact gradient $\nabla f(x)$  is  expectation-valued and unavailable
in a closed form,   we assume  there exists a {\em stochastic first-order
oracle}   such that  for any $x,\xi$, a  sampled gradient $\nabla F (x,\xi)$ is returned,
assumed to be an unbiased estimator of  $\nabla f(x)$.  In
this section, we present three  first-order  stochastic  algorithms to
find the optimal solution to \eqref{Problem1}. Throughout the
paper, time is slotted at $k=0,1,2,\dots$ and an iterate at time $k$ is denoted by $x_k\in \mathbb{R}^m.$

\subsection{Gradient Descent Method}\label{Sec:Alg:SGD}
 We  recall  a variable sample-size stochastic gradient algorithm  (Algorithm \ref{Alg_1}) to solve \eqref{Problem1}, where at iteration $k $,   the unavailable exact gradient $\nabla f(x_k)$ is estimated via the average of  an increasing batch-size of  sampled gradients.
 \begin{algorithm}[H]  \caption {Variance reduced SGD}\label{Alg_1}
Given an arbitrary  initial value $x_0\in \mathbb{R}^m, $  a positive constant $\alpha>0,$ and a positive  integer sequence $\{N_k\}_{k\geq 0}.$ Then iterate the following equation  for $k\geq 0.$
\begin{align}\label{Alg1} x_{k+1}=x_k - \alpha \tfrac{ \sum_{j=1}^{N_k} \nabla F(x_k,\xi_{j,k})}{N_k},
\end{align}
 where  $ \alpha>0$ is the constant steplength, $N_k$  is the number of sampled gradients used   at   time  $k$, and $  \xi_{j,k}, j=1, \cdots, N_k,$ denote the  \blue{independent and identically distributed (i.i.d.)}
  realizations of $\xi$.
\end{algorithm}

If  the gradient observation   noise $w_{k,N_k}$ is defined as
\begin{align} \label{def-noise}
w_{k,N_k}\triangleq \tfrac{ \sum_{j=1}^{N_k} \nabla F (x_k,\xi_{j,k})}{N_k}- \nabla f(x_k),
\end{align}  the  update \eqref{Alg1} can be rewritten as follows:
\begin{align}\label{ref_alg1}
x_{k+1}=x_k - \alpha  (\nabla f(x_k)+w_{k,N_k}).
\end{align}
 Define  $\mathcal{F}_k \triangleq \sigma\{x_0, \xi_{  j,t}, 1\leq j \leq N_{t} ,     0\leq t\leq k-1  \}.$ Then by  the definition of Algorithm
 \ref{Alg_1},  $x_k$ is adapted to $\mathcal{F}_k$. We impose the following  conditions on the  objective function  and  the gradient noise.

  \begin{assumption}\label{ass-fun}
   (i)  $f(\cdot)$ is continuously differentiable in $x\in \mathbb{R}^m$ with Lipschitz continuous gradient, i.e., there exists a constant  $L>0$ such that for any $x,x'\in \mathbb{R}^m $, $ \| \nabla f(x) -\nabla f(x')  \|  \leq L \| x-x'\|.$ (ii) $f(\cdot) $ is $\eta$-strongly convex, i.e.,
    $ ( \nabla f(x) -\nabla f(x') )^T (x-x')  \geq \eta \| x-x'\|^2$ for any $ x,x'\in \mathbb{R}^m .$
(iii)   Suppose there exists  a constant $\nu>0$ such that  for any $x\in \mathbb{R}^m$ and $\xi \in \mathbb{R}^d$,
the stochastic first-order oracle produces  $ \nabla  F (x ,\xi) $ satisfying  the following  in an a.s. sense:
$$\mathbb{E}[ \nabla  F (x ,\xi) |x]= \nabla f(x ) ~  \mbox{~and~}~  \mathbb{E}[\| \nabla  F (x ,\xi)  - \nabla f(x)\|^2|x]\leq  \nu^2 .$$
    \end{assumption}

Since  $x_k$ is adapted to $\mathcal{F}_k$ and the   samples $  \xi_{j,k}, j=1, \cdots, N_k,$  are    independent, by Assumption \ref{ass-fun}(iii), we obtain that  for any $k\geq 0$,
\begin{equation}\label{ass-bd-noise}
\begin{split}
 \mathbb{E}[w_{k,N_k}| \mathcal{F}_k]&=\tfrac{ \sum_{j=1}^{N_k}   \mathbb{E}\big[ \nabla F(x_k,\xi_{j,k}) - \nabla f(x_k) | \mathcal{F}_k\big]  }{ N_k}=0,  \quad a.s.~,
\\  \mathbb{E}[\|w_{k,N_k}\|^2| \mathcal{F}_k]&= \tfrac{\mathbb{E} \Big[  \big\| \sum_{j=1} ^{N_k }   ( \nabla F(x_k,\xi_{j,k}) - \nabla f(x_k) \big\|^2 \Big|F_k  \Big]} {N_k^2}\\
 & =\tfrac{\sum_{j=1} ^{N_k} \mathbb{E}   [  \|  \nabla F(x_k,\xi_{j,k}) - \nabla f(x_k) \|^2 |F_k   ]}{N_k^2}  \leq \tfrac{\nu^2}{N_k} ,\quad a.s. ~.
\end{split}
\end{equation}

Since $f(\cdot)$ is strongly convex,  it has a unique optimal solution denoted by  $x^*$. Then by the optimality condition, $\nabla f(x^*)=0.$ We now introduce an inequality \cite[Eqn. (2.1.24)]{nesterov2013introductory} on  $f(\cdot)$  satisfying   Assumptions \ref{ass-fun}(i) and \ref{ass-fun}(ii):
\begin{align}\label{inequ-f}
(x-y)^T(\nabla f(x) -\nabla f(y)) \geq \tfrac{\eta L\|x-y\|^2}{\eta+L} +\tfrac{\|\nabla f(x)- \nabla f(y)\|^2}{\eta+L},\quad \forall x,y\in \mathbb{R}^m.
\end{align}
We then establish a  recursion on   the mean-squared estimation  error,
for which the proof can be found in Appendix \ref{app:lem1}.
\begin{lemma}\label{Lemma_1}
Let   Algorithm  \ref{Alg_1} be applied to   \eqref{Problem1}. Suppose Assumption \ref{ass-fun} holds and $\alpha \in (0, {2\over \eta+L}]$.
Then the following relation holds for any $k\geq 0$:
\begin{align}\label{recursion1}
\mathbb{E}[\| x_{k+1}-x^*\|^2  ] &   \leq
\left(1-\tfrac{2 \alpha \eta L}{\eta+L} \right)  \mathbb{E}[\| x_k - x^*\|^2]+\tfrac{ \alpha^2   \nu^2}{N_k} .
\end{align}
\end{lemma}
 \subsection{Accelerated Gradient Method.}\label{Sec:Alg:ACC}
Nesterov's  accelerated gradient   descent method converges  with a rate $\mathcal{O}(q^k)$ where $q \triangleq 1-\sqrt{\eta  / L}$ for $\eta$-strongly convex and $L$-smooth   functions  \cite{nesterov2013introductory}, and  with a rate $ \mathcal{O}(1/k^2 )$ for merely convex functions \cite{nesterov1983method}. Nesterov proved that this is the optimal possible rate
for any method based on first-order gradients. As such, we  combine the
Nestrov's  accelerated  method with  Algorithm  \ref{Alg_1}  and   propose an
accelerated variable sample-size stochastic gradient descent algorithm
(Algorithm \ref{Alg_2})  so as to improve the  convergence rate.  Such a
scheme has been employed for smooth strongly convex~\cite{jalilzadeh2018optimal}, smooth
convex~\cite{ghadimi2016accelerated}, and nonsmooth (but smoothable) convex~\cite{jalilzadeh2018optimal} stochastic
optimization problems with associated rates of convergence given by
$\mathcal{O}(q^k)$,  $\mathcal{O}(1/k^2)$, and $\mathcal{O}(1/k)$,
respectively.   The present paper  takes a crucial towards developing  CLTs and
confidence regions for the optimal solution.
  \begin{algorithm}[H]
\caption{Variance-reduced accel. SGD} \label{Alg_2}
Given   arbitrary  initial values  $x_0=y_0\in \mathbb{R}^m, $   positive constants $\alpha,\beta, $ and a positive integer sequence $\{N_k\}_{k\geq 0}.$
Then iterate the following equations for $k\geq 0.$
\begin{subequations}
  \begin{align} y_{k+1}& =x_k -\alpha     (\nabla f(x_k)+w_{k,N_k}) ,\label{Alg21}
  \\ x_{k+1}&=y_{k+1}+\beta  (y_{k+1}-y_k),\label{Alg22}
\end{align}
\end{subequations}
where $w_{k,N_k} $ is defined as in \eqref{def-noise}, $\alpha>0$ and $\beta>0 $ are constant steplengths.
\end{algorithm}

In the following,  we establish a bound on  the  expected  sub-optimality gap   of the iterates  generated by Algorithm  \ref{Alg_2}. The  proof of Lemma~\ref{lem-acc}   motivated by  \cite{nesterov2013introductory}
is given in Appendix \ref{app-lem-acc}.  This is an important preliminary result
to be used in the rate analysis of Algorithm  \ref{Alg_2}.
 \begin{lemma}\label{lem-acc}    Let Algorithm  \ref{Alg_2}
be applied to the problem \eqref{Problem1}. Suppose Assumption \ref{ass-fun} holds and   $   \alpha \in ( 0, {1\over L}]$.
\blue{Set} $\gamma = \sqrt{ \alpha \eta } $ and  $\beta = {1 -\gamma \over 1+ \gamma}.$ Then  for any $k\geq 1,$
\begin{align}\label{lem-fy}
\mathbb{E}[f(y_{k }) ]-f^* \leq     \tfrac{(\eta+L)(1-\gamma)^{k }}{2}     \mathbb{E}[\| x_0-x^*\|^2]  +  \nu^2 \left( \alpha +\tfrac{(1-\gamma)\gamma}{2\eta}\right)  \sum_{i=0}^{k-1}  \tfrac{(1-\gamma )^i}{N_{k-1-i}}.
\end{align}
\end{lemma}
\subsection{The Heavy Ball Method.}   \label{Sec:Alg:HB}
The classical heavy ball method of Polyak \cite{polyak1964some}
takes the following  form $$x_{k+1}  =x_k -\alpha     \nabla f(x_k)
+\beta(x_k-x_{k-1}) $$ with a constant steplength $\alpha>0$ and
a momentum parameter $\beta>0. $   It was  shown in
\cite{polyak1987introduction}  that for smooth strongly convex,
the heavy ball method  generates sequences    convergent  to
the optimal  solution at a faster rate  than the    gradient descent
method. Stochastic variants of  the  heavy ball method have been
employed widely in practice (cf. \cite{sutskever2013importance,szegedy2015going,krizhevsky2012imagenet}
for applications to machine learning). Recent efforts
\cite{yang2016unified,gadat2018stochastic}  have investigated
rate  for the stochastic variant, where \cite{yang2016unified} proved
a sublinear rate $ \mathcal{O}(  1/\sqrt{k} )$  for general Lipschitz
continuous convex objectives  with bounded variance, while a rate $ \mathcal{O}(
1/k^{\beta} )$ with $\beta\in(0,1)$ was provided   in
\cite{gadat2018stochastic}   for the  case of the quadratic  strongly  convex
functions.  While convergence properties have not been   addressed  so far  variable sample-size  schemes. As such,
 we  consider a variable sample-size  variant of the stochastic heavy ball method assuming constant steplengths, in contrast with   the diminishing steplengths utilized in \cite{yang2016unified,gadat2018stochastic}, and prove  the global linear convergence.  In  the developed Algorithm \ref{Alg_3},   we   add a momentum term $\beta(x_k-x_{k-1}) $ to the variable sample-size  stochastic gradient step \eqref{Alg1} and obtain the update \eqref{Alg3}.

\begin{algorithm} [H] \caption {Variance-reduced heavy-ball SGD}\label{Alg_3}
Given an arbitrary  initial value $x_0\in \mathbb{R}^m, $  two positive constants $\alpha,\beta>0,$ and a positive  integer sequence $\{N_k\}_{k\geq 0}.$  Set $x_{-1}=x_0$. Then iterate the following equation  for $k\geq 0.$
\begin{align}\label{Alg3} x_{k+1}& =x_k -\alpha     (\nabla f(x_k)+w_{k,N_k}) +\beta(x_k-x_{k-1}),
\end{align}
where $w_{k,N_k} $ is defined as in \eqref{def-noise}, $\alpha>0$ and $\beta>0 $ are constant steplengths.
\end{algorithm}

In Section \ref{sec:CLT-geo}, we derive a geometric rate
and asymptotic normality statements for
Algorithm \ref{Alg_3}  under geometrically  increasing batch-sizes,
whereas the corresponding  results    of Algorithm \ref{Alg_3} under
polynomially increasing batch-sizes will be  established  in
Section \ref{sec:CLT-poly}. In the following,  we  give a  simple
recursion on   the  expected mean-squared error of the generated
sequence $\{x_k\}$.  \blue{The proof can be found in Appendix \ref{App-Lem-HB}.}
\begin{lemma} \label{Lem-HB}
  Suppose that Assumption \ref{ass-fun} holds and   $f(\cdot)$ is
twice continuously differentiable. Consider Algorithm  \ref{Alg_3}, where     $   \alpha \in (0,4/L)$ and  $\beta \triangleq  \max\{|1-\sqrt{\alpha \eta}|^2,|1-\sqrt{\alpha L}|^2\}<1.$  Then  the following    holds.
 \begin{align}\label{lem-hb}
 \mathbb{E} \left [ \left \| \begin{pmatrix}
  x_{k+1}-x^*
\\ x_k-x^* \end{pmatrix} \right\|^2   \right]\leq \beta   \mathbb{E} \left [ \left\| \begin{pmatrix}
 x_{k }-x^*
\\ x_{k-1}-x^* \end{pmatrix} \right \|^2 \right]+{ \alpha^2 \nu^2\over N_k}  ,\quad \forall k\geq 0.
\end{align}
\end{lemma}

 \section{Central Limit Theorems  under geometrically increasing $N_k$}  \label{sec:CLT-geo}
 In this section, we  establish CLTs for   Algorithms  \ref{Alg_1},   \ref{Alg_2},  and   \ref{Alg_3} when the number of sampled gradients,  denoted by $N_k$, increases at a geometric rate.

 \subsection{ Rate and Oracle Complexities} \label{sec:geo_rate}
\blue{Based on Lemmas \ref{Lemma_1}, \ref{lem-acc}, and \ref{Lem-HB}, we can  establish the geometric rate of convergence  along with  the iteration and oracle complexity  of the iterates generated by Algorithms~\ref{Alg_1}, ~\ref{Alg_2}, and \ref{Alg_3}.} \blue{The proofs of Propositions \ref{prp1}, \ref{prp2}, and \ref{prp-alg3} are similar to that of \cite[Theorem 4.2 and Corollary 4.7]{lei18asynchronous}. Related results for linear convergence for stochastic gradient methods can be found in~\cite{byrd12sample,friedlander12hybrid,shanbhag2015budget,jofre19variance,pasupathy18sampling} while a linear rate for the accelerated variants has been provided in~\cite{jalilzadeh2018optimal}.}

 \begin{proposition}[Rate and Oracle Complexity   for Algorithm  \ref{Alg_1}]  \label{prp1} Let Algorithm  \ref{Alg_1}
be applied to  \eqref{Problem1}. Let Assumption \ref{ass-fun} hold  and $   \alpha \in (0, {2\over \eta+L}]$.
  Suppose  that for any $k\geq 0$,  $N_{k } \triangleq \lceil \rho_1^{-(k+1)} \rceil $  with   $\rho_1 \in (q,1)$ and  $q \triangleq 1-{2 \alpha \eta L \over \eta+L}$. Then
  \begin{align}\label{recursion2}
\mathbb{E}[\| x_{k}-x^*\|^2  ]  \leq\rho_1^{k} \left(\mathbb{E}[ \| x_0 - x^*\|^2]+\tfrac{\alpha^2 \nu^2}{1- q/\rho_1}\right) ,\quad \forall k\geq 1.
\end{align}
  The  iteration and oracle  complexity for computing an $\epsilon-$solution, defined as $\mathbb{E}[\|x_k-x^*\|^2] \leq \epsilon$,  are  repectively bounded by  $\mathcal{O}\left (  \kappa    \ln(1/\epsilon)   \right)$ and   $\mathcal{O}\left ({ \kappa    / \epsilon} \right)$, where  $\kappa \triangleq {L \over \eta}$ denotes  the condition number.
 \end{proposition}
 {\bf Proof.}
 By  substituting $N_k=\lceil \rho_1^{-(k+1)} \rceil $  into \eqref{recursion1},     using $q= 1-{2 \alpha \eta L \over \eta+L} $ and $\rho_1 \in (q,1)$,  one obtains
 \begin{align*}
&\mathbb{E}[\| x_{k+1}-x^*\|^2  ]    \leq q  \mathbb{E}[\| x_k - x^*\|^2]+\alpha^2\nu^2 \rho_1^{k+1}
   \leq q^{k+1} \mathbb{E}[ \| x_0 - x^*\|^2]+\alpha^2\nu^2 \sum_{t=0}^{k} q^t \rho_1^{k+1-t}    \\
& = q^{k+1}\mathbb{E}[ \| x_0 - x^*\|^2] +\alpha^2\nu^2  \rho_1^{k+1}\sum_{t=0}^{k} ( q/\rho_1)
  \leq \left(\mathbb{E}[ \| x_0 - x^*\|^2]+{\alpha^2 \nu^2 \over 1- q/\rho_1}\right) \rho_1^{k+1}.
\end{align*}

  Suppose  we set $\alpha\triangleq{2\over \eta+L}$ and  $\rho_1 \triangleq \left({\kappa   \over \kappa+1}\right)^2 <q=\left({\kappa-1 \over \kappa+1}\right)^2$.   From  \eqref{recursion2} it follows that
$\mathbb{E}[\| x_{k }-x^*\|^2  ] \leq \epsilon $ for any $ k\geq K(\epsilon) $, where
\begin{align}\label{iter-bd} K(\epsilon) \triangleq {  \ln \left(\mathbb{E}[ \| x_0 - x^*\|^2]+{\alpha^2 \nu^2 \over 1- q/\rho_1}\right)  +\ln \left ({1/\epsilon} \right)  \over 2 \ln  \left (1+{1\over \kappa} \right)} .\end{align}
It is noticed that $ \ln (1+{1/ \kappa}) \approx 1/\kappa$ for large $\kappa$, then  the number of iterations   required  to obtain an  $\epsilon $-optimal solution in a mean-squared sense,  i.e. $\mathbb{E}[\| x  -x^*\|^2  ] \leq \epsilon$,  is $  \mathcal{O} \big(\kappa \ln \left ({1/\epsilon} \right)\big)$.
Finally, based on the iteration complexity bound \eqref{iter-bd},   we achieve the following   {oracle complexity}  bound, measured by the number of sampled gradients,  for obtaining an  $\epsilon $-optimal solution:
\begin{align*}
\sum_{k=0}^{K(\epsilon) -1 } N_k  & \leq \sum_{k=1}^{K(\epsilon)  } \rho_1^{-k}\leq \int_{1}^{K(\epsilon)+1} \rho_1^{-t} dt
 \leq { \rho_1^{-(K(\epsilon)+1)} \over \ln(1/\rho_1)} ={\mathbb{E}[ \| x_0 - x^*\|^2]+{\alpha^2 \nu^2 \over  1- q/\rho_1} \over  2 \epsilon \rho_1 \ln(1+1/\kappa)}  =\mathcal{O}\left ({ \kappa    \over \epsilon} \right) .  \qquad \Box
\end{align*}

 Suppose we set $\alpha\triangleq{2\over \eta+L}$ in Algorithm \ref{Alg_1}.     Then   by \eqref{recursion2}, $\mathbb{E}[\| x_{k}-x^*\|^2  ] =\mathcal{O}\left(\rho_1^{k}\right) $ with $ \rho_1 \in \left (\big({\kappa-1 \over \kappa+1}\big)^2,1\right).$ This matches the best known rate for gradient descent method with optimal tuning, see \cite[Proposition 12]{lessard2016analysis}.
Next, we provide similar statements for Algorithm~\ref{Alg_2}.
 \begin{proposition}[Rate and Oracle Complexity for  Algorithm  \ref{Alg_2}]  \label{prp2} Let Assumption \ref{ass-fun} hold. Consider Algorithm \ref{Alg_2},  where  $   \alpha \in ( 0, {1\over L}]$, $\gamma\triangleq \sqrt{\alpha \eta}$,   $\beta \triangleq {1 -\gamma \over 1+ \gamma},$ and $N_{k } \triangleq \lceil \rho_2^{-(k+1)} \rceil $  with   $\rho_2 \in (1-\gamma,1)$.  Then for any $k\geq 0 $:
\begin{align}\label{lem-fy2}
\mathbb{E}[f(y_{k }) ]-f^* \leq      \rho_2^k  \left( \tfrac{\eta+L}{2}    \mathbb{E}[\| x_0-x^*\|^2] +  \tfrac{\rho_2 \nu^2}{\rho_2-(1-\gamma)}\left( \alpha +\tfrac{(1-\gamma)\gamma}{2\eta}\right)  \right) ,
\end{align}
\blue{and $\mathbb{E}[\|x_k-x^*\|^2] \leq c\rho_2^k  $ for some constant $c>0$. In addition,}
 the  iteration and  oracle  complexity for obtaining  an $\epsilon-$solution  are         $\mathcal{O}\left (  \sqrt{ \kappa } \ln(1/   \epsilon) \right)$ and   $\mathcal{O}\left ({ \sqrt{ \kappa }  / \epsilon} \right)$, respectively.
\end{proposition}
{\bf Proof.}  By  using $\rho_2 \in (1-\gamma,1)$ and   substituting $N_{k }=\lceil \rho_2^{-(k+1)} \rceil $  into  Eqn. \eqref{lem-fy}, we obtain that
\begin{align*}
\mathbb{E}[f(y_{k }) ]-f^* & \leq         \rho_2^{k }   \tfrac{\eta+L}{2}    \mathbb{E}[\| x_0-x^*\|^2]   +  \nu^2 \left( \alpha +\tfrac{(1-\gamma)\gamma}{2\eta}\right)  \rho_2^k \sum_{i=0}^{k-1}  \left(\tfrac{  1-\gamma}{\rho_2} \right)^i  ,
\end{align*}
which leads to \eqref{lem-fy2} by using  the bound $ \sum_{i=0}^{k-1}  \left({  1-\gamma   \over \rho_2} \right)^i \leq {1 \over 1-{  1-\gamma   \over \rho_2}}= {\rho_2   \over \rho_2-(1-\gamma)}$.

  Because $f(x)$ is $\eta$-strongly convex and $\nabla f(x^*)=0,$  we have $f(x)-f (x^*) \geq {\eta \over 2} \|x-x^*\|^2.$ Thus, $ y_k $ generated by Algorithm  \ref{Alg_2}  satisfies
$\|y_k-x^*\|^2 \leq {2 \over \eta} (f(y_{k })  -f^* ).$
Then from \eqref{lem-fy2} it follows that   $$\mathbb{E}[\|y_k-x^*\|^2] \leq c \rho_2^k  {~~\rm with ~~} c\triangleq \tfrac{2}{\eta}\left(\tfrac{\eta+L}{2}    \mathbb{E}[\| x_0-x^*\|^2] + \tfrac{\rho_2 \nu^2}{\rho_2-(1-\gamma)}\left( \alpha +\tfrac{(1-\gamma)\gamma}{2\eta}\right)\right)  .$$  Note  by \eqref{Alg22} that  $\|x_{k}-x^*\| \leq \| (1+\beta) (y_{k}-x^*) \| +\| \beta (y_{k-1}-x^*)\| $,
 and hence
\begin{align}\label{bd-x-alg2}\mathbb{E}[\|x_{k}-x^*\|^2] \notag & \leq 2(1+\beta)^2 \mathbb{E}[\|y_k-x^*\|^2] +2\beta^2 \mathbb{E}[\|y_{k-1}-x^*\|^2]   \\
& \leq c\left( 2(1+\beta)^2+2\beta^2 \rho_2^{-1} \right)\rho_2^k .
\end{align}
Therefore,  sequences $\{x_k\} $ and $ \{y_k\}$ converge to the optimal solution  $x^*$ at a geometric rate $\mathcal{O}(\rho_2^k)$ in the mean-squared sense. Suppose we set $\alpha\triangleq{1\over L}$. Then
$\gamma = \sqrt{   \eta  \over L } ={1\over \sqrt{\kappa}}$ and $\beta=  {\sqrt{\kappa}-1\over  \sqrt{\kappa}+1}$.
Select $\rho_2 \triangleq 1-{1\over 2 \sqrt{\kappa}} $ such that $\rho_2>1-\gamma$.
Similarly to the proof of Proposition \ref{prp1}(ii), we can show that  the number of iterations   required  to obtain an  $\epsilon $-optimal solution in a mean-squared sense  is  $  \mathcal{O} \left (  {   \ln \left ({1/\epsilon} \right)  \over   \ln  \left ( 1/\rho_2 \right)} \right) = \mathcal{O} \left (   \sqrt{\kappa}    \ln \left ({1/\epsilon} \right)   \right)$  since $\ln \left(1 \over 1-{1/(2 \sqrt{\kappa})}  )\right)\approx {1\over 2 \sqrt{\kappa}} $ for large  $\kappa,$ and hence the oracle complexity $
\sum_{k=0}^{K(\epsilon) -1 } N_k  =\mathcal{O}\left ({\sqrt{ \kappa } / \epsilon} \right) .  $
 \hfill $\Box$

 Similar statements  for Algorithm  \ref{Alg_3}  are given in the next proposition.

\begin{proposition}[Rate and Oracle Complexity for Algorithm  \ref{Alg_3}]  \label{prp-alg3}
 Suppose   Assumption \ref{ass-fun} holds and that $f(\cdot)$ is  twice continuously differentiable.
Consider  Algorithm  \ref{Alg_3},  where  $   \alpha \triangleq { 4 \over (\sqrt{\eta}+\sqrt{L})^2}$,
$\beta \triangleq   \left(  { \sqrt{\kappa} -1 \over \sqrt{\kappa} +1}\right)^2,$  and  $N_{k } \triangleq \lceil \rho_3^{-(k+1)} \rceil $  with   $\rho_3 \in \left(\beta,1\right)$.  Then
 \begin{align}\label{lem-hb2}
 \mathbb{E} \left [ \left \| \begin{pmatrix}
  x_{k+1}-x^*
\\ x_k-x^* \end{pmatrix} \right\|^2 \right]\leq  \left(2   \mathbb{E}[ \| x_0 - x^*\|^2]+\tfrac{\alpha^2  \nu^2}{1-  \beta/\rho_3}\right) \rho_3^{k+1}.  \end{align}
  In addition, the  iteration and  oracle  complexity required for obtaining  an $\epsilon-$solution in a mean-squared sense are      respectively   $\mathcal{O}\left (  \sqrt{ \kappa } \ln(1/   \epsilon) \right)$ and   $\mathcal{O}\left ({ \sqrt{ \kappa }  / \epsilon} \right)$.
\end{proposition}
{\bf Proof.}
By substituting $   \alpha ={ 4 \over (\sqrt{\eta}+\sqrt{L})^2}$ into $\beta =  \max\{|1-\sqrt{\alpha \eta}|,|1-\sqrt{\alpha L}| \}^2,$   there holds     $\beta = \left(1-{2 \over \sqrt{\kappa} +1}\right)^2.$  Then the result of  Lemma \ref{Lem-HB} holds. Therefore, by using  \eqref{lem-hb},   $N_{k  }=\lceil \rho_3^{- (k+1)} \rceil $, $x_{-1}=x_0$,  and $\rho_3 \in (  \beta ,1)$,  we obtain the following sets of inequalities for any $k\geq 0$:
\begin{align*}
\mathbb{E}& \left [ \left \| \begin{pmatrix}
  x_{k+1}-x^*
&\\ x_k-x^* \end{pmatrix} \right\|^2    \right]   \leq   \beta \mathbb{E} \left [ \left \| \begin{pmatrix}
  x_{k }-x^* \\ x_{k-1}-x^* \end{pmatrix} \right\|^2    \right]+\alpha^2 \nu^2   \rho_3^{k+1}
 \\ &    \leq 2   \beta ^{ k+1 } \mathbb{E}[ \| x_0 - x^*\|^2]+\alpha^2 \nu^2  \sum_{t=0}^{k}  \beta^t \rho_3^{k+1-t} .
\end{align*}
This together with  $\sum_{t=0}^{k}  \beta^{t} \rho_3^{k+1-t} =\rho_3^{k+1}\sum_{t=0}^{k}  (\beta/\rho_3)^{t} \leq{\rho_3^{k+1}\over 1-\beta/\rho_3}$ proves \eqref{lem-hb2}.

 By \eqref{lem-hb2}, $ \mathbb{E} [\| x_k-x^* \|]\leq  c \rho_3^{k}  $ for some constant $c>0.$ Suppose $\rho_3 = \left(1-{1 \over \sqrt{\kappa} +1}\right)^2> \beta$. Akin  to the proof of   Proposition \ref{prp1}(ii), we can show that  the number of iterations   required  to obtain an  $\epsilon $-optimal solution satisfying   $\mathbb{E}[\|x_k-x^*\|^2] \leq \epsilon$   is   $  \mathcal{O} \left (  {   \ln \left ({1/\epsilon} \right)  \over   \ln  \left ( 1/\rho_3 \right)} \right) = \mathcal{O} \left (   \sqrt{\kappa}    \ln \left ({1/\epsilon} \right)   \right)$  since $\ln \left(1 \over 1-{1/(  \sqrt{\kappa}+1)}  )\right)\approx {1\over   \sqrt{\kappa}+1} $ for large  $\kappa,$ and the oracle complexity $
\sum_{k=0}^{K(\epsilon) -1 } N_k  =\mathcal{O}\left ({\sqrt{ \kappa } / \epsilon} \right) .  $
\hfill $\Box$

\begin{remark}
(i) From the  complexity  statements in Propositions \ref{prp1} and \ref{prp2},
the  dependency on the condition  number is improved from $\kappa$  (in Algorithm \ref{Alg_1})
to $\sqrt{\kappa}$  by the accelerated Algorithm \ref{Alg_2}.

\noindent (ii) The rate and   oracle  complexity for the heavy ball method  (Algorithm \ref{Alg_3}) is
 similar  to that of  Algorithm \ref{Alg_2} but {Algorithm \ref{Alg_3} requires stronger assumptions than     Algorithms~\ref{Alg_1}--\ref{Alg_2} in that $f(\cdot)$ is  twice continuously differentiable.}

\noindent (iii) Propositions \ref{prp1}, \ref{prp2}, \ref{prp-alg3} imply that a smaller $\kappa$ leads to smaller constants in the oracle complexity bounds.
\end{remark}

\begin{remark}[Convergence  in Probability]\label{rem-con-pro}
Because  the mean-squared convergence implies converges in probability,    on the basis of   Propositions \ref{prp1}-\ref{prp-alg3}, we may conclude that the sequences $\{x_k\}$   and $\{y_k\}$  generated by  Algorithms   \ref{Alg_1}-\ref{Alg_3}  converge in probability to the optimal solution $ x^* $.
\end{remark}
\vspace{-0.2in}
\subsection{Preliminary Lemmas}
  Before establishing CLTs for Algorithms \ref{Alg_1}-\ref{Alg_3},  we first introduce a preliminary   CLT on doubly-indexed random variables \cite[Lemma 3.3.1]{chen2006stochastic}. We state it as Lemma \ref{lem-CLT} (proof in~\cite[Chapter 12]{hall2014martingale}).

  \begin{lemma}\label{lem-CLT}
  Let $\xi_{kt},t= 1,\cdots, k$ be $m$-dimensional random vectors. Define
 \begin{align} \label{def-RS}
 & \mathbf{S}_{kt}  \triangleq \mathbb{E}[\xi_{kt}\xi_{kt}^T], \quad  \mathbf{R}_{kt}\triangleq \mathbb{E}[\xi_{kt}\xi_{kt}^T| \xi_{k1},\cdots, \xi_{k,t-1}],
  {\rm ~and ~} \mathbf{S}_k\triangleq \sum_{t=1}^k  \mathbf{S}_{kt} \\
& \mbox{\rm Assume that \ }   \mathbb{E}[\xi_{kt}| \xi_{k1},\cdots, \xi_{k,t-1}]=0, ~ \sup_{k\geq 1} \sum_{t=1}^k \mathbb{E}[\| \xi_{kt}\|^2]<\infty, \label{cond1}\\
   &  \lim_{k\to \infty}  \mathbf{S}_k =\mathbf{S},~ \lim_{k \to \infty} \sum_{t=1}^k \mathbb{E}[\|  \mathbf{S}_{kt}- \mathbf{R}_{kt}\|]=0 , \label{cond2}\\
 & {\rm and~ }    \lim_{k\to \infty}  \sum_{t=1}^k \mathbb{E}[\| \xi_{kt}\|^2I_{[\|\xi_{kt}\|\geq \epsilon]}]=0, \quad \forall \epsilon>0.\label{cond3}
  \end{align}
  Then   $ \sum_{t=1}^k \xi_{kt}   \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\mathbf{S}).$
  \end{lemma}

To establish  the CLTs,   we further require  the following conditions.
 \begin{assumption}\label{ass-CLT}
 (i)  $ \nabla f(x):\mathbb{R}^m \to \mathbb{R}^m$  is differentiable at $x^*$ with Hessian   matrix $ \mathbf{H}$,
 \begin{align}\label{Ass2}
 \nabla f(x) = \mathbf{H}(x-x^*)+\delta(x),   {\rm where~} \delta(x^*)=0 {\rm~and~} \delta(x)=o(\| x-x^*\|).
 \end{align}
 (ii)  The noise  sequence $\{w_{k,N_k}\} $ further satisfies:
 \begin{align}
&\lim_{k\to \infty} \Big( N_k \mathbb{E}\left[ w_{k,N_k}w_{k,N_k}^T | \mathcal{F}_{k } \right]\Big) =\lim_{k\to \infty} \Big( N_k \mathbb{E}\left[ w_{k,N_k}w_{k,N_k}^T  \right]\Big) =\mathbf{S}_0,\quad a.s., \label{Ass3}
 \\ &\label{Ass32} {\rm and~}
\lim_{r\to \infty} \sup_k \mathbb{E}\left[ \left \|\sqrt{N_{k}} w_{k,N_{k}} \right \|^2 I_{\left[\|\sqrt{N_{k}} w_{k,N_{k}}\|> r\right]}\right]=0,
 \end{align}
 where  $I_{[a>b]}=1$ if $a>b$, and $I_{[a>b]}=0$, otherwise.
  \end{assumption}

\begin{remark}  Since  the Lindeberg condition   condition   \eqref{Ass32} is  less easily verified,  we   provide a   sufficient condition for the ease of understanding.  Suppose there exists a constant $\delta>0$ and   a finite value $b>0$ such that $  \mathbb{E}[\| \nabla f(x ,\xi)  - \nabla f(x)\|^{2+\delta} ]\leq  b .$ Similarly to \eqref{ass-bd-noise}, we can  show that $\mathbb{E}\big[   \|\sqrt{N_k} w_{k,N_k}   \|^{2+\delta}\big]  \leq b $ for any $ k\geq 0.$  Therefore,

\begin{align*}
& \quad \mathbb{E}\left[ \left \|\sqrt{N_{k }} w_{k,N_k} \right \|^2 I_{\left[\|\sqrt{N_{k }} w_{k,N_k}\|> r\right]}\right] \\
 & \leq  {1\over r^{\delta}}\mathbb{E}\left[ \left \|\sqrt{N_{k }} w_{k,N_k} \right \|^{2+\delta} I_{\left[\|\sqrt{N_k } w_{k,N_k}\|> r\right]}\right]
 \leq  {1\over r^{\delta}}\mathbb{E}\left[ \left \|\sqrt{N_{k }} w_{k,N_k} \right \|^{2+\delta}  \right] \leq  {b\over r^{\delta}}.
\end{align*}
\begin{align*}
\Longrightarrow \sup_k \mathbb{E}\left[ \left \|\sqrt{N_{k }} w_{k,N_k} \right \|^2 I_{\left[\|\sqrt{N_{k }} w_{k,N_k}\|> r\right]}\right]& \leq  {b\over r^{\delta}}  \quad \Rightarrow \eqref{Ass32} {\rm~holds}.
\end{align*}
\end{remark}
In the following, we  establish  the  central limit theorem of  a    linear  recursion,  the proof of which is provided  in Appendix \ref{app-lem-CLT2}. This result will be applied in the proof of Theorems \ref{thm-CLT1}-\ref{thm-CLT1-Alg3}.

  \begin{lemma}\label{lem-CLT2}  Suppose  that  $ \mathbf{P}$ is a stable matrix (i.e., $\|\mathbf{P} \| <
\varrho$ for some $\varrho<1$),    $N_{k }=\lceil \rho^{-(k+1)} \rceil$ for any $ k\geq 0, $ and that  $\{ w_{k,N_k} \} $ satisfies Assumption \ref{ass-CLT}(ii).  Let  $\{e_k\}$ be generated by the following  recursion
\begin{align}\label{recursion-e}
e_{k+1}& =    \mathbf{P}  e_k -\alpha \rho^{-(k+1)/2} \mathbf{G}  w_{k,N_k}+\zeta_{k+1},~\mathbb{E}[\| e_0\|^2]<\infty,
\end{align}
where $\zeta_{k+1}    \xlongrightarrow [k \rightarrow \infty]{P} 0 $.
Then  the following holds:
 \begin{align*}
  &\alpha^{-1}   e_k\xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}),
\quad {\rm where} \quad   \bm{\Sigma}  \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}  \mathbf{P} ^{t}  \mathbf{G} \mathbf{S}_0    \mathbf{G}^T \left( \mathbf{P}^t \right)^T  .\end{align*}
\end{lemma}

  \subsection{Central Limit Theorems }\label{sec:geo_CLT}
Based on  Proposition \ref{prp1} and  Lemma \ref{lem-CLT2}, using Assumption \ref{ass-CLT}, we are now  ready to formulate  the central limit theorem for  Algorithm \ref{Alg_1} with geometrically increasing  batch-sizes.

 \begin{theorem}[CLT of Algorithm \ref{Alg_1} with Geometrically increasing $N_k$]\label{thm-CLT1} Let Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold. Consider Algorithm \ref{Alg_1},  where  $   \alpha \in (0, {2\over \eta+L}]$.    Set $q\triangleq 1-{2 \alpha \eta L \over \eta+L} $, $N_{k } \triangleq \lceil \rho_1^{-(k+1)} \rceil $ with $\rho_1 \in (q,1)$,
 and  $\mathbf{P}_1  \triangleq  \rho_1^{-1/2}(\mathbf{I}_m-\alpha  \mathbf{H})$.  Then
 \begin{align}  \label{Geo-CLT}
  &\alpha^{-1} \rho_1^{-k/2}( x_{k}-x^*)  \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_1),
~{\rm where}~  \bm{\Sigma}_1 \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}  \mathbf{P}_1 ^{t} \mathbf{S}_0    \mathbf{P}_1 ^{t}  .\end{align}
 \end{theorem}
 {\bf Proof. }   By definition of $\delta(x)$ in \eqref{Ass2},  Eqn.  \eqref{ref_alg1} can be rewritten as:
 \begin{align}\label{recursion3}
x_{k+1}-x^*& =x_k-x^* -\alpha  \mathbf{H}(x_k-x^*)-\alpha  (\nabla f(x_k)- \mathbf{H}(x_k-x^*)+w_{k,N_k}) \notag \\
& = (\mathbf{I}_m-\alpha  \mathbf{H})(x_k-x^*) -\alpha  (\delta(x_k)+w_{k,N_k}).
\end{align}
 Note by  \eqref{Ass2} if $\mathbf{D}_k\triangleq { \delta(x_k)(x_k-x^*)^T \over \| x_k-x^*\|^2}$,
 \begin{align}\label{def-delta} \delta(x_k)=\left( \delta(x_k) {(x_k-x^*)^T \over \| x_k-x^*\|^2} \right) (x_k-x^*)=\mathbf{D}_k (x_k-x^*).
 \end{align}
 Since  $ \delta(x_k)=o(\| x_k-x^*\|)$ and    $ x_k  \xlongrightarrow [k \rightarrow \infty]{P}  x^* $ (since   $ x_k  \xlongrightarrow [k \rightarrow \infty]{L^2}  x^* $),  we obtain that  $\mathbf{D}_k \xlongrightarrow [k \rightarrow \infty]{P} 0$. Define $e_{k} \triangleq \rho_1^{-k/2}(x_{k}-x^*)$. By multiplying  both sides of  \eqref{recursion3} by  $\rho_1^{-(k+1)/2}$,   using \eqref{def-delta}  and the definition $  \mathbf{P}_1  \triangleq  \rho_1^{-1/2}(\mathbf{I}_m-\alpha  \mathbf{H}) $,
 \begin{align}\label{recursion4}
e_{k+1}& =   \rho_1^{-1/2}(\mathbf{I}_m-\alpha  \mathbf{H})\rho_1^{-k/2}(x_{k}-x^*) -\alpha \rho_1^{-(k+1)/2} (\mathbf{D}_k (x_k-x^*)+w_{k,N_k}) \notag
\\& =    \mathbf{P}_1   e_k -\alpha \rho_1^{-(k+1)/2} w_{k,N_k}+\zeta_{k+1},
\end{align}
where $\zeta_{k+1} \triangleq - \alpha \rho_1^{-(k+1)/2} \mathbf{D}_k (x_k-x^*).$  From \eqref{recursion2},  it follows  that for any $k\geq 0,$ $\mathbb{E}[\| x_{k}-x^*\|^2 /\rho_1^k ]  \leq \mathbb{E}[\|x_0-x^*\|^2]+{\alpha^2 \nu^2 \over  1-q/\rho_1}  .$  Thus,  the random variable $ e_k={\rho_1^{-k/2} (x_{k}-x^* )}$ is almost surely (a.s.)  finite.
This is because $\mathbb{E}[X^2]<\infty$ implies that $P(\omega: X(\omega)=\infty)=0.$
Otherwise,   $P(\omega: X(\omega)=\infty)>0$    leads to a contradiction $\mathbb{E}[X^2]=\infty$. Hence
\begin{align}\label{limit-zeta}
\zeta_{k+1} = -\alpha \rho_1^{-1/2} \mathbf{D}_k  e_k  \xlongrightarrow [k \rightarrow \infty]{P} 0 {\rm~~by~~}  \mathbf{D}_k \xlongrightarrow [k \rightarrow \infty]{P} 0.
\end{align}

Since $\mathbf{H}$ is the Hessian matrix of    $f(x)$ at $x=x^*,$ by Assumptions \ref{ass-fun}(i) and  \ref{ass-fun}(ii),
we conclude  that $\mathbf{H}$ has eigenvalues $\lambda_1, \cdots, \lambda_m$ that satisfy $0<\eta \leq \lambda_m\leq \lambda_{m-1} \leq \dots\leq \lambda_2\leq \lambda_1\leq L.$  Then  $\|  \mathbf{I}_m-\alpha  \mathbf{H}  \|_2 \leq   \max\{ |1-\alpha \eta|,  |1-\alpha L| \} .$ We first show that  \blue{$ |1-\alpha \eta| \geq  |1-\alpha L| $}. It is easily seen that  $|1-\alpha \eta| \leq  |1-\alpha L|$  when $\alpha\in (0,1/L]$.
While for any $\alpha\in [1/L,{2\over \eta +L}]$, $| 1-\alpha \eta|=1-\alpha \eta,  |1-\alpha L|=\alpha L -1 $, and hence $|1-\alpha \eta| - |1-\alpha L| =2-\alpha (\eta+L) \geq 0$. Then for any     $\alpha \in (0,{2\over \eta +L}] $,
$$\|  \mathbf{I}_m-\alpha  \mathbf{H}  \|  \leq   | 1-\alpha \eta | =1-\alpha \eta \leq \sqrt{1-\tfrac{2 \alpha \eta L}{\eta+L}}  .
  $$
 This is because  $ \alpha \leq {2\over \eta+L}\Rightarrow {2L \over \eta+L}+\alpha\eta-2\leq 0\Rightarrow {2\alpha\eta L \over \eta+L}+(\alpha\eta)^2-2\alpha\eta\leq 0 \Rightarrow (1-\alpha \eta)^2 \leq  1-{2 \alpha \eta L \over \eta+L}.  $  Since $q=1-{2 \alpha \eta L \over \eta+L}$ and   $\rho_1 \in (q,1)$,  we have that
\begin{align}\label{bd-P}
\|  \mathbf{P}_1 \| =\rho_1^{-1/2}\| (\mathbf{I}_m-\alpha  \mathbf{H}) \|  \leq     (q/\rho_1)^{1/2}<1.
\end{align}
Then  by  invoking  Lemma \ref{lem-CLT2},   setting $
\mathbf{G}=\mathbf{I}_m$ and $ \mathbf{P}= \mathbf{P}_1$
(symmetric),  we conclude that the sequence $\alpha^{-1}e_k$
generated by  \eqref{recursion4}  converges in distribution to a
normally distributed random variable with  zero mean and  covariance
$\bm{\Sigma}_1$ defined as in \eqref{Geo-CLT}. The result follows by  $
e_k=\rho_1^{-k/2}( x_{k}-x^*)$.   \hfill $\Box$

 \begin{remark}
Suppose we set  $\alpha=\tfrac{2}{\eta+L}$ and $\rho_1 \triangleq \left(\tfrac{\kappa}{\kappa+1}\right)^2 $in  Algorithm \ref{Alg_1}.    If $V \sim  N(0,\mathbf{I}_m)$, it follows from \eqref{Geo-CLT}  that
\begin{align}
& \tfrac{ \eta+L}{2} \left(  \tfrac{ \kappa+1}{\kappa }\right)^k( x_{k}-x^*)  \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_1)  \notag \\  \implies & x_k  \overset {\tiny D}{\approx } x^* +\tfrac{2}{\eta+L }\left(1- \tfrac{ 1}{\kappa+1 }\right)^k\bm{\Sigma}_1^{1/2} V \mbox{~for large~}k. \label{Geo-CLT2}
\end{align}
This result has several implications.  \\
(i)  Expression  \eqref{Geo-CLT2} implies that the sequence $ \{x_{k}\} $ converges in distribution to the optimal solution $x^* $ at    rate 	   $(\tfrac{\kappa}{\kappa+1 })^k $, and  $\{x_k\}$ is asymptotically normally distributed  for large $k$.   This provides the possibility of assessing   confidence regions of the estimate  from the normal distribution.  \\
(ii) The estimation error for large $k$ depends on the     structure of the studied problem  (including $\eta, L$, and the Hessian matrix $\mathbf{H}$) and probability distribution  of the gradient   noise,   measured through  the coupling   matrix $\bm{\Sigma}_1.$  Thus, the problem's difficulty is largely characterized by  the covariance matrix $\bm{\Sigma}_1.$
\end{remark}

Based on the CLT established in Theorem \ref{thm-CLT1}, we
proceed to  use the Delta method~\cite{shapiro2009lectures} to    derive the
asymptotic distribution  of the gradient  map  $\nabla f (\cdot) $
and  the objective function value $f(\cdot)$  based on  the estimation  sequence  $\{x(k)\}$.

\begin{corollary} Consider Algorithm \ref{Alg_1} and suppose that all  conditions listed in Theorem \ref{thm-CLT1} hold. Then

\noindent (i) $\alpha^{-1} \rho_1^{-k/2}  \nabla f (x_{k})    \xlongrightarrow [k \rightarrow \infty]{d}   N(0,  \mathbf{H} \bm{\Sigma}_1 \mathbf{H} ).$

\noindent  (ii) $\alpha^{-2} \rho_1^{-k}  (f (x_{k}) -f(x^*))   \xlongrightarrow [k \rightarrow \infty]{d}   {1\over 2 }N(0,  \bm{\Sigma}_1  )^T \mathbf{H}N(0,  \bm{\Sigma}_1  ).$
\end{corollary}
{\bf Proof.} (i)  Note  by Assumption \ref{ass-CLT}(i) that  the gradient
mapping  $ \nabla f(x):\mathbb{R}^m \to \mathbb{R}^m$  is differentiable at
$x^*$ with Hessian   matrix $ \mathbf{H}$.  By  using \eqref{Geo-CLT} and the
Delta theorem~\cite[Eqn. (7.182)]{shapiro2009lectures}, we have that
$$\alpha^{-1} \rho_1^{-k/2} \big( \nabla f (x_{k})- \nabla f (x^*)\big)
\xlongrightarrow [k \rightarrow \infty]{d}   \mathbf{H}  N(0,  \bm{\Sigma}_1
).  $$
Then  by  the fact that $\nabla f (x^*)=0 $  and $ \mathbf{H} $ is symmetric,  the assertion (i) holds.

 (ii) By  using \eqref{Geo-CLT}, $\nabla f (x^*)=0,$  and the second-order Delta theorem~\cite[Theorem 7.70]{shapiro2009lectures},   one obtains
$$\alpha^{-2} \rho_1^{-k}  \big(f (x_{k}) -f(x^*)\big)   \xlongrightarrow [k \rightarrow \infty]{d}   {1\over 2 } f_{u}^{''}(x^*)  {~\rm with ~}u=N(0,  \bm{\Sigma}_1  ), $$
  where $f_{u}^{''}(x^*) $  denotes  the second order directional  derivative at   $x^*$ along the direction $u$.   Since   $f:\mathbb{R}^m \to \mathbb{R} $ is twice continuously differentiable, $f_{u}^{''}(x^*)=u^T \mathbf{H} u$.
Then   result (ii) follows from $u=N(0,  \bm{\Sigma}_1  )$.
\hfill $\Box$

 Akin  to Theorem \ref{thm-CLT1},  based on Proposition \ref{prp2}  and
Lemma \ref{lem-CLT2},   we    may also establish  the CLT for
Algorithm \ref{Alg_2} with geometrically increasing batch-sizes.

 \begin{theorem}[CLT for Algorithm \ref{Alg_2} with  Geometrically Increasing $N_k$]\label{thm-CLT1-Alg2} Suppose Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold. Consider Algorithm \ref{Alg_2}, where  $   \alpha \in (0, {1\over  L}]$,  $\gamma \triangleq \sqrt{ \alpha \eta } $,   $\beta \triangleq {1 -\gamma \over 1+ \gamma}$, and  $N_{k }\triangleq \lceil \rho_2^{-(k+1)} \rceil $  with   $\rho_2 \in (1-\gamma,1)$. Then
 \begin{align*}
  &\alpha^{-1} \rho_2^{-k/2} \begin{pmatrix}
y_{k }-x^* \\
y_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_2) {~\rm with~}\bm{\Sigma}_2 \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}  \mathbf{P}_2 ^{t} \begin{pmatrix} \mathbf{S}_0       & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}  \left( \mathbf{P}_2 ^{t} \right)^T ,
\\ &    \quad {\rm and}~   \mathbf{P}_2 \triangleq  \rho_2^{-1/2}\begin{pmatrix}
(1+\beta) (\mathbf{I}_m-\alpha  \mathbf{H})  &~ -\beta (\mathbf{I}_m-\alpha  \mathbf{H})
\\  \mathbf{I}_m   &~\mathbf{0}_m \end{pmatrix}.\end{align*}
 \end{theorem}
{\bf Proof.} Similarly to \eqref{recursion3}, by using  \eqref{Ass2},  Eqn.  \eqref{Alg21}   can be rewritten as: \begin{align}\label{Alg21-y}
y_{k+1}-x^*&  = (\mathbf{I}_m-\alpha  \mathbf{H})(x_k-x^*) -\alpha  (\delta(x_k)+w_{k,N_k}).
\end{align}
Note by  \eqref{Alg22}  that $ x_{k }-x^* = (1+\beta) (y_{k }-x^*)-\beta (y_{k-1}-x^*).$
This together with \eqref{Alg21-y} produces
  \begin{align*}
& y_{k+1}-x^* = (\mathbf{I}_m-\alpha  \mathbf{H})( (1+\beta) (y_{k }-x^*)-\beta (y_{k-1}-x^*)) -\alpha  (\delta(x_k)+w_{k,N_k})
\\&=  (1+\beta) (\mathbf{I}_m-\alpha  \mathbf{H})  (y_{k }-x^*)-\beta (\mathbf{I}_m-\alpha  \mathbf{H}) (y_{k-1}-x^*) -\alpha  (\delta(x_k)+w_{k,N_k})    .
\end{align*}
Define $z_{k+1} \triangleq     \begin{pmatrix}
y_{k+1}-x^* \\
y_{k}-x^* \end{pmatrix} $ and $\mathbf{H}_2 \triangleq  \begin{pmatrix}
(1+\beta) (\mathbf{I}_m-\alpha  \mathbf{H})  & ~-\beta (\mathbf{I}_m-\alpha  \mathbf{H})
\\  \mathbf{I}_m   &  ~\mathbf{0}_m \end{pmatrix} $. Then based on  the above equation, we obtain the following recursion:
\begin{align}  \label{re-z}
z_{k+1}&  =\mathbf{H}_2 z_k   -  \alpha \begin{pmatrix}   \delta(x_k)+w_{k,N_k}    \\ 0\end{pmatrix}  .
\end{align}

The eigenvalue decomposition of $\mathbf{H}$ is given by $\mathbf{H}=\mathbf{U}\bm{\Lambda} \mathbf{U}^T$,  where $\bm{\Lambda} \triangleq {\rm diag}\{\lambda_1, \lambda_2,\cdots, \lambda_m\}$ and $\mathbf{U}$ is orthogonal. This allows us to rewrite  \blue{$\mathbf{H}_2$ as}
\begin{align*} \mathbf{H}_2 = \begin{pmatrix}\mathbf{U} & 0
\\  0  & \mathbf{U}
  \end{pmatrix}    \begin{pmatrix}
(1+\beta) (\mathbf{I}_m-\alpha   \bm{\Lambda})  & ~-\beta (\mathbf{I}_m-\alpha  \bm{\Lambda})
\\  \mathbf{I}_m   &  ~0 \end{pmatrix}  \begin{pmatrix}\mathbf{U} & 0
\\  0  & \mathbf{U}
  \end{pmatrix}^T .  \end{align*}
As a result,  the eigenvalues of  the matrix $ \mathbf{H}_2  $ are the roots of the following  characteristic  equations for $i = 1, \cdots, m$.
$$  \begin{vmatrix} v-   (1+\beta) (1-\alpha    \lambda_i)  & ~ \beta (1-\alpha   \lambda_i )
\\ -1   & ~ v \end{vmatrix} =v^2-(1+\beta) (1-\alpha    \lambda_i) v+   \beta (1-\alpha   \lambda_i ) =0.$$
By   $\gamma^2=\alpha \eta$ and $\beta = {1 -\gamma \over 1+ \gamma},$  we show that  the discriminant of the above  quadratic equation is nonpositive for each $i=1,\cdots,m$.
\begin{align*}
\Delta_i & =(1+\beta)^2 (1-\alpha    \lambda_i)^2 -4  \beta (1-\alpha   \lambda_i )=4(1-\alpha   \lambda_i ) \left(  {(1-\alpha   \lambda_i )  \over (1+ \gamma)^2}-{1 -\gamma \over 1+ \gamma} \right )
\\& ={ 4(1-\alpha   \lambda_i ) \over (1+ \gamma)^2}  \left(   \gamma^2-\alpha   \lambda_i \right )=
{ 4\alpha (1-\alpha   \lambda_i ) \over (1+ \gamma)^2}  \left(  \eta-  \lambda_i \right ) \leq 0,
\end{align*}
where the last inequality holds by the fact that for any $   i=1, \dots,m, $  $1-\alpha   \lambda_i  \geq  1-  \lambda_i /L\geq 0$ and  $   \eta -\lambda_i \leq 0. $  Consequently, the above equation has two complex roots ${ (1+\beta) (1-\alpha    \lambda_i) \pm \mathbf{i} \sqrt{-\Delta_i } \over 2}$,  where $\mathbf{i}$ denotes the imaginary part.  Thus,  the magnitude  of the roots is  $
\sqrt{ \left(\tfrac{ (1+\beta) (1-\alpha    \lambda_i)}{2} \right)^2 -\tfrac{\Delta_i}{4}}= \sqrt{\beta (1-\alpha   \lambda_i ) }  .$  Since $\gamma^2=\alpha \eta$, $\beta ={1 -\gamma \over 1+ \gamma},$  $\rho_2 \in (1-\gamma,1)$, and $\lambda_i  \in [\eta, L]$,  we obtain that
\begin{align} \label{sp-h}
\| \mathbf{H}_2 \| & =\max_i  \sqrt{\beta (1-\alpha   \lambda_i ) } \leq \sqrt{\beta (1-\alpha  \eta) }
=\sqrt{\tfrac{(1 -\gamma)(1-\gamma^2)}{1+ \gamma}  } \\
	\notag & =1-\gamma <\rho_2 <\rho_2^{1/2}.
\end{align}
Define $\varepsilon_k \triangleq \rho_2^{-k/2} z_k.$  By multiplying  both sides of \eqref{re-z} by $\rho_2^{-(k+1)/2}$,    we obtain the following recursion \blue{with $\varsigma_{k+1} \triangleq -\alpha  \rho_2^{-(k+1)/2}  \delta(x_k ):$}
\begin{equation*}
\begin{split}
\varepsilon_{k+1}&  =\rho_2^{-1/2}\mathbf{H}_2 \varepsilon_k   -  \alpha  \rho_2^{-(k+1)/2} \begin{pmatrix}   \delta(x_k)+w_{k,N_k}    \\ 0\end{pmatrix}
\\&  = \mathbf{P}_2 \varepsilon_k   -  \alpha  \rho_2^{-(k+1)/2} \begin{pmatrix}    w_{k,N_k}    \\ 0\end{pmatrix}+\begin{pmatrix}   \varsigma_{k+1}    \\ 0\end{pmatrix} .
\end{split}
\end{equation*}
 By \blue{Proposition \ref{prp2}},   $\mathbb{E}[\rho_2^{-k}\|x_{k}-x^*\|^2]  $   is uniformly bounded.
 Hence    $ \rho_2^{-k/2}  x_k  $ is a.s. finite.   Recall from  Remark  \ref{rem-con-pro}   that $ x_k  \xlongrightarrow [k \rightarrow \infty]{P}  x^*$. Similarly to \eqref{limit-zeta},  we can   show  that  $\varsigma_{k+1} =-\alpha  \rho_2^{-(k+1)/2}  \delta(x_k ) \xlongrightarrow [k \rightarrow \infty]{P} 0 $.
 Because $\|\mathbf{P}_2\|   = \rho_2^{-1/2} \| \mathbf{H}_2\|<1$  by \eqref{sp-h},   by
  \blue{setting $\mathbf{G}= \begin{pmatrix} \mathbf{I}_m    \\ \mathbf{0}_m\end{pmatrix} $ and} using Lemma \ref{lem-CLT2}, we obtain that $ \alpha^{-1}  \varepsilon_k \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_2) $
with $\bm{\Sigma}_2  \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}  \mathbf{P}_2 ^{t}  \mathbf{G} \mathbf{S}_0    \mathbf{G}^T \left( \mathbf{P}_2^t \right)^T   .$
Then by $\mathbf{G} \mathbf{S}_0    \mathbf{G}^T= \begin{pmatrix} \mathbf{S}_0      & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}$ and $\varepsilon_k=\rho_2^{-k/2} z_k=\rho_2^{-k/2}  \begin{pmatrix}
y_{k }-x^* \\ y_{k-1}-x^* \end{pmatrix}$, we prove the result.
\hfill  $\Box$

\begin{remark}
 Suppose $\alpha={1\over L}$ in Algorithm \ref{Alg_2} and   $\rho_2
=1-{1\over a\sqrt{\kappa}} $ for some $a>1$. Then it may be noticed that
$1-\gamma=1-{1\over \sqrt{\kappa}}<\rho_2.$ Then Theorem \ref{thm-CLT1-Alg2}
holds, and hence  we obtain the following with  $V \sim  N(0,\mathbf{I}_{2m}).$  \begin{align*}
&L  \left(1-{1\over a\sqrt{\kappa}}\right)^{-k/2} \begin{pmatrix}
y_{k }-x^* \\
y_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_2)   \\  \Longrightarrow \quad & \begin{pmatrix}
y_{k }  \\
y_{k-1} \end{pmatrix} \overset {\tiny D}{\approx } \begin{pmatrix}
 x^* \\
 x^* \end{pmatrix} +{1\over L }\left( 1-{1\over  a\sqrt{\kappa}}\right)^{k/2}\bm{\Sigma}_2^{1/2} V \mbox{~for large~}k
\end{align*}
\end{remark}

 In a similar fashion,   we establish  the central limit theorem for the heavy ball method (Algorithm \ref{Alg_3}) with geometrically increasing batch-sizes  based on Proposition~\ref{prp-alg3}  and   Lemma \ref{lem-CLT2}.

 \begin{theorem}[CLT  for  Algorithm \ref{Alg_3} with Geometrically Increasing $N_k$]\label{thm-CLT1-Alg3} Suppose Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold and  $f(\cdot)$ is  twice continuously differentiable.
  Let Algorithm \ref{Alg_3} be applied to the problem \eqref{Problem1}, where  $   \alpha \triangleq { 4 \over (\sqrt{\eta}+\sqrt{L})^2}$,  $\beta \triangleq   \left(1-{2 \over \sqrt{\kappa} +1}\right)^2,$  and   $N_{k }\triangleq \lceil \rho_3^{- (k+1)} \rceil $  with   $\rho_3 \in \left( \beta,1\right)$ for any $k\geq 0.$ Then
 \begin{align*}
  &\alpha^{-1} \rho_3^{-k/2} \begin{pmatrix}
x_{k }-x^* \\
x_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_3) {~\rm with~}\bm{\Sigma}_3 \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}  \mathbf{P}_3 ^{t} \begin{pmatrix} \mathbf{S}_0       & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}  \left( \mathbf{P}_3 ^{t} \right)^T,\end{align*}
where $\mathbf{P}_3 \triangleq  \rho_3^{-1/2}\begin{pmatrix} (1+\beta) \mathbf{I}_m -\alpha    \mathbf{H} &~-\beta  \mathbf{I}_m\\
\mathbf{I}_m &~\mathbf{0}_m
\end{pmatrix}.$
 \end{theorem}
{\bf Proof.}  From  \eqref{Alg3} and Assumption \ref{ass-CLT}(a),  we observe  that
\begin{align*}
x_{k+1}-x^*& =x_k-x^* +\beta(x_k-x_{k-1}) -\alpha    \mathbf{H}   ( x_k -x^*) -\alpha \delta(x_k)  -\alpha  w_{k,N_k}
\end{align*}Then we may rewrite the recursion in a  matrix form as follows:
\begin{align*}
\begin{pmatrix}
 x_{k+1}-x^*
\\ x_k-x^* \end{pmatrix} =\begin{pmatrix} (1+\beta) \mathbf{I}_m -\alpha   \mathbf{H} &-\beta  \mathbf{I}_m\\
\mathbf{I}_m &\mathbf{0}_m
\end{pmatrix} \begin{pmatrix}
 x_{k }-x^*
\\ x_{k-1}-x^* \end{pmatrix}  -\alpha\begin{pmatrix}   \delta(x_k)+  w_{k,N_k} \\ 0\end{pmatrix} .
\end{align*}
Define $\varepsilon_k \triangleq \rho_3^{-k/2} \begin{pmatrix}
 x_{k }-x^*
\\ x_{k-1}-x^* \end{pmatrix} .$  By multiplying   the above equation with $\rho_3^{-(k+1)/2}$,  one obtains
\blue{the following with $\varsigma_{k+1} \triangleq -\alpha  \rho_3^{-(k+1)/2}  \delta(x_k )$:}
\begin{equation*}
\begin{split}
\varepsilon_{k+1}&  =  \mathbf{P}_3   \varepsilon_k   -  \alpha  \rho_3^{-(k+1)/2} \begin{pmatrix} \mathbf{I}_m  \\ \mathbf{0}_m\end{pmatrix}   w_{k,N_k}  +\begin{pmatrix}   \varsigma_{k+1}    \\ 0\end{pmatrix}  .
\end{split}
\end{equation*}
 We recall that  $ x_k  \xlongrightarrow [k \rightarrow \infty]{P}  x^* {\rm ~and ~}   \rho_3^{-k/2}  x_k \mbox{ ~is a.s. finite}  $  by  Proposition \ref{prp-alg3}. Akin to \eqref{limit-zeta},  we  have  that  $\varsigma_{k+1}   \xlongrightarrow [k \rightarrow \infty]{P} 0 $.   Since  $\|\mathbf{P}_3\|= \rho_3^{-1/2}  \sqrt{\beta}<1$  by \eqref{eigen-alg3} ,   the result follows by invoking Lemma \ref{lem-CLT2}.
  \hfill $\Box$

\begin{remark}
 Suppose we set $\rho_3 = \left(1-{1 \over \sqrt{\kappa} +1}\right)^2$ in Algorithm \ref{Alg_3},  then $\rho_3>\beta $ and  Theorem \ref{thm-CLT1-Alg3} holds. If $V \sim  N(0,\mathbf{I}_{2m})$, then
\begin{align*}
& \tfrac{   (\sqrt{\eta}+\sqrt{L})^2}{4}  \left(1-\tfrac{1}{\sqrt{\kappa} +1}\right)^{-k} \begin{pmatrix}
x_{k }-x^* \\
x_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_3)   \\  \Longrightarrow \quad & \begin{pmatrix}
x_{k }  \\
x_{k-1} \end{pmatrix} \overset {\tiny D}{\approx } \begin{pmatrix}
 x^* \\
 x^* \end{pmatrix} +\tfrac{ 4}{(\sqrt{\eta}+\sqrt{L})^2} \left( 1-\tfrac{1}{\sqrt{\kappa} +1}\right)^{k }\bm{\Sigma}_3^{1/2} V \mbox{~for large~}k.
\end{align*}
\end{remark}

   \section{Central Limit Theorems on Polynomial  Batch-size}\label{sec:CLT-poly}

 There are many settings where a geometric increase in $N_k$
is impractical. For instance, this is true when the
generation of a sampled gradient is computationally
expensive; an example of this commonly arises in simulation
optimization problems in the context of large manufacturing
or queueing simulations. To this end, we consider the use of
polynomial increases in sample-size, an avenue that allows
for more gentle growth and proceed to   investigate  the
rate  of  convergence and the associated central limit
statements for Algorithms~\ref{Alg_1}--\ref{Alg_3}. \us{Polynomial rates have been studied in~\cite{lei18asynchronous} as well as~\cite{pasupathy18sampling} but we remain unaware of rate and complexity statements as well as the ensuing CLTs in accelerated and heavy-ball regimes.}
\vspace{-0.2in}
 \subsection{ Rate of Convergence}
 We first recall  some  preliminary results     that find utility in  the
 rate analysis of the proposed algorithms with  the polynomially increasing batch-sizes.
\begin{lemma} \label{lem-recur}  (i)   \cite[Eqn. (17)]{lei18game}  For any $q\in (0,1)$ and $v>0$, there  holds   \[\sum_{t=1}^{k  } q^{k  -t} t^{-v}    \leq  q^{k  } {e^{2v}q^{-1}-1 \over 1-q}+{2  k^{-v}\over q \ln(1/q) }.\]

\noindent  (ii) \cite[Lemma 4]{lei18game}  For any $q\in (0,1)$ and $v>0$,  $q^x \leq c_{q,v}  {x^{-v}}$ for all $x >0$ where $c_{q,v} \triangleq   e^{- v  }{\left(\tfrac{v}{ \ln(1/q)}\right)^{ v  }}.$
\end{lemma}
Based on Lemmas \ref{Lemma_1}-\ref{Lem-HB}  in Section \ref{Sec:Alg} and Lemma \ref{lem-recur}, we can establish polynomial  rates of convergence of the iterates generated by the three proposed methods.
 \begin{proposition}[Rate statement for Algorithm  \ref{Alg_1} under polynomially increasing $N_k$]  \label{prp3}    Suppose Assumption \ref{ass-fun} holds  and that   $N_{k }   \triangleq \lceil (k+1)^v\rceil $ for some $v>0$. Consider  Algorithm  \ref{Alg_1} with $  \alpha \in (0, {2\over \eta+L}]$.   Define  $q\triangleq 1-{2 \alpha \eta L \over \eta+L} $ and $c_{q,v} \triangleq   e^{- v  }{\left(\tfrac{v}{ \ln(1/q)}\right)^{ v  }}.$   Then   the following holds  for any $k\geq 1$.
\[ \mathbb{E}[\| x_{k}-x^*\|^2  ]  \leq  \left(  c_{q,v} \mathbb{E}[ \| x_0 - x^*\|^2] + \tfrac{\alpha^2\nu^2 c_{q,v}(e^{2v}q^{-1}-1)}{1-q}+\tfrac{2\alpha^2\nu^2}{q \ln(1/q) }\right) k^{-v}.\]
 \end{proposition}
{\bf Proof.} By substituting   $N_k   \triangleq \lceil (k+1)^v\rceil $ into \eqref{recursion1}, using   $q= 1-{2 \alpha \eta L \over \eta+L} $ and Lemma \ref{lem-recur}(i), one obtains
 \begin{align*}
 \mathbb{E}[\| x_{k}-x^*\|^2  ]  & \leq  q \mathbb{E}[\| x_{k-1} - x^*\|^2]+\alpha^2\nu^2 k^{-v} \\& = q^{k} \mathbb{E}[ \| x_0 - x^*\|^2] +\alpha^2\nu^2 \sum_{m=1}^{k} q^{k-m} m^{-v} \\ & =  q^k \mathbb{E}[ \| x_0 - x^*\|^2]  +\alpha^2\nu^2  \left(q^{k } \tfrac{e^{2v}q^{-1}-1}{1-q}+\tfrac{2   k^{-v}}{q \ln(1/q) }\right) ,
\end{align*}
which together with  Lemma \ref{lem-recur}(ii) proves  the lemma.
\hfill $\Box$

Thus, the   number of iterations   required    to  obtain an $\epsilon $-accurate optimal solution in the mean-squared sense (i.e. $\mathbb{E}[\| x  -x^*\|^2  ] \leq \epsilon$) is   $  \mathcal{O}  \big( \left ({1/\epsilon} \right)^{1/v} \big)$.

 \begin{proposition}[Rate statement	 for  Algorithm  \ref{Alg_2} under polynomially increasing $N_k$]  \label{prp4}Let Algorithm  \ref{Alg_2}
be applied to   \eqref{Problem1}, where $N_{k }   \triangleq \lceil (k+1)^v\rceil $ with some   $v>0$. Suppose Assumption \ref{ass-fun} holds and   $   \alpha \in ( 0, {1\over L}]$.
Define  $\gamma\triangleq \sqrt{\alpha \eta}$ and  $\beta \triangleq {1 -\gamma \over 1+ \gamma}.$  Then there exists a constant $C (v) >0$ such that \[
   \mathbb{E}[f(y_{k }) ]-f^* \leq   C (v)  k^{-v},\quad \forall k\geq 1 .\]
\end{proposition}
{\bf Proof.}  It is noticed   by $t = k-i$ and Lemma \ref{lem-recur}(i),
\begin{align}\label{bd-rk}
  \sum_{i=0}^{k-1}  (1-\gamma )^i (k-i)^{-v} &   =    \sum_{t=1}^{k } (1-\gamma )^{k-t} t^{-v}\notag \\
	&  \leq      \tfrac{(1-\gamma )^{k  }(e^{2v} (1-\gamma )^{-1}-1)}{\gamma }+\tfrac{2  k^{-v}}{(1-\gamma ) \ln(1/ (1-\gamma )) } .
\end{align}
By substituting $N_{k }=\lceil (k+1)^v \rceil $  into  Eqn. \eqref{lem-fy},    we obtain that
\begin{align*}
\mathbb{E}[f(y_{k }) ]-f^* & \leq          \tfrac{(1-\gamma)^{k }(\eta+L)\mathbb{E}[\| x_0-x^*\|^2]}{2}       +  \nu^2  \left( \alpha +\tfrac{(1-\gamma)\gamma}{2\eta}\right)  \sum_{i=0}^{k-1}  (1-\gamma )^i (k-i)^{-v} .
\end{align*}
This  together   with   Lemma \ref{lem-recur}(ii) and \eqref{bd-rk}   proves the  required result.
 \hfill $\Box$

\begin{remark} \label{rem-poly}   Since $f(x)-f (x^*) \geq {\eta \over 2} \|x-x^*\|^2,$   we have that     $\mathbb{E}[\|y_k-x^*\|^2]= \mathcal{O}(k^{-v}).$  Note  by \eqref{Alg22} that   $$\mathbb{E}[\|x_{k}-x^*\|^2] \leq 2(1+\beta)^2\mathbb{E}[\|y_k-x^*\|^2] +2\beta^2 \mathbb{E}[\|y_{k-1}-x^*\|^2]  =  \mathcal{O}(k^{-v}).$$
Thus,  the mean-squared error of Algorithm  \ref{Alg_2} also  displays  the polynomial rate of convergence similar to that shown in Proposition \ref{prp3} for Algorithm  \ref{Alg_1}. Because   the mean-squared convergence implies converges in probability,    the sequences   $\{x_k\}$ and $\{y_k\}$ generated by Algorithm   \ref{Alg_2} satisfy    $ x_k  \xlongrightarrow [k \rightarrow \infty]{P}  x^* $ and $ y_k  \xlongrightarrow [k \rightarrow \infty]{P}  x^* $ when the conditions of Proposition \ref{prp4} hold.
\end{remark}

 {
\begin{proposition}[Rate Statement for   Algorithm  \ref{Alg_3} under polyomially increasing $N_k$]  \label{prp2-alg3}
Suppose  Assumption \ref{ass-fun} holds and  $f(\cdot)$ is  twice continuously differentiable.  Consider Algorithm  \ref{Alg_3}, where $   \alpha \triangleq \tfrac{ 4}{(\sqrt{\eta}+\sqrt{L})^2}$,  $\beta \triangleq    ( \tfrac{\sqrt{\kappa} -1}{\sqrt{\kappa} +1} )^2 $,  and  $N_{k }   \triangleq \lceil (k+1)^v\rceil $ with  $v>0$. Then   there exists a constant $C(v)>0$ such that
\[ \mathbb{E} \left [ \left \| \begin{pmatrix}
  x_{k+1}-x^*
\\ x_k-x^* \end{pmatrix} \right\|^2  \right]\leq  C(v) (k+1)^{-v},\quad \forall k\geq 0.  \]
\end{proposition}
{\bf Proof.}
By substituting $N_k=\lceil (k+1)^v\rceil $  into \eqref{lem-hb} and using  Lemma \ref{lem-recur}(i),  we obtain that
 \begin{align*}
\mathbb{E}& \left [ \left \| \begin{pmatrix}
  x_{k+1}-x^*
&\\ x_k-x^* \end{pmatrix} \right\|^2    \right]   \leq   \beta \mathbb{E} \left [ \left \| \begin{pmatrix}
  x_{k }-x^*
\\ x_{k-1}-x^* \end{pmatrix} \right\|^2    \right]+\alpha^2 \nu^2  (k+1)^{-v}
 \\
&    \leq  2   \beta ^{k+1} \mathbb{E}[ \| x_0 - x^*\|^2]+\alpha^2 \nu^2  \sum_{m=1}^{k+1}  \beta^{k+1-m }   m^{-v} .
\\&  \leq  2  \beta ^{ k+1  } \mathbb{E}[ \| x_0 - x^*\|^2]+\alpha \nu  \left(\beta^{ k+1} {e^{2 v} \beta^{-1}-1 \over 1-\beta}+{    2(k+1)^{-v}\over  \beta  \ln(1/\beta) }\right) ,
\end{align*}
which together with  Lemma \ref{lem-recur}(ii) proves  the lemma.
\hfill $\Box$

\subsection{Central Limit Theorems under polynomially increasing $N_k$}
We now prove a preliminary lemma  that will be used in establishing the CLT.

\begin{lemma}\label{lem6}    For any $k\geq 1$, define $\mathbf{A}_0=\mathbf{A}$ and $ \mathbf{A}_k \triangleq   {\left(k+1 \over k\right)}^{v/2} \mathbf{A} $.
Suppose there exists $\varrho\in (0,1)$ such that  $\| \mathbf{A}\|_2 \leq \varrho.$   Define $\bm{\Phi}_{k,j} \triangleq \mathbf{A}_k \cdots \mathbf{A}_j$ with $\bm{\Phi}_{j,j+1} \triangleq \mathbf{I}_m$.    Then for any $a> 0,$  the following holds:
  \begin{align}\label{limit-Phi}
  \sup_{k\geq 1}\sum_{t=0}^k \|\bm{\Phi}_{k,t+1}\|^a \leq   c_{\varrho ^{a } , av/2}{e^{av}\varrho ^{-a  }-1 \over 1-\varrho ^{a } }+{2   \over a \varrho ^{a } \ln(1/\varrho  ) } .
  \end{align}
  \end{lemma}
{\bf Proof.}
 By  the definition  of  $\bm{\Phi}_{k,j} $ and $ \mathbf{A}_k   $, we have the  following bound:
 $$\bm{\Phi}_{k,t+1}  = {\left(k+1 \over t+1\right)}^{v/2} \mathbf{A}^{k -t} \quad \Longrightarrow \|\bm{\Phi}_{k,t+1}\| \leq {\left(k+1 \over t+1\right)}^{v/2} \varrho^{k-t}.$$
By defining  $\tilde{q}\triangleq \varrho ^{a},\tilde{v}\triangleq av/2,$  and using Lemma \ref{lem-recur}(i),
we obtain that for any $k\geq 1 :$
 \begin{align*}
& \sum_{t=0}^{k} \|\bm{\Phi}_{k,t+1}\|^a
\leq   \sum_{t=0}^{k}  {\left(k+1 \over t+1\right)}^{av/2} \varrho ^{a(k+1-(t+1)) }\leq
    \left(k+1 \right)^{\tilde{v}} \sum_{t=1}^{k+1}  t^{-\tilde{v}}  \tilde{q} ^{k+1-t}
 \\&  \leq   \left(k+1 \right)^{\tilde{v}}  \left( \tilde{q}^{k+1 } {e^{2\tilde{v}}\tilde{q}^{-1}-1 \over 1-\tilde{q}}+{2   (k+1)^{-\tilde{v}}\over \tilde{q} \ln(1/\tilde{q}) } \right)
  \\& \leq   \left(k+1 \right)^{\tilde{v}} \tilde{q}^{k+1 } {e^{2\tilde{v}}\tilde{q}^{-1}-1 \over 1-\tilde{q}}+{2   \over \tilde{q} \ln(1/\tilde{q}) }
  \leq   c_{\tilde{q},\tilde{v}}{e^{2\tilde{v}}\tilde{q}^{-1}-1 \over 1-\tilde{q}}+{2   \over \tilde{q} \ln(1/\tilde{q}) }. \end{align*}
Therefore,  the result holds.
\hfill $\Box$

Next, we  establish  the  asymptotic normality of  a     time-varying linear  recursion and   provide   the proof  in Appendix \ref{app-lem-CLT3}.  This result will be applied in  {proving} Theorems \ref{thm-CLT3}-\ref{thm-CLT2-Alg3}.

\begin{lemma}\label{lem-CLT3}  Suppose  that    $\| \mathbf{A}\| \leq \varrho $  for some $\varrho\in (0,1)$,    $N_k=\lceil ( k+1)^v \rceil$ with $v>0$, and $\{ w_{k,N_k} \} $ satisfies Assumption \ref{ass-CLT}(ii).  Let the sequence  $\{e_k\}$ be generated by the following:
\begin{align}\label{poly-recursion4}
e_{k+1}& =  \mathbf{A}_k    e_k -\alpha(k+1)^{v/2}  \mathbf{G}w_{k ,N_k }+ \zeta_{k+1}, ~~\mathbb{E} [\|e_0\|^2]<\infty,
\end{align}
where  $\mathbf{A}_0=\mathbf{A}$, $ \mathbf{A}_k \triangleq   {\left(k+1 \over k\right)}^{v/2} \mathbf{A}$  for any $k\geq 1 $, and $\zeta_{k+1}    \xlongrightarrow [k \rightarrow \infty]{P} 0 $.
Then
 \begin{align*}
  &\alpha^{-1}   e_k\xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}),
~ {\rm where} ~  \bm{\Sigma}  \triangleq \lim_{k\to \infty}\sum_{t=1}^{k} {\left(k  \over  t\right)}^{v}  \mathbf{A} ^{k-t}  \mathbf{G} \mathbf{S}_0    \mathbf{G}^T \left( \mathbf{A}^T \right)^{k-t}  .\end{align*}
\end{lemma}

Based on  Proposition \ref{prp3} and  Lemma \ref{lem-CLT3},  by using
Assumption \ref{ass-CLT}, we are now  ready to derive the associated  central
limit theorem  for Algorithm \ref{Alg_1} with polynomially increasing
batch-sizes.

 \begin{theorem}[CLT  of Algorithm \ref{Alg_1} with  Polynomially increasing    Batch-sizes]\label{thm-CLT3} Suppose Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold.  Consider Algorithm \ref{Alg_1}, where   $ \alpha \in (0, {2\over \eta+L}]$ and  $N_{k }   \triangleq \lceil (k+1)^v\rceil $ with some $v>0$. Define $\mathbf{A} \triangleq  \mathbf{I}_m-\alpha  \mathbf{H}$. Then
 \[\alpha^{-1} k^{v/2}( x_{k}-x^*)  \xlongrightarrow [k \rightarrow \infty]{d}   N(0,  \bm{\Sigma}_4) ~
{\rm with} ~ \bm{\Sigma}_4\triangleq \lim_{k\to \infty}\sum_{t=1}^k  {\left(k  \over  t\right)}^{v} \mathbf{A}^{k-t}  \mathbf{S}_0 \mathbf{A}^{k-t}  . \]
 \end{theorem}
 {\bf Proof. }    We begin by noting that  $x_k\xlongrightarrow [k \rightarrow \infty]{P} x^*$  and
 \eqref{recursion3}   holds.
Define $e_0\triangleq x_0-x^* $, $e_{k} \triangleq k^{v/2}(x_{k}-x^*)$ for any $k\geq 1$, and $\zeta_{k+1} \triangleq -\alpha(k+1)^{v/2} \delta (x_k )$ for any $k\geq 0$. By multiplying  both sides of  \eqref{recursion3} with $(k+1)^{v/2}$,   and using   $ \mathbf{A}_k =  {\left(k+1 \over k\right)}^{v/2}(\mathbf{I}_m-\alpha  \mathbf{H}) $, we obtain the following recursion:
 \begin{align} \label{def-ek}
e_{k+1}& =  {\left(k+1 \over k\right)}^{v/2} (\mathbf{I}_m-\alpha  \mathbf{H}) k^{v/2}(x_{k}-x^*) -\alpha(k+1)^{v/2}(\delta (x_k )+w_{k,N_k}) \notag
\\& =   \mathbf{A}_k  e_k - \alpha(k+1)^{v/2}w_{k,N_k}+\zeta_{k+1} , \quad \forall k\geq 1  .
\end{align}
 By setting  $k=0$ in Eqn.  \eqref{recursion3},     and using $ \mathbf{A}_0=\mathbf{A}=\mathbf{I}_m-\alpha  \mathbf{H}$,  we see that
 \begin{align*} e_1=x_1-x^*=\mathbf{A}_0(x_0-x^*)-\alpha  w_{1,N_1}-\alpha  \delta (x_0).
 \end{align*}
Then by $\zeta_1=-\alpha  \delta (x_0), $  we see that the
\eqref{def-ek} also  holds for $k=0.$ Hence the recursion \eqref{def-ek}  holds for any $k\geq 0.$
From  \eqref{bd-P} it follows  that  the symmetric matrix $\mathbf{A}$ satisfies
 $$\|\mathbf{A}\|_2 =\| \mathbf{I}_m-\alpha  \mathbf{H}  \|_2= \rho^{1/2} \| \mathbf{P}_1\|_2 \leq q^{1/2} = \left(  1-{2 \alpha \eta L \over \eta+L}\right)^{1/2}\triangleq \varrho.$$

We conclude from  Proposition \ref{prp3}    that $\mathbb{E}[k^v\| x_{k}-x^*\|^2   ]  $ is uniformly bounded by some constant.  Thus,  the random variable $e_{k} = k^{v/2}(x_{k}-x^*)$ is a.s.  finite. Then  by recalling that  $\delta(x)=o(\| x-x^*\|)$  and $x_k \xrightarrow[k \to \infty]{P} x^*$,   we obtain
\begin{align}\label{limit-zeta2}
\zeta_{k+1} =    -\alpha {\left (1+{1 \over k}\right)}^{v/2} \left( \delta(x_k) {(x_k-x^*)^T \over \| x_k-x^*\|^2} \right) e_k  \xlongrightarrow [k \rightarrow \infty]{P} 0 .
\end{align}
 Therefore,  by using Lemma \ref{lem-CLT3}  with
$\mathbf{G}=\mathbf{I}_m, $ we conclude that \blue{
$\alpha^{-1}e_k \xrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_4) $.}
 By \blue{the definition  $  e_k=k^{v/2}( x_{k}-x^*)$,}  we obtain the result.  \hfill $\Box$

Similarly to Theorem \ref{thm-CLT3},  based on Proposition \ref{prp4}
and   Lemma \ref{lem-CLT3},   we    now  establish  the central
limit theorem for Algorithm \ref{Alg_2}  under the assumption of
 polynomially increasing  batch-sizes.

 \begin{theorem}[CLT of Algorithm \ref{Alg_2} under Polynomially increasing  Batch-sizes]\label{thm-CLT2-Alg2} Suppose Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold.  Let Algorithm \ref{Alg_2} be applied to  \eqref{Problem1}, where  $   \alpha \in (0, {1\over  L}]$ and $N_{k }=\lceil (k+1)^v\rceil $  with   some $v>0$.  Define $\gamma \triangleq \sqrt{ \alpha \eta } $, $\beta \triangleq {1 -\gamma \over 1+ \gamma},$  and $\mathbf{H}_2 \triangleq  \begin{pmatrix}
(1+\beta) (\mathbf{I}_m-\alpha  \mathbf{H})  & ~-\beta (\mathbf{I}_m-\alpha  \mathbf{H})
\\  \mathbf{I}_m   &~\mathbf{0}_m \end{pmatrix}.$ Then
 \begin{align*}
  &\alpha^{-1} k^{v/2} \begin{pmatrix}
y_{k }-x^* \\
y_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0, \bm{\Sigma}_5) ,
 {\rm~where~}\bm{\Sigma}_5 \triangleq \lim_{k\to \infty} \sum_{t=0}^{k}   {\left(k  \over  t\right)}^{v}   \mathbf{H}_2 ^{k-t} \begin{pmatrix} \mathbf{S}_0       & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}  \left( \mathbf{H}_2 ^T\right)^{k-t}  .\end{align*}
 \end{theorem}
{\bf Proof.}     Define $z_{k+1} \triangleq     \begin{pmatrix}
y_{k+1}-x^* \\
y_{k}-x^* \end{pmatrix} $, $\varepsilon_0 \triangleq  z_0$, and $\varepsilon_k \triangleq k^{v/2} z_k$ for any $  k\geq 1.$  Therefore, by defining $ \varsigma_{k+1} \triangleq     -  \alpha   (k+1)^{v/2}   \delta(x_k) $ for any $k\geq 0,$ and  multiplying  both sides of  \eqref{re-z} by  $(k+1)^{v/2}$,     we obtain that for any $k\geq 1:$
\begin{equation*}
\begin{split}
\varepsilon_{k+1}&  = {\left(k+1 \over k\right)}^{v/2}\mathbf{H}_2 \varepsilon_k   -  \alpha  (k+1)^{v/2} \begin{pmatrix}  \mathbf{I}_m  \\ \mathbf{0}_m\end{pmatrix}   w_{k ,N_k }  +  \begin{pmatrix}  \varsigma_{k+1}  \\ 0 \end{pmatrix}  .
\end{split}
\end{equation*}
By defining    $\mathbf{G}\triangleq  \begin{pmatrix} \mathbf{I}_m    \\ \mathbf{0}_m\end{pmatrix}$, $\mathbf{A}_0  \triangleq \mathbf{H}_2 ,   $  and $ \mathbf{A}_k \triangleq   {\left(k+1 \over k\right)}^{v/2}  \mathbf{H}_2 $, there holds
\begin{equation}  \label{thm4-vare}
\begin{split}
\varepsilon_{k+1}&  = \mathbf{A}_k \varepsilon_k   -  \alpha  (k+1)^{v/2}  \mathbf{G}  w_{k ,N_k } +\begin{pmatrix}  \varsigma_{k+1}  \\ 0 \end{pmatrix},\quad \forall k\geq 1.
\end{split}
\end{equation}
From  \eqref{re-z} it is  seen  that    \eqref{thm4-vare}    holds for $k=0 $ as well.
Thus,  the recursion \eqref{thm4-vare} holds for any $k\geq 0.$

 Note from  Remark \ref{rem-poly} that  $\mathbb{E}[\| k^{v/2}(x_{k}-x^*) \|^2]$ is uniformly bounded. This implies that   $x_k\xlongrightarrow [k \rightarrow \infty]{P} x^*$  and that the  random variable $ k^{v/2}(x_{k}-x^*)$ is a.s.  finite. As shown in  \eqref{limit-zeta2},  we {may} show that  $\varsigma_{k+1}   \xlongrightarrow [k \rightarrow \infty]{P} 0 $. Then by     using Lemma \ref{lem-CLT3} and  \eqref{sp-h}, we obtain that
\[\alpha^{-1}  \varepsilon_k \xlongrightarrow [k \rightarrow \infty]{d}   N(0, \bm{\Sigma}_5),
 {\rm~where~}   \bm{\Sigma}_5    \triangleq \lim_{k\to \infty}\sum_{t=0}^{k} {\left(k  \over  t\right)}^{v}  \mathbf{H}_2 ^{k-t}  \mathbf{G} \mathbf{S}_0    \mathbf{G}^T \left( \mathbf{H}_2^T \right)^{k-t}  .\]
Then by the fact that $\mathbf{G} \mathbf{S}_0    \mathbf{G}^T= \begin{pmatrix} \mathbf{S}_0      & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}$ and $\varepsilon_k=k^{v/2}  \begin{pmatrix}
y_{k }-x^* \\ y_{k-1}-x^* \end{pmatrix}$,  we prove the result.
\hfill  $\Box$

 \begin{remark}
(i) Suppose we set $\alpha={2\over  L+\eta}$ in  the unaccelerated Algorithm \ref{Alg_1}.  Then $q=1-{2 \alpha \eta L \over \eta+L}=\left({\kappa-1\over \kappa+1} \right)^2$ and  by \eqref{bd-P} it is seen that  the matrix $\mathbf{A}$, defined in Theorem \ref{thm-CLT3}, satisfies  $\|  \mathbf{A}  \|_2  \leq      q ^{1/2}={\kappa-1\over \kappa+1}=1-{2\over \kappa+1}.$\\
 (ii) Suppose $\alpha={1\over  L}$ in Algorithm \ref{Alg_2}.  Then   $\gamma=1/\sqrt{\kappa} $ and by \eqref{sp-h}  we know that $\mathbf{H}_2$, defined in Theorem \ref{thm-CLT1-Alg2}, satisfies $\| \mathbf{H}_2 \|_2  =1-\gamma =1-{1\over \sqrt{\kappa}}.$\\
(iii) From Theorems \ref{thm-CLT3} and \ref{thm-CLT2-Alg2}, we note
that both  the unaccelerated gradient method  (Algorithm \ref{Alg_1})
and its    accelerated  counterpart (Algorithm \ref{Alg_2}) have convergence rates given by  $\|x_k-x^*\|=\mathcal{O}(k^{-v/2})$;  however the   accelerated scheme  has a smaller constant than its unaccelerated counterpart in settings where $\kappa$ is large,   i.e.  $\| \mathbf{H}_2 \|_2<\| \mathbf{A} \|_2$ for
large $\kappa.$ \end{remark}

 Akin  to the proof of Theorem \ref{thm-CLT1-Alg2}, based on Proposition \ref{prp-alg3}  and   Lemma \ref{lem-CLT2},   we    now  state  the central limit theorem for  Algorithm \ref{Alg_3} with polynomially  increasing batch-size (proof omitted).

 \begin{theorem}[CLT of Algorithm \ref{Alg_3} with  polynomially increasing $N_k$]\label{thm-CLT2-Alg3} Suppose Assumptions \ref{ass-fun}  and \ref{ass-CLT} hold, and that  $f(\cdot)$ is  twice continuously differentiable.
Consider  Algorithm \ref{Alg_3} , where  $   \alpha \triangleq { 4 \over (\sqrt{\eta}+\sqrt{L})^2}$,  $\beta \triangleq   \left({\sqrt{\kappa} -1 \over \sqrt{\kappa} +1}\right)^2,$  and   $N_{k }=\lceil (k+1)^v\rceil $  with   $ v>0 $. Set $\mathbf{H}_3 \triangleq   \begin{pmatrix} (1+\beta) \mathbf{I}_m -\alpha    \mathbf{H} &~-\beta  \mathbf{I}_m\\ \mathbf{I}_m &~\mathbf{0}_m \end{pmatrix}.$ Then
 \begin{align*}
    &\alpha^{-1} k^{v/2} \begin{pmatrix}
x_{k }-x^* \\
x_{k-1}-x^* \end{pmatrix} \xlongrightarrow [k \rightarrow \infty]{d}   N(0, \bm{\Sigma}_6),
~{ \rm where~} \bm{\Sigma}_6 \triangleq \lim_{k\to \infty}\sum_{t=0}^{k}   {\left(k  \over  t\right)}^{v}   \mathbf{H}_3 ^{k-t} \begin{pmatrix} \mathbf{S}_0       & \mathbf{0}_m   \\ \mathbf{0}_m &\mathbf{0}_m\end{pmatrix}  \left( \mathbf{H}_3 ^T\right)^{k-t}.\end{align*}
 \end{theorem}

\section{Confidence Regions of the Optimal Solution} \label{sec:CI}

 A crucial motivation for developing CLTs lies in developing
confidence statements. In this section,  we  proceed to construct  the
confidence regions  for  the optimal solution $x^*$. Note that the limiting covariance matrix is dependent on
the Hessian at the solution, which is unavailable. Furthermore, we do
not have a  consistent estimate of the covariance matrix.  Yet, in
the absence of such an estimate, we proceed to develop rigorous
confidence statements, adopting an   approach developed in
\cite{hsieh2002recent}.

 Since Algorithms \ref{Alg_1}, \ref{Alg_2}, and  \ref{Alg_3} lead to  similar central limit results,
we show how to construct  confidence regions  merely for the sequence
$\{x_k\}$ generated by Algorithm \ref{Alg_1}.  The  simulation
framework  in \cite{hsieh2002recent}  lies in generating $n$
independent  replications of Algorithm \ref{Alg_1}, leading to $n$ copies of the  random iterate  $ x_k $,  denoted by
 $x_{1k},\cdots, x_{nk}.$
Then the sample mean  and  the covariance estimator    are respectively defined as
\begin{align}\label{est-covariance}
\bar{\mathbf{x}}_k={1\over n}\sum_{i=1}^n x_{ik} \mbox{ ~and } \quad \mathbf{S}_k={1\over n-1}\sum_{i=1}^n (x_{ik}-\bar{\mathbf{x}}_k) (x_{ik}-\bar{\mathbf{x}}_k)^T.
\end{align}
  Based on  \cite[Theorem 4 and Corollary 1]{bodnar2016singular} and  \cite[Theorem 2]{hsieh2002recent}, we   may prove  the following result.

\begin{proposition}\label{cf-clt}
Consider Algorithm \ref{Alg_1}. Suppose the conditions of Theorem  \ref{thm-CLT1}  hold and  $n\geq m+1$. Then \\
(i) $\sqrt{n} \alpha^{-1} \rho_1^{-k/2}( \bar{\mathbf{x}}_k-x^*)  \xlongrightarrow [k \rightarrow \infty]{d}   N(0,\bm{\Sigma}_1).$ \\
(ii)  $n ( \bar{\mathbf{x}}_k-x^*)^T\mathbf{S}_k^{-1}   ( \bar{\mathbf{x}}_k-x^*)  \xlongrightarrow [k \rightarrow \infty]{d} {m(n-1) \over n-m}F(m,n-m), $ where $F(m,n-m)$ denotes the  $F$-distribution  with $(m,n-m)$ degrees of freedom.
\end{proposition}
{\bf Proof.} Based on Theorem \ref{thm-CLT1} it is seen that
\begin{align}\label{clt} e_{ik} \triangleq \alpha^{-1} \rho_1^{-k/2}( x_{ik}-x^*)  \xlongrightarrow [k \rightarrow \infty]{d}   Y_i \sim N(0,\bm{\Sigma}_1),\quad  i=1,\cdots, n,\end{align}
where $Y_1,\cdots, Y_n$   are  $n$   i.i.d.  random vectors with distribution $N(0,\bm{\Sigma}_1)$. Define
\begin{align*}
\bar{ \mathbf{e}}  \triangleq{1\over n} \sum_{i=1}^n  Y_i {\rm~and~}
\mathcal{\bm{S}}   \triangleq  {1\over n-1}  \sum_{i=1}^n\big (Y_i-
\bar{ \mathbf{e}}  \big)\big (Y_i- \bar{ \mathbf{e}}  \big)^T  .
\end{align*} Recall from \cite[Theorem 4]{bodnar2016singular}  that  $ \bar{\mathbf{e}} $ and  $\mathcal{\bm{S}} $ are independently distributed,  and
\begin{align}\label{est-es}  \bar{ \mathbf{e}}  \sim N(0,\bm{\Sigma}_1/n) {\rm~~and ~~}(n-1) \mathcal{\bm{S}}   \sim \mathcal{\bm{W}}_m(\bm{\Sigma}_1,n-1),
\end{align}
where $ \mathcal{\bm{W}}_m(\bm{\Sigma}_1,n-1)$ denotes the  $m$-dimensional Wishart distribution with $n-1$  degrees of freedom and the matrix parameter   $\bm{\Sigma}_1.$  Because $\bm{\Sigma}_1$ is invertible and $n\geq m+1$, the random  matrix  $\mathcal{\bm{S}} $ is
almost surely invertible.

 (i)   Denote by    $\mathbf{e}=col\{e_1,\cdots,e_n\} \triangleq (e_1^T,\cdots, e_n^T)^T \in \mathbb{R}^{mn}$ with $e_i \in \mathbb{R}^m,~i=1,\cdots, n$.  Note that    $g_1(\mathbf{e} )\triangleq  {1\over n}\sum_{i=1}^n e_i  $   is a continuous function.

 Since $e_{ik},~i=1,\cdots,m$ are mutually indepedent,    from \eqref{clt}  it follows that
\begin{align}\label{lim-concated-e}
col\{e_{1k},\cdots,e_{nk}\} \xlongrightarrow [k \rightarrow \infty]{d}   col\{ Y_1,\cdots, Y_n\}  .
\end{align} Then  by the continuous mapping theorem \cite[Theorem 1.14]{dasgupta2008asymptotic} and \eqref{est-es},
\begin{align}\label{clt2} &\sqrt{n} \alpha^{-1} \rho_1^{-k/2}( \bar{\mathbf{x}}_k-x^*)   =
\sqrt{n} {1\over n} \sum_{i=1}^n e_{ik}\notag
\\&  \xlongrightarrow [k \rightarrow \infty]{d} \sqrt{n} {1\over n} \sum_{i=1}^n  Y_i =\sqrt{n} \bar{ \mathbf{e}}    \sim N(0,\bm{\Sigma}_1).
\end{align}

(ii)    Note that  $g_2(\mathbf{e})   \triangleq    {1\over n-1}  \sum_{i=1}^n\left ( e_i-{1\over n}\sum_{i=1}^n e_i \right)\left ( e_i-{1\over n}\sum_{i=1}^n e_i \right)^T  $ is a continuous function.  Then by the continuous mapping theorem \cite[Theorem 1.14]{dasgupta2008asymptotic},
   \eqref{est-es}, and \eqref{lim-concated-e},  we obtain  that
\begin{equation} \label{clt3} \begin{split} &(n-1) \alpha^{-2} \rho_1^{-k }  \mathbf{S}_k = \sum_{i=1}^n \left ( e_{ik}-{1\over n}\sum_{i=1}^n e_{ik} \right)    \left ( e_{ik}-{1\over n}\sum_{i=1}^n e_{ik} \right) ^T\xlongrightarrow [k \rightarrow \infty]{d}
\\& \sum_{i=1}^n\left (Y_i-{1\over n}\sum_{i=1}^nY_i  \right)\left (Y_i-{1\over n}\sum_{i=1}^nY_i  \right)^T    = (n-1) \mathcal{\bm{S}}   \sim \mathcal{\bm{W}}_m(\bm{\Sigma}_1,n-1).
\end{split}\end{equation}
Because $rank(\bm{\Sigma}_1)=m,$   $ \bar{\mathbf{e}} \in \mathbb{R}^m$ and  $\mathcal{\bm{S}} $ are independently distributed with $ \bar{\mathbf{e}}^T\bm{\Sigma}_1$ being non-zero with probability one,  from    \cite[Corollary 1]{bodnar2016singular} it follows that
\begin{align} \label{chi-square}
 { \bar{ \mathbf{e}}^T\bm{\Sigma}_1^{-1}   \bar{ \mathbf{e}}   \over \bar{ \mathbf{e}} ^T \left( (n-1) \mathcal{\bm{S}}  \right)^{-1}  \bar{ \mathbf{e}}}\sim  \chi ^2(n-m)
\end{align}
 is  independent of $ \bar{ \mathbf{e}}  $ with  $ \chi ^2(n-m)$ denoting the   chi-squared distribution with $n-m$ degrees of freedom.  Because the matrix inverse functional is continuous in
a neighborhood of any non-singular matrix, and  $ \sum_{i=1}^n\left ( Y_i-{1\over n}\sum_{i=1}^n Y_i \right)\left ( Y_i-{1\over n}\sum_{i=1}^n Y_i \right)^T$  is almost surely invertible from  \eqref{est-es}, we conclude that  $  \big(g_2(\mathbf{e}) \big)^{-1}  $  is almost surely continuous in a neighborhood of   $col\{ Y_1,\cdots, Y_n\}.$   Hence,  $g(\mathbf{e}) =g_1(\mathbf{e})^T \big(g_2(\mathbf{e}) \big)^{-1} g_1(\mathbf{e})$  is almost surely continuous in a neighborhood of   $col\{Y_1,\cdots, Y_n\}.$  Since   $rank(\bm{\Sigma}_1)=m$ and $\sqrt{n} \bar{ \mathbf{e}}    \sim N(0,\bm{\Sigma}_1) $  by \eqref{clt2}, we obtain $(\sqrt{n} \bar{ \mathbf{e}})^T\bm{\Sigma}_1^{-1} (\sqrt{n} \bar{ \mathbf{e}}) \sim  \chi^2(m) . $ Therefore,    by the continuous mapping theorem \cite[Theorem 1.14]{dasgupta2008asymptotic}, \eqref{clt2}, and \eqref{clt3},    we have
\begin{align*} & n ( \bar{\mathbf{x}}_k-x^*)^T\mathbf{S}_k^{-1}   ( \bar{\mathbf{x}}_k-x^*)  \\&
=n \left({\sum_{i=1}^n e_{ik}\over n}\right)^T
\Big({1\over n-1} \sum_{i=1}^n \big ( e_{ik}-{1\over n}\sum_{i=1}^n e_{ik} \big)    \big ( e_{ik}-{1\over n}\sum_{i=1}^n e_{ik} \big) ^T\Big)^{-1} {\sum_{i=1}^n e_{ik}\over n}
 \\& \xlongrightarrow [k \rightarrow \infty]{d} n(n-1) \bar{ \mathbf{e}}^T \left( (n-1) \mathcal{\bm{S}}  \right)^{-1}  \bar{ \mathbf{e}} \sim   (n-1){ n\bar{ \mathbf{e}}^T\bm{\Sigma}_1^{-1}   \bar{ \mathbf{e}}   \over   { \bar{ \mathbf{e}}^T\bm{\Sigma}_1^{-1}   \bar{ \mathbf{e}}   \over \bar{ \mathbf{e}} ^T \big( (n-1) \mathcal{\bm{S}}  \big)^{-1}  \bar{ \mathbf{e}}}}    \sim {m(n-1) \over n-m} F(n,n-m),
\end{align*}
where the last one holds because $ { \bar{ \mathbf{e}}^T\bm{\Sigma}_1^{-1}   \bar{ \mathbf{e}}   \over \bar{ \mathbf{e}} ^T \left( (n-1) \mathcal{\bm{S}}  \right)^{-1}  \bar{ \mathbf{e}}} \sim  \chi ^2(n-m) $ and $n\bar{ \mathbf{e}}^T\bm{\Sigma}_1^{-1}   \bar{ \mathbf{e}} \sim  \chi^2(m)$ are independent,  and $F(d_1, d_2) $ arises as the ratio of two appropriately scaled chi-squared variates \cite{degroot2012probability}. \footnote{ If the random variables $U_1\sim \chi^2(d_1)$ and $U_2\sim \chi^2(d_2)$  are independent, then ${U_1 /d_1 \over U_2 /d_2} \sim F(d_1, d_2).$}
\hfill $\Box$

Proposition \ref{cf-clt}  can be used to construct the   confidence region  of the optimal solution.   Define
\begin{align}\label{def-conf-region}
X_{mk}(z) \triangleq \left \{ x:  n ( \bar{\mathbf{x}}_k-x )^T\mathbf{S}_k^{-1}   ( \bar{\mathbf{x}}_k-x ) \leq  {m(n-1) \over n-m}z \right \},
\end{align}
where $z$ is selected such that $\mathbb{P}(F(m,n-m) \leq z) \geq 1-\delta$ with  some $\delta\in(0,1)$.   Then we have the  following corollary.

\begin{corollary} (\cite[Proposition 3]{hsieh2002recent}) Consider Algorithm \ref{Alg_1}. Suppose   $n\geq m+1$ and the conditions of Theorem  \ref{thm-CLT1}  hold.  Then the  confidence region    $X_{mk}(z) $   defined in  \eqref{def-conf-region}  is asymptotically correct, i.e.,
$$\lim_{k \rightarrow \infty}\mathbb{P}\big(x^*\in X_{mk}(z) \big) = 1-\delta.$$
\end{corollary}
The above result asserts that  the  estimated confidence region $X_{mk}(z) $     asymptotically covers the  optimal solution  $x^*$  with probability  $100(1-\delta)\%$. The approach is easily implementable  because it   merely  requires $n$ independent   replicatications of  Algorithm \ref{Alg_1},
 while without  requiring  a consistent estimator  of the covariance  matrix of the stationary normal distribution.
The   confidence regions of  Algorithm \ref{Alg_2} and Algorithm \ref{Alg_3}  can be constructed in a similar way.

\section{Numerical Simulations} \label{sec:simu}

   In this section, we carry out  simulations for    the    parameter estimation problem.   We aim to   estimate  the    unknown  $m $-dimensional parameter $x^*$ based on the  gathered scalar  measurements $ \{d_k\}_{k\geq 1}$ given by
$  d_{k}=u_{k}^Tx^* +\nu_{k},$  where   $u_{k} \in  \mathbb{R}^m$ denotes the regression vector and $\nu_{ k}  \in \mathbb{R}$ denotes   the local observation noise.
Assume that $\{u_k\}$ and $\{\nu_k\}$ are mutually independent   i.i.d.  Gaussian sequences with  distributions $N(\mathbf{0}, R_u)$ and    $N(0,\sigma_{\nu}^2)$, respectively. Suppose   the   covariance matrix $  R_u$    is  positive definite.  Then we might  model the parameter  estimation problem  as   the following stochastic optimization problem:
 \begin{equation}\label{filter1}
  \min_{x\in \mathbb{R}^m } ~  f(x)\triangleq  \mathbb{E}   \big[\| d_{ k}-u_{k}^Tx\|^2\big].
  \end{equation}
Thus,   $f(x)=(x-x^*)^T R_u(x-x^*)+\sigma_{\nu}^2 $ and $\nabla f(x)= R_u(x-x^*).$
 Because the Hessian matrix $R_u$ of the objective function  $f(\cdot)$ is positive  definite, $x^*$ is the unique optimal solution to  \eqref{filter1}.
Suppose  we are able to   observe the  regressor  $u_{k}$ and  the measurement $d_{k}$,  then the noisy observation  of the   gradient  $\nabla f(x) $ might  be  constructed as $ u_ku_k^T x-d_k u_k$. Set the dimension of $x^*$ as $m=5.$ We run Algorithm \ref{Alg_1} with  $\alpha={2\over L+\eta}$, Algorithm \ref{Alg_2} with $\alpha={1 \over L}$  and  $\beta = {\sqrt{\kappa} -1 \over\sqrt{\kappa}+1},$ and  Algorithm \ref{Alg_3}    with   $\alpha = { 4 \over (\sqrt{\eta}+\sqrt{L})^2}$ and $\beta =  \left(  { \sqrt{\kappa} -1 \over \sqrt{\kappa} +1}\right)^2,$    where the  initial values  $x_{0}=y_0=0$ and the batch-size  $N_k=\lceil \rho^{-k}\rceil$  with $\rho={\kappa^2 \over (\kappa+1)^2}$.

 {\bf Convergence rate, iteration and oracle complexity. } We run
Algorithms  \ref{Alg_1},    \ref{Alg_2},   \ref{Alg_3},  and
the standard SGD algorithm $ x_{k+1}=x_k - \alpha_k  \nabla f(x_k,\xi_{j,k})$
with  $\alpha_k=R_u^{-1}/k$   setting to be the optimal tuning steplength, and
terminate the schemes when   $ { \mathbb{E}[\| x_k-x^*\|_2] \over
\|x^*\|_2}\leq 10^{-3}. $   We then examine their  empirical rate  of
convergence,   iteration and oracle complexity. Here the empirical mean is
calculated by averaging across  100 trajectories. The convergence rate of   the
relative error $ { \mathbb{E}[\| x_k-x^*\|_2] \over \|x^*\|_2} $   is   shown
in Figure \ref{GOne},  which demonstrates   that the iterates   generated by
Algorithms  \ref{Alg_1}-\ref{Alg_3}  converge  in mean to the  optimal solution
at a linear rate. It is noticed  that the accelerated  scheme  (Algorithm
\ref{Alg_2}) has the fastest convergence rate,  while the heavy ball method
(Algorithm  \ref{Alg_3}) displays the   most unsteadiness where the  error goes
up at the first few iterations.  The empirical relation between the accuracy
$\epsilon$ and $K(\epsilon)$ is shown in Figure \ref{Giter}, where
$K(\epsilon)$ denotes  the   number of iterations  required   to make  $  {
\mathbb{E}[\| x_k-x^*\|_2] \over \|x^*\|_2} <\epsilon.$ It is seen that  the
standard SGD algorithm requires much more iterations than the proposed  variance reduced schemes  for obtaining  an approximate solution with the same accuracy.  The empirical relation between $\epsilon$ and $O(\epsilon)$ is shown
in Figure \ref{GTwo}, where  $O(\epsilon)$ denotes  the   number of    sampled
gradients required   to make  $  { \mathbb{E}[\| x_k-x^*\|_2] \over \|x^*\|_2}
<\epsilon.$ We can see that for obtaining an estimate with the same accuracy,
the accelerated  scheme  (Algorithm \ref{Alg_3}) requires the smallest number of sampled gradients,
and  the variable sample-size SGD method (Algorithm \ref{Alg_1})   requires
more  sampled gradients than  the  standard SGD.

 \begin{figure}[htbp]
 \centering
 \begin{minipage}{0.33\linewidth}
       \centering
  \includegraphics[width=2.2in]{rate}
  \caption{Linear Rate  }
   \label{GOne}
    \end{minipage
 \begin{minipage}{0.03\linewidth}
\end{minipage
  \begin{minipage}{0.33\linewidth}
       \centering
 \includegraphics[width=2.2in]{iter_comp}
 \caption{Iter. Complexity}     \label{Giter}
    \end{minipage}
  \begin{minipage}{0.33\linewidth}
       \centering
 \includegraphics[width=2.2in]{oracle_comp}
 \caption{Oracle Complexity}     \label{GTwo}
    \end{minipage
\end{figure}

  {\bf Display the limiting  distributions.}   We run     Algorithms  \ref{Alg_1},   \ref{Alg_2}, and   \ref{Alg_3}  by  $1000$  independent sample paths and terminate the schemes at $k=50$.  The  empirical largest eigenvalue of the   covariance matrix $\mathbf{S}_k$ (estimated by \eqref{est-covariance})  at $k=50$
are  $1.648\times 10^{-5} ,  1.649\times 10^{-6} , 1.825\times 10^{-6} $, respectively.     This might  imply that the   accelerated   scheme (Algorithm  \ref{Alg_2}) has the best  performance measured by the covariance of the stationary distribution.  Since the unknown
parameter $x^*$ is multi-dimensional, we merely show  the limiting  distribution of one component of the   rescaled error    $\rho^{-k/2}(x_{ k}-x^*).$ The histograms  of    $\rho^{-k/2}(x_{ k}^5-x^*)$  at    $k=50 $ are shown in   Figure  \ref{GThree}   along with the fitted normal distribution (the red curve), where $x_k^5 $ denotes the fifth component of $x_k.$    It is also seen that the rescaled error approximately satisfy the normal distribution and among them the  accelerated   scheme  has the smallest variance. In addition,    the histograms of the rescaled  suboptimality gap $\rho^{-k}(f(x_k)-f(x^*))$ at $k=50$  are shown in Figure \ref{Ggap}.  We can further conclude that the  empirical sub-optimality gap of the   accelerated   scheme  is the best among  the proposed variable sample-size schemes. We further run  the SGD method  $ x_{k+1}=x_k - R_u^{-1} \nabla f(x_k,\xi_{j,k})/k$.  The histogram  along with the fitted normal distribution     of   $\sqrt{N}(x_N^5-x^*)$     is   displayed in Figure  \ref{GFour},  while the  empirical sub-optimality gap $N(f(x_N)-f(x^*))$  is shown in Figure  \ref{Gsub},  where $N\triangleq \sum_{k=1}^{50} N_k $ denotes the number of sampled gradients utilized by    Algorithms  \ref{Alg_1}-\ref{Alg_3}.

\begin{figure}[htbp]
       \centering
  \includegraphics[width=6.3in]{hist}
 \caption{ Histograms of   $\rho^{-k/2}(x_k^5-x^*) $ at $k=50$ along fitted normal distributions }      \label{GThree}
\vspace{-0.2in}
\end{figure}

\begin{figure}[htbp]
\centering
 \includegraphics[width=6.3in]{sub_opt_vss}
 \caption{ Histograms of $\rho^{-k}(f(x_k)-f(x^*))$ at $k=50$ }     \label{Ggap}
\end{figure}

 \begin{figure}[htbp]
  \begin{minipage}{0.66\linewidth}
       \centering
   \includegraphics[width=3in]{hist_sgd}
 \caption{ Hist.    along fitted normal for $\sqrt{N}(x_N^5-x^*)$ }     \label{GFour}
    \end{minipage
 \begin{minipage}{0.02\linewidth}
\end{minipage
\begin{minipage}{0.3\linewidth}
       \centering
   \includegraphics[height=1.8in,width=2.1in]{sub_opt_sgd}
 \caption{ Hist.     of  $N(f(x_N)-f(x^*))$ }     \label{Gsub}
    \end{minipage
\vspace{-0.2in}
\end{figure}

  {\bf Coverage probability of the constructed  confidence region.} In one single  replication,  we run     Algorithms  \ref{Alg_1},    \ref{Alg_2}  and   \ref{Alg_3}  by  $n$  independent sample paths and terminate them when the total  number of sampled gradients used reaches $N_{\max}.$
Then we can  construct the $95\%$ confidence region  defined by  \eqref{def-conf-region} and check whether the true parameter lies in the  constructed  confidence region.  We estimate the  coverage probability (i.e., the proportion of  replications  that the confidence region  contains the true value)  by   1000 replications.  The estimated coverage probability for  Algorithms \ref{Alg_1}-\ref{Alg_3}  with different selection of  sample paths $n$ and  simulation budget $N_{\max} $ is shown in Table \ref{Tab1}. It can be  seen that the coverage probabilities are getting closer to the nominal level $95\%$ when the sample paths $n$ grows larger.  It seems  that the simulation budget $N_{\max}$ does not  significantly impact  the coverage probability, however,   larger  $N_{\max}$ does  narrow  the confidence region.

  \begin{table} [htbp]
 \centering
 \begin{tabular}{|c|c|c|c|c|c|c|c|c|}
 \hline   \multicolumn{2}{|c|}{ } & VSS-SGD &VSS-ACC& VSS-HB  \\  \hline
{\multirow{3}{*}{$N_{\max}=10^3$}}  &  $n=6$ &  $    0.671\pm  0.221$ &$ 0.652	\pm 0.2259$  &$  0.665	\pm
    0.223$ \\   \cline{2-5}
 &  $n=8 $     & $  0.904	\pm  0.0869 $  & $  0.872\pm0.1117$ &$   0.893	\pm   0.0995$    \\ \cline{2-5}
 &  $n=10 $     & $ 0.946	\pm   0.0511$  & $  0.923\pm 0.0591$ &$ 0.951	\pm   0.0466$     \\ \cline{2-5}
 &  $ n=15$    &  $    0.973	\pm   0.0263$ &$  0.966\pm     0.0329$  &$ 0.961\pm   0.0375$    \\ \hline
{\multirow{3}{*}{$N_{\max}=10^4$}}  &  $n=6$ &  $   0.664\pm   0.2233$ &$ 0.62	\pm    0.2358$  &$   0.646 \pm     0.2289 $ \\   \cline{2-5}
 &  $n=8 $     & $  0.872	\pm  0.117 $  & $  0.877\pm  0.108$ &$   0.896	\pm   0.0933$    \\ \cline{2-5}
 &  $n=10 $     & $ 0.929\pm    0.066$  & $      0.944 \pm 0.0592$ &$ 0.949	\pm   0.0484$     \\ \cline{2-5}
 &  $ n=15$    &  $       0.965	\pm     0.0338$ &$      0.946 \pm       0.0511$  &$ 0.961\pm   0.0375$    \\ \hline
{\multirow{3}{*}{$N_{\max}=10^5$}}  &  $n=6$ &  $ 0.651	\pm 0.2274$ &$ 0.656	\pm 0.2259$  &$ 0.664	\pm  0.2233$ \\   \cline{2-5}
 &  $n=8 $     & $ 0.87	\pm 0.1132$  & $ 0.871\pm0.1125$ &$   0.896	\pm   0.0933$    \\ \cline{2-5}
 &  $n=10 $     & $ 0.912	\pm 0.083$  & $ 0.937\pm 0.0591$ &$ 0.919	\pm 0.0745$    \\ \cline{2-5}
 &  $ n=15$    &  $ 0.952	\pm 0.0457$ &$ 0.964	\pm 0.0347$  &$ 0.956	\pm 0.0421$    \\ \hline
 \end{tabular}
  \caption{The estimated coverage probability  of  Algorithms \ref{Alg_1}-\ref{Alg_3}. The ideal coverage probability is 0.95. }  \label{Tab1}
\vspace{-0.2in}
\end{table}

  {\bf Polynomially increasing batch-sizes. } We run Algorithms  \ref{Alg_1}-\ref{Alg_3} with the batch-size  increasing  at a polynomial rate $N_k=\lceil k^{v}\rceil ,v=2 $.  The convergence rate of   the relative error $ { \mathbb{E}[\| x_k-x^*\|_2] \over \|x^*\|_2} $   is   shown in Figure \ref{GPoly1}, while the    histograms  of    $k^{ v/2}(x_{ k}^5-x^*)$  at    $k=100 $  along with the fitted normal distribution are shown in   Figure  \ref{GPoly2}.
It can be seen that the  accelerated  scheme  has the best performance because it displays  fastest convergence rate and the rescaled error has smallest variance.

 \begin{figure}[htbp]
     \begin{minipage}{0.34\linewidth}
       \centering
 \includegraphics[width=2.2in]{poly_rate}
 \caption{Poly. Batch-size}     \label{GPoly1}
    \end{minipage
 \begin{minipage}{0.01\linewidth}
\end{minipage
  \begin{minipage}{0.64\linewidth}
       \centering
   \includegraphics[width=4.4in]{poly_hist}
 \caption{ Histograms        of  $k^{ v/2}(x_{ k}^5-x^*)$  at $k=100$ }     \label{GPoly2}
    \end{minipage
\vspace{-0.2in}
\end{figure}

\section{Conclusions} In this work, we considered the
strongly convex stochastic optimization problem and
proposed three classes of  variance reduced stochastic
gradient descent algorithms (unaccelerated, accelerated, and
heavy ball methods), where the  unavailable exact gradient
is approximated by an increasing batch of sampled gradients.
We then establish convergent rate and complexity guarantees.
Further, we establish amongst the first formal central limit
theorems for  all the three schemes when the batch-size
increased at either a geometric or a polynomial  rate.  The
covariance matrix specifies how problem   structure
(including the strong convexity parameter,  the Lipschitz
constant, and the Hessian matrix) and distribution  of
gradient   noise influences the algorithm  performance.  In
addition,  we provide an avenue to construct  the
confidence region  of the  optimal solution based on the
central limit theorems. Finally, we apply  the proposed
schemes to the stochastic  parameter estimation problem to
validate our theoretical  findings.

\def\cprime{$'$}