# Online Statistical Inference for Stochastic Optimization via Kiefer-
Wolfowitz Methods

Xi Chen111Stern School of Business, New York University, Email:
<EMAIL_ADDRESS>Zehua Lai222Committee on Computational and Applied
Mathematics, University of Chicago, Email<EMAIL_ADDRESS>He Li333Stern
School of Business, New York Universitwy, Email<EMAIL_ADDRESS>Yichen
Zhang444Krannert School of Management, Purdue University, Email:
<EMAIL_ADDRESS>
###### Abstract

This paper investigates the problem of online statistical inference of model
parameters in stochastic optimization problems via the Kiefer-Wolfowitz
algorithm with random search directions. We first present the asymptotic
distribution for the Polyak-Ruppert-averaging type Kiefer-Wolfowitz (AKW)
estimators, whose asymptotic covariance matrices depend on the distribution of
search directions and the function-value query complexity. The distributional
result reflects the trade-off between statistical efficiency and function
query complexity. We further analyze the choice of random search directions to
minimize certain summary statistics of the asymptotic covariance matrix. Based
on the asymptotic distribution, we conduct online statistical inference by
providing two construction procedures of valid confidence intervals.

Keywords: Asymptotic normality, Kiefer-Wolfowitz stochastic approximation,
online inference, stochastic optimization

##  1 Introduction

Stochastic optimization algorithms, introduced by Robbins and Monro (1951);
Kiefer and Wolfowitz (1952), have been widely used in statistical estimation,
especially for large-scale datasets and online learning where the sample
arrives sequentially (e.g., web search queries, transactional data). The
Robbins-Monro algorithm (Robbins and Monro, 1951), often known as the
stochastic gradient descent, is perhaps the most popular algorithm in
stochastic optimization and has found a wide range of applications in
statistics and machine learning. Nevertheless, in many modern applications,
the gradient information is not available. For example, the objective function
may be embedded in a black box and the user can only access the noisy
objective value for a given input. In such cases, the Kiefer-Wolfowitz
algorithm (Kiefer and Wolfowitz, 1952) becomes a natural choice as it is
completely free of gradient computation. Despite being equipped with an
evident computational advantage to avoid gradient measurements, the Kiefer-
Wolfowitz algorithm has been historically out of practice as compared to the
Robbins-Monro counterpart. Nonetheless, heralded by the big data era, there
has been a restoration of the interest of gradient-free optimization in a wide
range of applications in recent years (Conn et al., 2009; Nesterov and
Spokoiny, 2017). We briefly highlight a few of them to motivate our paper.

  * •

In some bandit problems, one may only have black-box access to individual
objective values but not to their gradients (Flaxman et al., 2005; Shamir,
2017). Other examples include graphical models and variational inference
problems, where the objective is defined variationally (Wainwright and Jordan,
2008), and the explicit differentiation can be difficult.

  * •

In some scenarios, the computation of gradient information is possible but
very expensive. For example, in the online sensor selection problem (Joshi and
Boyd, 2008), evaluating the stochastic gradient requires the inverse of
matrices, which generates $\mathcal{O}(d^{3})$ computation cost per iteration,
where $d$ is the number of sensors in the network. In addition, the storage
for gradient calculation also requires an $O(d^{3})$ memory, which could be
practically infeasible.

  * •

In some statistical problems such as quantile regression and its variants
(Koenker, 2005), the objective function is not differentiable. Extending the
gradient definition to nonsmooth functions is generally nontrivial, and
techniques of defining sets of local differential characteristics suffer from
the incompleteness of chain rule in complex problems (Nesterov, 2005).

This paper aims to study the asymptotic properties of the Kiefer-Wolfowitz
stochastic optimization and conduct online statistical inference. In
particular, we consider the problem,

$\displaystyle\boldsymbol{\theta}^{\star}=\mathrm{argmin}\leavevmode\nobreak\
F(\boldsymbol{\theta}),\quad\text{where}\;F(\boldsymbol{\theta}):=\mathbb{E}_{\mathcal{P}_{\boldsymbol{\zeta}}}\left[f(\boldsymbol{\theta};\boldsymbol{\zeta})\right]=\int
f(\boldsymbol{\theta};\boldsymbol{\zeta})\text{d}\mathcal{P}_{\boldsymbol{\zeta}},$
(1)

where $f(\boldsymbol{\theta};\boldsymbol{\zeta})$ is a convex _individual loss
function_ for a data point $\boldsymbol{\zeta}$, $F(\boldsymbol{\theta})$ is
the _population loss_ function, and $\boldsymbol{\theta}^{\star}$ is the true
underlying parameter of a fixed dimension $d$. Let $\boldsymbol{\theta}_{0}$
denote any given initial point. Given a sequentially arriving online sample
$\\{\boldsymbol{\zeta}_{n}\\}$, the Robbins and Monro (1951) algorithm (RM),
also known as the stochastic gradient descent (SGD), iteratively updates,

  * (RM)

$\boldsymbol{\theta}_{n}^{\tt{(RM)}}=\boldsymbol{\theta}_{n-1}^{\tt{(RM)}}-\eta_{n}g(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$,
(2)

where $\\{\eta_{n}\\}$ is a positive non-increasing step-size sequence, and
$g(\boldsymbol{\theta};\boldsymbol{\zeta})$ denotes the stochastic gradient,
i.e., $g(\boldsymbol{\theta};\boldsymbol{\zeta})=\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})$. In the scenarios that direct
gradient measurements are inaccessible to practitioners, the Kiefer and
Wolfowitz (1952) algorithm (KW) becomes the natural choice, as

  * (KW)

$\boldsymbol{\theta}_{n}^{\tt{(KW)}}=\boldsymbol{\theta}_{n-1}^{\tt{(KW)}}-\eta_{n}\widehat{g}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$,
(3)

where $\widehat{g}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$ is an
estimator of $g(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$. Under the
univariate framework ($d=1$), Kiefer and Wolfowitz (1952) considered the
finite-difference approximation

$\displaystyle\widehat{g}(\theta_{n-1};\zeta_{n})=\frac{f(\theta_{n-1}+h_{n};\zeta_{n})-f(\theta_{n-1};\zeta_{n})}{h_{n}},$
(4)

where $h_{n}$ is be a positive deterministic sequence that goes to zero. Blum
(1954) later extended the algorithm to the multivariate case and proved its
almost sure convergence. This pioneering work extended in various directions
of statistics and control theory (see, e.g., Fabian (1967, 1978); Hall and
Heyde (1980); Ruppert (1982); Chen (1988); Polyak and Tsybakov (1990); Spall
(1992); Chen et al. (1999); Spall (2000); Hall and Molchanov (2003); Dippon
(2003); Mokkadem and Pelletier (2007); Broadie et al. (2011)). In the
optimization literature, the Kiefer-Wolfowitz (KW) algorithm is often referred
to as the gradient-free stochastic optimization, or zeroth-order SGD (Agarwal
et al., 2010; Agarwal et al., 2011; Jamieson et al., 2012; Ghadimi and Lan,
2013; Duchi et al., 2015; Shamir, 2017; Nesterov and Spokoiny, 2017; Wang et
al., 2018, among others).

For the (RM) algorithm in ((RM)), Ruppert (1988) and Polyak and Juditsky
(1992) characterize the limiting distribution and statistical efficiency of
the _averaged iterate_
$\overline{\boldsymbol{\theta}}_{n}^{{\tt{(RM)}}}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{\theta}_{i}^{{\tt{(RM)}}}$
by

$\displaystyle\sqrt{n}\left(\overline{\boldsymbol{\theta}}_{n}^{{\tt{(RM)}}}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}SH^{-1}\right),$
(5)

where $H=\nabla^{2}F(\boldsymbol{\theta}^{\star})$ is the Hessian matrix of
$F(\boldsymbol{\theta})$ at $\boldsymbol{\theta}=\boldsymbol{\theta}^{\star}$,
and $S=\mathbb{E}[\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}]$ is the Gram matrix
of the stochastic gradient. Under a well-specified model, this asymptotic
covariance matrix matches the inverse Fisher information and the averaged (RM)
estimator is asymptotically efficient. Based on the limiting distribution
result (5), there are many recent research efforts devoted to statistical
inference for (RM). A brief survey is conducted at the end of the
introduction.

For the (KW) scheme, we can similarly construct the averaged Kiefer-Wolfowitz
(AKW) estimator

  * (AKW)

$\overline{\boldsymbol{\theta}}_{n}^{\tt{(KW)}}=\frac{1}{n}\sum\limits_{i=1}^{n}\boldsymbol{\theta}_{i}^{\tt{(KW)}}$.
(6)

As compared to well-established asymptotic properties of (RM), study of the
asymptotics of (AKW) is limited, particularly with a random sampling direction
in multivariate (KW). In this paper, we study the (KW) algorithm ((KW)) with
random search directions
$\\{\boldsymbol{v}_{i}\\}_{i=1}^{n}\overset{i.i.d.}{\sim}\mathcal{P}_{\boldsymbol{v}}$,
i.e., at each iteration $i=1,2,\dots,n$, a random direction
$\boldsymbol{v}_{i}$ is sampled independently from
$\mathcal{P}_{\boldsymbol{v}}$, and the (KW) gradient

$\displaystyle\widehat{g}_{h_{n},\boldsymbol{v}_{n}}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})=\frac{f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n};\boldsymbol{\zeta}_{n})-f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})}{h_{n}}\boldsymbol{v}_{n}.$
(7)

Compared to the (RM) scheme, (KW) introduces additional randomness into the
stochastic gradient estimator through $\\{\boldsymbol{v}_{n}\\}$. Indeed, as
one can see from our main result in Theorem 3.3, (AKW) is no longer
statistically efficient and its asymptotic covariance structure depends on the
distribution $\mathcal{P}_{\boldsymbol{v}}$. It opens the room for the
investigation on the impact of $\mathcal{P}_{\boldsymbol{v}}$ (see Section 3.1
for details). We further extend the estimator to utilize multiple function-
value queries per step and establish an online statistical inference
framework. We summarize our main results and contributions as follows,

  * •

First, we quantify the asymptotic covariance structure of (AKW) in Theorem
3.3. Since the asymptotic distribution depends on the choice of the direction
variable $\boldsymbol{v}$, we provide an introductory analysis on the
asymptotic performance for different choices of random directions for
constructing (AKW) estimators (see Section 3.1).

  * •

The efficiency loss of (AKW) is due to the information constraint as one
evaluates only _two_ function values at each iteration. We analyze the (AKW)
estimators in which multiple function queries can be assessed at each
iteration, and show that the asymptotic covariance matrix decreases as the
number of function queries $m+1$ increases (see Section 3.2). Moreover, (AKW)
achieves asymptotic statistical efficiency as $m\rightarrow\infty$. We further
show that when $\boldsymbol{v}$ is sampled without replacement from
$\mathcal{P}_{\boldsymbol{v}}$ with a discrete uniform distribution of any
orthonormal basis, (AKW) achieves asymptotic statistical efficiency with $d+1$
function queries per iteration.

  * •

Based on the asymptotic distribution, we propose two online statistical
inference procedures. The first one is using a plug-in estimator of the
asymptotic covariance matrix, which separately estimates the Hessian matrix
and Gram matrix of the (KW) gradients (with additional function-value queries,
see Theorem 4.3). The second procedure is to characterize the distribution of
intermediate (KW) iterates as a stochastic process and construct an
asymptotically pivotal statistic by normalizing the (AKW) estimator, without
directly estimating the covariance matrix. This inference procedure is
inspired by the “random scaling” method proposed in Lee et al. (2022) that
considers the online inference for the (RM) scheme. These two procedures have
their advantages and disadvantages: the plug-in approach leads to better
empirical performance but requires additional function-value queries to
estimate the Hessian matrix, while the other one is more efficient in both
computation and storage, though its finite-sample performance is inferior in
practice when the dimension is large. A practitioner may choose the approach
suitable to her computational resources and requirement of the inference
accuracy.

Lastly, we provide a brief literature survey on the recent works for
statistical inference for the (RM)-type SGD algorithms. Chen et al. (2020)
developed a batch-means estimator of the limiting covariance matrix
$H^{-1}SH^{-1}$ in (5), which only uses the stochastic gradient information
(i.e., without estimating any Hessian matrices). Zhu et al. (2021) further
extended the batch-means method in Chen et al. (2020) to a fully online
covariance estimator. Lee et al. (2022) extended the results in Polyak and
Juditsky (1992) to a functional central limit theorem and utilize it to
propose a novel online inference procedure that allows for efficient
implementation. Fang et al. (2017) presented a perturbation-based resampling
procedure for inference. Su and Zhu (2018) proposed a tree-structured
inference scheme, which splits the SGD into several threads to construct
confidence intervals. Liang and Su (2019) introduced a moment-adjusted method
and its corresponding inference procedure. Toulis and Airoldi (2017)
considered the implicit SGD, and investigate the statistical inference problem
under the variant. Duchi and Ruan (2021) studied the stochastic optimization
problem with constraints and investigate its optimality properties. Chao and
Cheng (2019) proposed a class of generalized regularized dual averaging (RDA)
algorithms and make uncertainty quantification possible for online
$\ell_{1}$-penalized problems. Shi et al. (2021) developed an online
estimation procedure for high-dimensional statistical inference. Chen et al.
(2021) studied statistical inference of online decision-making problems via
SGD in a contextual bandit setting.

###  1.1 Notations and organization of the paper

We write vectors in boldface letters (e.g., $\boldsymbol{\theta}$ and
$\boldsymbol{v}$) and scalers in lightface letters (e.g., $\eta$). For any
positive integer $n$, we use $[n]$ as a shorthand for the discrete set
$\\{1,2,\cdots,n\\}$. Let $\\{\boldsymbol{e}_{k}\\}_{k=1}^{d}$ be the standard
basis in $\mathbb{R}^{d}$ with the $k$-th coordinate as $1$ and the other
coordinates as $0$. Denote $I_{d}$ as the identity matrix in
$\mathbb{R}^{d\times d}$. Let $\|\cdot\|$ denote the standard Euclidean norm
for vectors and the spectral norm for matrices. We use $A_{k\ell}$ and
$A_{n,k\ell}$ to denote the $(k,\ell)$-th element of matrices
$A,A_{n}\in\mathbb{R}^{d\times d}$, respectively, for all $k,\ell\in[d]$.
Furthermore, we denote by $\mathrm{diag}(\boldsymbol{v})$ a matrix in
$\mathbb{R}^{d\times d}$ whose main diagonal is the same as the vector
$\boldsymbol{v}$ and off-diagonal elements are zero, for some vector
$\boldsymbol{v}\in\mathbb{R}^{d}$. With a slight abuse of notation, for a
matrix $M\in\mathbb{R}^{d\times d}$, we also let $\mathrm{diag}(M)$ denote a
$\mathbb{R}^{d\times d}$ diagonal matrix with same diagonal elements as matrix
$M$. We use the standard Loewner order notation $A\succeq 0$ if a matrix $A$
is positive semi-definite. We use $\boldsymbol{\theta}^{\tt{(RM)}}$ and
$\boldsymbol{\theta}^{\tt{(KW)}}$ to denote the iterates generated by the (RM)
scheme and the (KW) scheme, respectively. We use
$\widehat{\boldsymbol{\theta}}^{\mathtt{(ERM)}}$ for the offline empirical
risk minimizer, i.e.,
$\widehat{\boldsymbol{\theta}}^{\mathtt{(ERM)}}=\mathrm{argmin}_{\boldsymbol{\theta}}\frac{1}{n}\sum_{i=1}^{n}f(\boldsymbol{\theta};\boldsymbol{\zeta}_{i})$.
As we focus on the (KW) scheme in this paper, we sometimes omit the
superscript ${\tt{(KW)}}$ in the estimator to make room for the other
notations. In derivations of the (KW) estimator, we denote the finite
difference of $f(\cdot)$ as,

$\displaystyle\Delta_{h,\boldsymbol{v}}f(\boldsymbol{\theta};\boldsymbol{\zeta})=f(\boldsymbol{\theta}+h\boldsymbol{v};\boldsymbol{\zeta})-f(\boldsymbol{\theta};\boldsymbol{\zeta}),$
(8)

for some spacing parameter $h\in\mathbb{R}_{+}$ and search vector
$\boldsymbol{v}\in\mathbb{R}^{d}$. We use $\mathbb{E}_{n}$ to denote the
conditional expectation with respect to the natural filtration, i.e.,

$\displaystyle\mathbb{E}_{n}[\boldsymbol{\theta}_{n+1}]:=\mathbb{E}[\boldsymbol{\theta}_{n+1}|\mathcal{F}_{n}],\;\;\mathcal{F}_{n}:=\sigma\\{\boldsymbol{\theta}_{k},\boldsymbol{\zeta}_{k}|k\leq
n\\}.$

We use the $\mathcal{O}(\cdot)$ notation to hide universal constants
independent of the sample size $n$.

The remainder of the paper is organized as follows. In Section 2, we describe
the Kiefer-Wolfowitz algorithm with random search directions along with three
illustrative examples of the classical regression problems. We also provide a
technical lemma to characterize the limiting behavior of the (KW) gradient,
which leads to the distributional constraint of the random direction vector.
In Section 3, we first introduce the technical assumptions before we present
the finite-sample rate of convergence of the (KW) estimator. We further
provide the asymptotic distribution of the (AKW) estimator, accompanied by
discussions on the statistical (in)efficiency. We highlight a comparison of
the choices of the direction distributions in Section 3.1, and further extend
the theoretical analysis to multi-query settings of the (KW) algorithm in
Section 3.2. Based on the established asymptotic distribution results, we
propose two types of online statistical inference procedures in Section 4. A
functional extension of the distributional analysis of (KW) as a stochastic
process is also provided. Numerical experiments in Section 5 lend empirical
support to our theory. Further discussions are provided in Section 6 and all
proofs are relegated to the Appendix.

##  2 Kiefer-Wolfowitz Algorithm

In this section, we introduce the general form of the Kiefer-Wolfowitz (KW)
gradient estimator and the corresponding iterative algorithm
$\boldsymbol{\theta}_{n}=\boldsymbol{\theta}_{n-1}-\eta_{n}\widehat{g}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$.
In the seminal work by Blum (1954), the (KW) gradient estimator
$\widehat{g}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$ is constructed
by approximating the stochastic gradient
$g(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$ using the canonical
basis of $\mathbb{R}^{d}$,
$\\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\dots,\boldsymbol{e}_{d}\\}$, as
search directions. In particular, given any
$\boldsymbol{\theta}\in\mathbb{R}^{d}$ and
$\boldsymbol{\zeta}\sim\mathcal{P}_{\boldsymbol{\zeta}}$, the $k$-th
coordinate of the (KW) gradient estimator

$\displaystyle\big{(}\widehat{g}_{h,\boldsymbol{e}}(\boldsymbol{\theta};\boldsymbol{\zeta})\big{)}_{k}=\frac{f(\boldsymbol{\theta}+h\boldsymbol{e}_{k};\boldsymbol{\zeta})-f(\boldsymbol{\theta};\boldsymbol{\zeta})}{h},\qquad\text{for}\quad
k=1,2,\ldots,d,$ (9)

where $h$ is a spacing parameter for approximation. At each iteration, (9)
queries $d+1$ function values from $d$ fixed directions
$\\{\boldsymbol{e}_{k}\\}_{k=1}^{d}$. To reduce the query complexity, a random
difference becomes a natural choice. Koronacki (1975) introduced a random
version of the (KW) algorithm using a sequence of random unit vectors that are
independent and uniformly distributed on the unit sphere or unit cube. Spall
(1992) also provided a random direction version of the (KW) algorithm, named
as the simultaneous perturbation stochastic approximation (SPSA) algorithm and
later extended to several variants Chen et al. (1999); Spall (2000); He et al.
(2003). These random direction methods can reduce the bias in gradient
estimates as compared to their non-random counterparts. In the following, we
write the (KW) algorithm with general random search directions, as in (7),

$\displaystyle\boldsymbol{\theta}_{n}$
$\displaystyle=\boldsymbol{\theta}_{n-1}-\eta_{n}\widehat{g}_{h_{n},\boldsymbol{v}_{n}}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n}),$
$\displaystyle\text{where
}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta}):=\frac{1}{h}\Delta_{h,\boldsymbol{v}}f(\boldsymbol{\theta};\boldsymbol{\zeta})\boldsymbol{v}=\frac{f(\boldsymbol{\theta}+h\boldsymbol{v};\boldsymbol{\zeta})-f(\boldsymbol{\theta};\boldsymbol{\zeta})}{h}\boldsymbol{v}.$
(10)

Here $\\{\boldsymbol{v}_{n}\\}$ is sampled from an underlying distribution
$\mathcal{P}_{\boldsymbol{v}}$ satisfying certain conditions (see Assumption 4
in Section 3). At each iteration $n$, the algorithm samples a direction vector
$\boldsymbol{v}_{n}$ independently from $P_{\boldsymbol{v}}$, and makes two
solitary function-value queries,
$f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$ and
$f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n};\boldsymbol{\zeta}_{n})$.
We refer to the (KW) gradient estimator
$\widehat{g}_{h_{n},\boldsymbol{v}_{n}}(\boldsymbol{\theta}_{n-1},\boldsymbol{\zeta}_{n})$
in (10) as a _two-query_ finite-difference approximation of the stochastic
gradient. If one is allowed to make additional function-value queries, an
averaging of the function values from multiple directions generates a _multi-
query_ stochastic gradient estimator with reduced variance. In particular, at
each iteration $n$, the practitioner makes $m+1$ queries
$\\{f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n}),f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n}^{(j)};\boldsymbol{\zeta}_{n})\\}_{1\leq
j\leq m}$ via $m$ random directions
$\big{\\{}\boldsymbol{v}_{n}^{(j)}\big{\\}}$ sampled from
$\mathcal{P}_{\boldsymbol{v}}$. If $\mathcal{P}_{\boldsymbol{v}}$ is a finite
distribution, practitioners may choose to sample _with_ or _without
replacement_. In summary, an ($m+1$)_-query_ (KW) algorithm constructs a
stochastic gradient estimator

$\overline{g}_{n}^{(m)}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})=\frac{1}{m}\sum_{j=1}^{m}\widehat{g}_{h_{n},\boldsymbol{v}_{n}^{(j)}}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})=\frac{1}{mh_{n}}\sum_{j=1}^{m}\Delta_{h_{n},\boldsymbol{v}_{n}^{(j)}}f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})\boldsymbol{v}_{n}^{(j)},$
(11)

at each iteration $n$, and updates
$\boldsymbol{\theta}_{n}=\boldsymbol{\theta}_{n-1}-\eta_{n}\overline{g}_{n}^{(m)}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})$.
Here we restrict the procedure to sampling from the same distribution
$\mathcal{P}_{\boldsymbol{v}}$ independently across different iterations. We
use $\boldsymbol{\theta}_{n}^{(m)}$ to denote the final (KW) estimator using
the above $(m+1)$-query finite-difference approximation.

We now provide some illustrative examples of the two-query (KW) estimator
$\widehat{g}_{h_{n},\boldsymbol{v}_{n}}$ in (10) used in popular statistical
models, and we will refer to these examples throughout the paper. A multi-
query extension of the examples can be constructed accordingly.

######  Example 2.1 (Linear Regression).

Consider a linear regression model
$y_{i}=\boldsymbol{x}_{i}^{\top}\boldsymbol{\theta}^{\star}+\epsilon_{i}$
where
$\\{\boldsymbol{\zeta}_{i}=(\boldsymbol{x}_{i},y_{i}),\leavevmode\nobreak\
i=1,2,\dots,n\\}$ is an i.i.d. sample of
$\boldsymbol{\zeta}=(\boldsymbol{x},y)$ and the noise
$\epsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$. We use a quadratic loss function
$f(\boldsymbol{\theta};\boldsymbol{\zeta})=(y-\boldsymbol{x}^{\top}\boldsymbol{\theta})^{2}$.
Therefore, the stochastic gradient $\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})=\left(\boldsymbol{x}^{\top}\boldsymbol{\theta}-y\right)\boldsymbol{x}$,
and the (KW) gradient estimator
$\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\\{\boldsymbol{x},y\\})$
in (10) becomes

$\displaystyle\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\\{\boldsymbol{x},y\\})$
$\displaystyle=\frac{1}{h}\left[\big{(}y-\boldsymbol{x}^{\top}(\boldsymbol{\theta}+h\boldsymbol{v})\big{)}^{2}-\big{(}y-\boldsymbol{x}^{\top}\boldsymbol{\theta}\big{)}^{2}\right]\boldsymbol{v}=\boldsymbol{v}\boldsymbol{v}^{\top}\big{(}\boldsymbol{x}^{\top}\boldsymbol{\theta}-y\big{)}\boldsymbol{x}+h(\boldsymbol{x}^{\top}\boldsymbol{v})^{2}\boldsymbol{v}.$

######  Example 2.2 (Logistic Regression).

Consider a logistic regression model with a binary response
$y_{i}\in\\{-1,1\\}$ generated by
$\Pr(y_{i}|\boldsymbol{x}_{i})=\left(1+\exp\left(-y_{i}\boldsymbol{x}_{i}^{\top}\boldsymbol{\theta}^{\star}\right)\right)^{-1}$
. The individual loss function
$f(\boldsymbol{\theta};\boldsymbol{\zeta})=\log\left(1+\exp(-y\boldsymbol{x}^{\top}\boldsymbol{\theta})\right)$.
The stochastic gradient $\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})=-y\boldsymbol{x}\left(1+\exp(y\boldsymbol{x}^{\top}\boldsymbol{\theta})\right)^{-1}$,
and the (KW) gradient estimator
$\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\\{\boldsymbol{x},y\\})$
in (10) becomes

$\displaystyle\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\\{\boldsymbol{x},y\\})$
$\displaystyle=\frac{\boldsymbol{v}}{h}\left[\log\big{(}1+\exp(-y\boldsymbol{x}^{\top}(\boldsymbol{\theta}+h\boldsymbol{v}))\big{)}-\log\big{(}1+\exp(-y\boldsymbol{x}^{\top}\boldsymbol{\theta})\big{)}\right]$
$\displaystyle=\frac{-y\boldsymbol{v}\boldsymbol{v}^{\top}\boldsymbol{x}}{1+\exp(y\boldsymbol{x}^{\top}\boldsymbol{\theta})}+\frac{y^{2}(\boldsymbol{x}^{\top}\boldsymbol{v})^{2}\exp(y\boldsymbol{x}^{\top}\boldsymbol{\theta})h\boldsymbol{v}}{2(1+\exp(y\boldsymbol{x}^{\top}\boldsymbol{\theta}))^{2}}+\mathcal{O}(h^{2}),\;\;\text{as
}h\rightarrow 0_{+},$

under some regularity conditions on $\boldsymbol{\theta}$ and the distribution
of $\boldsymbol{x}$.

######  Example 2.3 (Quantile Regression).

Consider a quantile regression model
$y_{i}=\boldsymbol{x}_{i}^{\top}\boldsymbol{\theta}^{\star}+\epsilon_{i}$
where
$\\{\boldsymbol{\zeta}_{i}=(\boldsymbol{x}_{i},y_{i}),\leavevmode\nobreak\
i=1,2,\dots,n\\}$ is an _i.i.d._ sample of
$\boldsymbol{\zeta}=(\boldsymbol{x},y)$ and the noise satisfies
$\Pr(\epsilon_{i}\leq 0|\boldsymbol{x}_{i})=\tau$. The individual loss
$f(\boldsymbol{\theta};\zeta)=\rho_{\tau}(y-\boldsymbol{x}^{\top}\boldsymbol{\theta})$,
where $\rho_{\tau}(z)=z(\tau-1_{\\{z<0\\}})$. Although $\rho_{\tau}$ is non-
differentiable, the (KW) gradient estimator $\widehat{g}_{h,\boldsymbol{v}}$
is well-defined and takes the following form,

$\displaystyle\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\\{\boldsymbol{x},y\\})=$
$\displaystyle\frac{\boldsymbol{v}}{h}\left[\rho_{\tau}\big{(}y-\boldsymbol{x}^{\top}(\boldsymbol{\theta}+h\boldsymbol{v})\big{)}-\rho_{\tau}\big{(}y-\boldsymbol{x}^{\top}\boldsymbol{\theta}\big{)}\right]$
$\displaystyle=$
$\displaystyle\boldsymbol{v}\boldsymbol{v}^{\top}\boldsymbol{x}\big{(}\tau-1_{\\{y-\boldsymbol{x}^{\top}\boldsymbol{\theta}<0\\}}\big{)},\quad\text{for
}0<h<\left|\frac{y-\boldsymbol{x}^{\top}\boldsymbol{\theta}}{\boldsymbol{x}^{\top}\boldsymbol{v}}\right|.$

We note that for the (RM) scheme with differentiable loss functions, the
stochastic gradient is an unbiased estimator of the population gradient under
very mild assumption, i.e.,
$\mathbb{E}_{\boldsymbol{\zeta}}g(\boldsymbol{\theta};\boldsymbol{\zeta})=\nabla
F(\boldsymbol{\theta})$. In contrast, the (KW) gradient estimator is no longer
an unbiased estimator of $\nabla F(\boldsymbol{\theta})$. In the following
lemma, we precisely quantifies the bias incurred by the (KW) gradient
estimator.

######  Lemma 2.4.

We assume that the population loss function $F(\cdot)$ is twice continuously
differentiable and $L_{f}$-smooth, i.e.,
$\nabla^{2}F(\boldsymbol{\theta})\preceq L_{f}I_{d}$ for any
$\boldsymbol{\theta}\in\mathbb{R}^{d}$. Given any fixed parameter
$\boldsymbol{\theta}\in\mathbb{R}^{d}$, suppose the random direction vector
$\boldsymbol{v}$ is independent from $\boldsymbol{\zeta}$, we have

$\displaystyle\left\|\mathbb{E}\;\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta})\right\|$
$\displaystyle\leq\left\|\mathbb{E}\big{(}\boldsymbol{v}\boldsymbol{v}^{\top}-I_{d}\big{)}\nabla
F(\boldsymbol{\theta})\right\|+\frac{h}{2}L_{f}\mathbb{E}\|\boldsymbol{v}\|^{3},$

where the expectation in
$\mathbb{E}\;\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})$
takes over both the randomness in $\boldsymbol{v}$ and $\boldsymbol{\zeta}$.

The proof of Lemma 2.4 is provided in Appendix A. To reduce the bias in the
(KW) gradient, Lemma 2.4 indicates that one should choose the random direction
$\boldsymbol{v}_{n}$ that satisfies the distributional constraint
$\mathbb{E}[\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}]=I_{d}$ (see
Assumption 4 in Section 3). We will further conduct a comprehensive analysis
in Section 3.1 on different choices of distributions
$\mathcal{P}_{\boldsymbol{v}}$ satisfying the condition
$\mathbb{E}[\boldsymbol{v}_{n}\boldsymbol{v}_{n}^{\top}]=I_{d}$. Despite the
existence of the bias, as the spacing parameter $h_{n}\rightarrow 0$, the bias
convergences to zero asymptotically.

##  3 Theoretical Results

We first introduce some regularity assumptions on the population loss
$F(\boldsymbol{\theta})$ and the individual loss
$f(\boldsymbol{\theta};\boldsymbol{\zeta})$. To simplify the illustration, we
present the following assumptions first, even though the main results remain
to hold under some weakened assumptions. The relaxation is discussed below,
and detailed theoretical derivation is relegated to Appendix A.2.

######  Assumption 1.

The population loss function $F(\boldsymbol{\theta})$ is twice continuously
differentiable. Moreover, there exists $L_{f}>\lambda>0$, such that, $\lambda
I_{d}\preceq\nabla^{2}F(\boldsymbol{\theta})\preceq L_{f}I_{d}$ for any
$\boldsymbol{\theta}\in\mathbb{R}^{d}$.

######  Assumption 2.

Assume $\mathbb{E}\left[\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta}_{n})\right]=\nabla
F(\boldsymbol{\theta})$ for any $\boldsymbol{\theta}\in\mathbb{R}^{d}$.
Moreover, for some $0<\delta\leq 2$, there exists $M>0$ such that
$\mathbb{E}\|\nabla f(\boldsymbol{\theta};\boldsymbol{\zeta}_{n})-\nabla
F(\boldsymbol{\theta})\|^{2+\delta}\leq
M\big{(}\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\|^{2+\delta}+1\big{)}.$

######  Assumption 3.

There are constants $L_{h},L_{p}>0$ such that for any
$\boldsymbol{\theta},\boldsymbol{\theta^{\prime}}\in\mathbb{R}^{d}$,

$\displaystyle\mathbb{E}\left\|\nabla^{2}f(\boldsymbol{\theta};\boldsymbol{\zeta}_{n})-\nabla^{2}f(\boldsymbol{\theta^{\prime}};\boldsymbol{\zeta}_{n})\right\|^{2}\leq
L_{h}\|\boldsymbol{\theta}-\boldsymbol{\theta^{\prime}}\|^{2},\qquad\mathbb{E}\left\|[\nabla^{2}f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta}_{n})]^{2}-H^{2}\right\|\leq
L_{p},$

where $H$ is the Hessian matrix of the population loss function $F(\cdot)$,
i.e., $H=\nabla^{2}F(\boldsymbol{\theta}^{\star})$.

######  Assumption 4.

We adopt i.i.d. random direction vectors $\\{\boldsymbol{v}_{n}\\}$ from some
common distribution $\boldsymbol{v}\sim\mathcal{P}_{\boldsymbol{v}}$ such that
$\mathbb{E}[\boldsymbol{v}\boldsymbol{v}^{\top}]=I_{d}$. Moreover, assume that
the $(6+3\delta)$-th moment of $\boldsymbol{v}$ is bounded.

We discuss the above assumptions and compare them with the standard conditions
in the literature of (RM)-type SGD inference. Assumption 1 requires the
population loss function $F(\cdot)$ to be $\lambda$-strongly convex and
$L_{f}$-smooth, which is widely assumed in the existing literature of
statistical inference on stochastic optimization (Polyak and Juditsky, 1992;
Chen et al., 2020). Note that it is possible to replace this assumption with a
weaker one that only assumes local strong convexity in the neighborhood of the
true parameter $\boldsymbol{\theta}^{\star}$(Su and Zhu, 2018; Duchi and Ruan,
2021). This weaker condition satisfies the setting of logistic regression
(Example 2.2). To highlight the main contributions of the paper, we discuss
its relaxed form of Assumption 1 in Appendix A.2. Assumption 2 introduces the
unbiasedness condition on the stochastic gradient $\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})$ when the individual loss function
$f(\boldsymbol{\theta};\boldsymbol{\zeta})$ is smooth. The $(2+\delta)$-th
moment condition is the classical Lyapunov condition used in the derivation of
asymptotic normality. Relaxation to this assumption can be made to handle
nonsmooth loss functions $f(\boldsymbol{\theta};\boldsymbol{\zeta})$, such as
the quantile regression as described in Example 2.3. We defer the discussion
of this weaker assumption to Appendix A.2. The statements in Assumption 3
introduce the Lipschitz continuity condition and the concentration condition
on the Hessian matrix. Assumption 4 guarantees that the (KW) gradient
$\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})$ is an
asymptotically unbiased estimator of $\nabla F(\boldsymbol{\theta})$ when the
spacing parameter $h_{n}$ decreases to $0$, as suggested by Lemma 2.4. The
moment condition of $\boldsymbol{v}$ in Assumption 4 is imposed for technical
simplicity and could be possibly weakened. We provide several examples of
$\mathcal{P}_{\boldsymbol{v}}$ in Section 3.1.

Before we derive the asymptotic distribution for (AKW), we first provide a
finite sample error bound for the final (KW) iterate
$\boldsymbol{\theta}_{n}$:

######  Proposition 3.1.

Assume Assumptions 1, 2, and 4 hold. Set the step size as
$\eta_{n}=\eta_{0}n^{-\alpha}$ for some constant $\eta_{0}>0$ and
$\alpha\in\left(\frac{1}{2},1\right)$ and the spacing parameter as
$h_{n}=h_{0}n^{-\gamma}$ for some constant $h_{0}>0$, and
$\gamma\in\left(\frac{1}{2},1\right)$. The (KW) iterate
$\boldsymbol{\theta}_{n}$ converges to $\boldsymbol{\theta}^{\star}$ almost
surely. Furthermore, for sufficiently large $n$, we have for $0<\delta\leq 2$,

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}-\boldsymbol{\theta}^{\star}\|^{2+\delta}$
$\displaystyle\leq Cn^{-\alpha(2+\delta)/2},$ (12)

where the constant $C$ depends on $d,\lambda,L_{f},\alpha,\gamma,$
$\eta_{0},h_{0}$.

The proof of Proposition 3.1 and the explicit dependency of the constant $C$
in (12) on the parameters and the initial value $\boldsymbol{\theta}_{0}$ are
provided in Remark A.1. A similar error bound is given by Duchi et al. (2015)
in terms of the function values for $\delta=0$. We generalize the result to
the $(2+\delta)$-moment error bound on the parameter $\boldsymbol{\theta}$,
where $\delta\in(0,2]$ is assumed in Assumption 2 for the purpose of
derivation of asymptotic normality. Proposition 3.1 suggests that the
asymptotic rate of the (KW) estimator matches the best convergence rate of the
(RM) estimator (Moulines and Bach, 2011) when the spacing parameter
$h_{n}=h_{0}n^{-\gamma}$ is a decreasing sequence with
$\gamma\in(\frac{1}{2},1)$.

Recall that to characterize the asymptotic behavior of (RM) iterates, we
denote by $S$, the Gram matrix of $\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})$ at the true parameter
$\boldsymbol{\theta}^{\star}$, i.e., $S:=\mathbb{E}\left[\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}\right]$. Analogously,
we define the limiting Gram matrix of the (KW) gradient estimator
$\widehat{g}_{h,\boldsymbol{v}}$ at $\boldsymbol{\theta}^{\star}$ as
$h\rightarrow 0$ to be $Q$. The following lemma proves that the limiting Gram
matrix takes the form of
$Q=\mathbb{E}\left[\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\right]$,
and it quantifies the distance between
$\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}$
and $Q$, as the spacing parameter $h\rightarrow 0$.

######  Lemma 3.2.

Under Assumptions 1 to 4, we have

$\displaystyle\left\|\mathbb{E}\big{[}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}\big{]}-Q\right\|\leq
Ch(1+h^{2}),\;\;Q=\mathbb{E}\big{[}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\big{]}.$

where $S=\mathbb{E}\left[\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\nabla
f(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}\right]$ is defined in
Assumption 2.

With Lemma 3.2 in place, we state our first main result that characterizes the
limiting distribution of the averaged (AKW) iterates defined in ((AKW)).

######  Theorem 3.3.

Let Assumptions 1 to 4 hold. Set the step size as
$\eta_{n}=\eta_{0}n^{-\alpha}$ for some constant $\eta_{0}>0$ and
$\alpha\in\left(\frac{1}{2},1\right)$, and the spacing parameter as
$h_{n}=h_{0}n^{-\gamma}$ for some constant $h_{0}>0$, and
$\gamma\in\left(\frac{1}{2},1\right)$. The averaged (KW) estimator
$\overline{\boldsymbol{\theta}}_{n}$ satisfies,

$\displaystyle\sqrt{n}\;\left(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}QH^{-1}\right),\qquad\text{as}\quad
n\rightarrow\infty,$ (13)

where $H=\nabla^{2}F(\boldsymbol{\theta}^{\star})$ is the population Hessian
matrix and
$Q=\mathbb{E}\left[\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\right]$
is defined in Lemma 3.2. Here $\Longrightarrow$ represents the convergence in
distribution.

We now compare the asymptotic covariance matrix of
$\overline{\boldsymbol{\theta}}_{n}$ with that of the (RM) counterpart in (5)
555Note that the asymptotic covariance $H^{-1}SH^{-1}$ in (5) is “optimal” in
the sense that it matches the asymptotic covariance for the empirical risk
minimizer $\widehat{\boldsymbol{\theta}}^{(\mathtt{ERM})}$ without online
computation and gradient information constraint.. As one can see, the
asymptotic covariance matrix of (AKW) estimator
$\overline{\boldsymbol{\theta}}_{n}$ exhibits a similar sandwich form as the
covariance matrix of (RM), but strictly dominates the latter, regardless of
the choice of random direction vectors
$\\{\boldsymbol{v}_{1},\boldsymbol{v}_{2},\dots,\boldsymbol{v}_{n}\\}$. In
fact, it is easy to check that

$\displaystyle
H^{-1}QH^{-1}-H^{-1}SH^{-1}=H^{-1}\mathbb{E}_{\boldsymbol{v}}\left[(\boldsymbol{v}\boldsymbol{v}^{\top}-I_{d})S(\boldsymbol{v}\boldsymbol{v}^{\top}-I_{d})\right]H^{-1}\succ
0,$ (14)

which suggests the (AKW) estimator suffers an inevitable loss of efficiency
compared to the $\widehat{\boldsymbol{\theta}}^{{\tt{(RM)}}}$. In Section 3.2,
we analyze (AKW) with multiple function-value queries at each iteration. With
the price of additional per-iteration computational complexity, one is able to
improve the statistical efficiency of (AKW) and achieve the optimal asymptotic
variance $H^{-1}SH^{-1}$.

######  Remark 3.4.

To complete the distributional analysis on (KW) iterates, we also provide the
asymptotic distribution of the $n$-th iterate
$\boldsymbol{\theta}_{n}^{{\tt{(KW)}}}$ of ((KW)) without averaging. Assume
the Hessian matrix has decomposition $H=P\Lambda P^{\top}$, where $P$ is an
orthogonal matrix and $\Lambda$ is a diagonal matrix. Using the proof in
Fabian (1968), we establish the following asymptotic distribution for
$\boldsymbol{\theta}_{n}^{{\tt{(KW)}}}$,

$\displaystyle
n^{\alpha/2}(\boldsymbol{\theta}_{n}^{{\tt{(KW)}}}-\boldsymbol{\theta}^{\star})\Longrightarrow\mathcal{N}(0,\Sigma),$
(15)

where each $(k,\ell)$-th entry of the covariance matrix $\Sigma$ is,

$\displaystyle\Sigma_{k\ell}=\eta_{0}\big{(}P^{\top}QP\big{)}_{kl}\big{(}\Lambda_{kk}+\Lambda_{\ell\ell}\big{)}^{-1},\quad
1\leq k,\ell\leq d.$

Here $\eta_{0}>0$ and $\alpha\in(\frac{1}{2},1)$ are specified in the step
size $\eta_{n}=\eta_{0}n^{-\alpha}$. As $\alpha<1$, the $n$-th iterate
$\boldsymbol{\theta}_{n}^{{\tt{(KW)}}}$ without averaging converges at a
slower rate $n^{-\alpha/2}$ than that of (AKW) in Theorem 3.3.

###  3.1 Examples: choices of direction distribution

By Theorem 3.3, the asymptotic covariance matrix of (AKW) estimator,
$H^{-1}QH^{-1}$, depends on the distribution of search direction
$\mathcal{P}_{\boldsymbol{v}}$ via
$Q=\mathbb{E}[\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}]$.
In this section, we compare the asymptotic covariance matrices of the (AKW)
estimator when the random directions $\\{\boldsymbol{v}_{i}\\}_{i=1}^{n}$ are
sampled from different $\mathcal{P}_{\boldsymbol{v}}$’s. Several popular
choices of $\mathcal{P}_{\boldsymbol{v}}$ are listed as follows,

  1. (G)

Gaussian: $\boldsymbol{v}\sim\mathcal{N}(0,I)$.

  2. (S)

Spherical: $\boldsymbol{v}$ is sampled from the uniform distribution on the
sphere $\|\boldsymbol{v}\|^{2}=d$.

  3. (I)

Uniform in the canonical basis: $\boldsymbol{v}$ is sampled from
$\big{\\{}\sqrt{d}\boldsymbol{e}_{1},\sqrt{d}\boldsymbol{e}_{2},\ldots,\sqrt{d}\boldsymbol{e}_{d}\big{\\}}$
with equal probability, where
$\left\\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\dots,\boldsymbol{e}_{d}\right\\}$
is the canonical basis of $\mathbb{R}^{d}$.

It is easy to verify that the above three classical choices of
$\mathcal{P}_{\boldsymbol{v}}$ satisfy Assumption 4, among which (G) and (S)
are continuous distributions, while (I) is a discrete distribution. In
particular, (I) is a discrete uniform distribution with equal probability
among the $d$ vectors of the standard basis of Euclidean space
$\mathbb{R}^{n}$, which can be generalized in the following two forms.

  1. (U)

Uniform in an arbitrary orthonormal basis $U$: $\boldsymbol{v}_{i}$ is sampled
uniformly from
$\left\\{\sqrt{d}\boldsymbol{u}_{1},\sqrt{d}\boldsymbol{u}_{2},\right.$
$\left.\ldots,\sqrt{d}\boldsymbol{u}_{d}\right\\}$, where
$\left\\{\boldsymbol{u}_{1},\boldsymbol{u}_{2},\dots,\boldsymbol{u}_{d}\right\\}$
is an arbitrary _orthonormal basis_ of $\mathbb{R}^{d}$, i.e., the matrix
$U=(\boldsymbol{u}_{1},\boldsymbol{u}_{2},\dots,\boldsymbol{u}_{d})$ is a
$d\times d$ orthonormal matrix such that $UU^{\top}=U^{\top}U=I$.

  2. (P)

Non-uniform in the canonical basis with probability
$(p_{1},p_{2},\dots,p_{d})$:
$\boldsymbol{v}=\sqrt{1/p_{k}}\,\boldsymbol{e}_{k}$ with probability
$p_{k}>0$, for $k\in[d]$ and $\sum_{k=1}^{d}p_{k}=1$.

The following proposition provides expressions of the matrix $Q$ for the above
five choices of $\mathcal{P}_{\boldsymbol{v}}$.

######  Proposition 3.5.

Under the assumptions in Theorem 3.3, for above examples of
$\mathcal{P}_{\boldsymbol{v}}$, we have

  1. (G)

Gaussian: $Q^{\tt{(G)}}=\left(2S+\mathrm{tr}(S)I_{d}\right)$.

  2. (S)

Spherical: $Q^{\tt{(S)}}=\frac{d}{d+2}\left(2S+\mathrm{tr}(S)I_{d}\right)$.

  3. (I)

Uniform in the canonical basis: $Q^{\tt{(I)}}=d\,\mathrm{diag}(S)$.

  4. (U)

Uniform in an arbitrary orthonormal basis $U$:
$Q^{\tt{(U)}}=d\,U\mathrm{diag}(U^{\top}SU)U^{\top}$.

  5. (P)

Non-uniform in a natural coordinate basis:
$Q^{\tt{(P)}}=\mathrm{diag}(S_{11}/p_{1},S_{22}/p_{2},\ldots,S_{dd}/p_{d})$.

From Proposition 3.5, one can see that any of the above choices of
$\mathcal{P}_{\boldsymbol{v}}$ leads to a $Q^{(\cdot)}$ that strictly
dominates $S$. Take $S=I_{d}$ as an example, we have $Q^{\tt{(G)}}=(d+2)I_{d}$
and $Q^{\tt{(S)}}=Q^{\tt{(I)}}=Q^{\tt{(U)}}=dI_{d}$ and
$Q^{\tt{(P)}}=\mathrm{diag}(p_{1}^{-1},p_{2}^{-1},\dots,p_{d}^{-1})\succ
I_{d}$ where $p_{1}+p_{2}+\dots+p_{d}=1$. Several additional findings and
implications of Proposition 3.5 are discussed in Appendix A.3. To briefly
mention a few, the Gaussian direction (G) is always inferior to the spherical
direction (S). Among the rest of the choices, there is _no domination_
relationship, and different optimality criterion in the experimental design
leads to different optimal choices of $\mathcal{P}_{\boldsymbol{v}}$.

###  3.2 Multi-query extension and statistical efficiency

We now consider the (AKW) estimator using $(m+1)$ function queries
$\overline{\boldsymbol{\theta}}_{n}^{(m)}$ in (11),

$\overline{\boldsymbol{\theta}}_{n}^{(m)}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{\theta}_{i}^{(m)},\qquad\text{where
}\boldsymbol{\theta}_{i}^{(m)}=\boldsymbol{\theta}_{i-1}^{(m)}-\eta_{i}\overline{g}_{n}^{(m)}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})=\boldsymbol{\theta}_{i-1}^{(m)}-\frac{\eta_{i}}{m}\sum_{j=1}^{m}\widehat{g}_{h_{i},\boldsymbol{v}_{i}^{(j)}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i}).$

Here we first consider using the same sampling distribution across $m$ queries
and $n$ iterations. In other words, $\boldsymbol{v}_{i}^{(j)}$ is sampled
_i.i.d._ from $\mathcal{P}_{\boldsymbol{v}}$ for $i=1,2,\dots,n$ and
$j=1,2,\dots,m$. Analogous to Theorem 3.3, we present the asymptotic
distribution of the multi-query (AKW),

######  Theorem 3.6.

Under the assumptions in Theorem 3.3, the $(m+1)$-query (AKW) estimator has
the following asymptotic distribution, as $n\rightarrow\infty$,

$\displaystyle\sqrt{n}\left(\overline{\boldsymbol{\theta}}_{n}^{(m)}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}Q_{m}H^{-1}\right),\;\;\text{where
}Q_{m}=\frac{1}{m}Q+\frac{m-1}{m}S.$

Theorem 3.6 illustrates a trade-off effect between the statistical efficiency
and computational efficiency. When $m=1$ and only two queries of function
evaluations are available, Theorem 3.6 reduces to Theorem 3.3, and $Q_{m}=Q$.
Conversely, as $m\rightarrow\infty$, we have $Q_{m}\rightarrow S$. Therefore,
the asymptotic covariance of $(m+1)$-query (AKW) estimator
$\overline{\boldsymbol{\theta}}_{n}^{(m)}$ approaches the optimal covariance
$H^{-1}SH^{-1}$ as $m$ approaches infinite. Nevertheless, the algorithm
requires $m$ function-value queries at each iteration, which significantly
increases the computation complexity.

For a finite $m$, a slight revision of the sampling scheme of the direction
vectors $\\{\boldsymbol{v}_{i}^{(j)}\\}_{j=1,2,\dots,m}$ provides a remedy to
achieve a smaller and indeed optimal asymptotic covariance matrix.
Particularly at the $i$-th iteration, one may sample $m$ direction vectors
$\\{\boldsymbol{v}_{i}^{(j)}\\}_{j=1,2,\dots,m}$ from a discrete distribution
(such as (I) and (U)) _without replacement_. In such settings, the direction
vectors
$\big{\\{}\boldsymbol{v}_{i}^{(1)},\boldsymbol{v}_{i}^{(2)},\dots,\boldsymbol{v}_{i}^{(m)}\big{\\}}$
are no longer independent but they have the same marginal distribution. The
asymptotic distribution of the multi-query (KW) algorithm sampling without
replacement is provided in the following theorem of its asymptotic
distribution.

######  Theorem 3.7.

Under the assumptions in Theorem 3.3, and the direction vectors in all
iterations $\big{\\{}\widetilde{V}_{i}\big{\\}}_{i=1}^{n}$ are _i.i.d._ from
$\mathcal{P}_{\boldsymbol{v}}$ such that
$\widetilde{V}_{i}=\big{(}\boldsymbol{v}_{i}^{(1)},\boldsymbol{v}_{i}^{(2)},\dots,\boldsymbol{v}_{i}^{(m)}\big{)}$
follows discrete sampling scheme in (I) and (U) WithOut Replacement (WOR), the
$(m+1)$-query (AKW) estimator, referred to as
$\overline{\boldsymbol{\theta}}_{n}^{(m,\tt{WOR})}$, has the following
asymptotic distribution, as $n\rightarrow\infty$,

$\displaystyle\sqrt{n}\left(\overline{\boldsymbol{\theta}}_{n}^{(m,\tt{WOR})}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}Q_{m}^{(\tt{WOR})}H^{-1}\right),\;\;\text{
where }Q_{m}^{(\tt{WOR})}=\frac{(d-m)}{m(d-1)}Q+\frac{d(m-1)}{m(d-1)}S.$

By comparing the asymptotic covariance matrices in Theorems 3.6 and 3.7,
$Q_{m}^{(\tt{WOR})}$ for sampling without replacement case is strictly smaller
than $Q_{m}$ in Theorems 3.6 when we consider multi-query evaluation ($m\geq
2$). Moreover, when $m=d$, it is easy to see that $Q_{m}^{(\tt{WOR})}=S$.
Therefore, the $(d+1)$-query (AKW) estimator
$\overline{\boldsymbol{\theta}}_{n}^{(m,\tt{WOR})}$ achieves the same limiting
covariance as that of the averaged (RM) estimator. Furthermore, when the model
is well-specified, the limiting covariance matrix $H^{-1}SH^{-1}=H^{-1}$
achieves the Cramér-Rao lower bound. This result indicates that the
$(d+1)$-query (AKW) estimator $\overline{\boldsymbol{\theta}}_{n}$ is
asymptotically efficient (van der Vaart, 2000).

##  4 Online Statistical Inference

In the previous section, we provide the asymptotic distribution for the (AKW)
estimator. For the purpose of conducting statistical inference of
$\boldsymbol{\theta}^{\star}$, we need a consistent estimator of the limiting
covariance $H^{-1}QH^{-1}$ in (13). A direct way is to construct a pair of
consistent estimators $\widehat{H}$ and $\widehat{Q}$ of $H$ and $Q$,
respectively, and estimate the asymptotic covariance by the _plug-in_
estimator $\widehat{H}^{-1}\widehat{Q}\widehat{H}^{-1}$. Offline construction
of those estimators is generally straightforward. However, as the (KW) scheme
typically applies to sequential data, it is ideal to estimate the asymptotic
covariance in an online fashion without storing the data. Therefore, one
cannot simply replace the true parameter $\boldsymbol{\theta}^{\star}$ by its
estimate $\overline{\boldsymbol{\theta}}_{n}$ in $Q$ and $H$ in an online
setting, since we can no longer access the data stream
$\\{\boldsymbol{\zeta}_{i}\\}_{i=1}^{n}$ after the estimator
$\overline{\boldsymbol{\theta}}_{n}$ is obtained. To address this challenge,
we first propose the following finite-difference Hessian estimator at each
iteration $n$:

$\displaystyle\widetilde{G}_{n}=\sum_{k=1}^{d}\sum_{\ell=1}^{d}\widetilde{G}_{n,kl}\boldsymbol{e}_{k}\boldsymbol{e}_{\ell}^{\top}=\frac{1}{h_{n}^{2}}\sum_{k=1}^{d}\sum_{\ell=1}^{d}\left[\Delta_{h_{n},\boldsymbol{e}_{k}}f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{e}_{\ell};\boldsymbol{\zeta}_{n})-\Delta_{h_{n},\boldsymbol{e}_{k}}f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})\right]\boldsymbol{e}_{k}\boldsymbol{e}_{\ell}^{\top},$
(16)

This construction can be viewed as a multi-query (with $d^{2}+1$ queries of
function values at each iteration) (KW) scheme with the (I) choice of the
random directions. Other choices of the search directions can be used as well,
and discussions are provided in Appendix B.1. Each additional function-value
query beyond the first one provides an estimate $\widetilde{G}_{n,kl}$ for the
$(k,l)$-th entry of the matrix $\widetilde{G}_{n}$. To reduce the
computational cost in $\widetilde{G}_{n}$, at each iteration, the algorithm
may compute a random subset of entries of $\widetilde{G}_{n}$ and partially
inhere the remaining entries from the previous estimator
$\widetilde{G}_{n-1}$. For example, each entry $\widetilde{G}_{n,k\ell}$ is
updated with probability $p\in(0,1]$. The procedure thus requires
$\mathcal{O}(pd^{2})$ function-value queries at each step. If we set
$p=\mathcal{O}(1/d^{2})$, then the query complexity is reduced to
$\mathcal{O}(1)$ per step. Since the construction of (16) does not guarantee
symmetry, an additional symmetrization step needs to be conducted, as

$\displaystyle\widetilde{H}_{n}=\frac{1}{n}\sum_{i=1}^{n}\frac{\widetilde{G}_{i}+\widetilde{G}_{i}^{\top}}{2}.$
(17)

The next lemma quantifies the estimation error of the Hessian estimator
$\widetilde{H}_{n}$ in (17) and the proof is provided in Appendix B.

######  Lemma 4.1.

Under the assumptions in Theorem 3.3, we have

$\displaystyle\mathbb{E}\|\widetilde{H}_{n}-H\|^{2}\leq
C_{1}n^{-\alpha}+C_{2}p^{-1}n^{-1}.$ (18)

From Lemma 4.1, as $n\rightarrow\infty$, the error rate is dominated by the
$C_{1}n^{-\alpha}$ term, where $\alpha$ is the parameter of the decaying step
sizes.

######  Remark 4.2.

In construction of the estimator of the limiting covariance matrix
$H^{-1}QH^{-1}$, it is necessary to avoid the possible singularity of
$\widetilde{H}_{n}$. A common practice is to adopt a thresholding version of
$\widetilde{H}_{n}$ in (17). Let $U\widetilde{\Lambda}_{n}U^{\top}$ be the
eigenvalue decomposition of $\widetilde{H}_{n}$, and define

$\displaystyle\widehat{H}_{n}=U\widehat{\Lambda}_{n}U^{\top},\;\;\widehat{\Lambda}_{n,kk}=\max\left\\{\kappa_{1},\widetilde{\Lambda}_{n,kk}\right\\},\;k=1,2,\dots,d,$
(19)

for any positive constant $\kappa_{1}<\lambda$ where $\lambda$ is defined in
Assumption 1. It is guaranteed by construction that $\widehat{H}_{n}$ is
strictly positive definite and thus invertible.

On the other hand, the estimator of Gram matrix $Q$ can be naturally
constructed as

$\displaystyle\widehat{Q}_{n}:=\frac{1}{n}\sum_{i=1}^{n}\widehat{g}_{h_{i},\boldsymbol{v}_{i}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})\,\widehat{g}_{h_{i},\boldsymbol{v}_{i}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})^{\top},$
(20)

where
$\widehat{g}_{h_{i},\boldsymbol{v}_{i}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})$
is the (KW) update in the $i$-th iteration obtained by (10). As both
$\widehat{H}_{n}$ in (19) and $\widehat{Q}_{n}$ in (20) can be constructed
sequentially without storing historical data666The sequence
$\widehat{Q}_{n}:=\frac{1}{n}\sum_{i=1}^{n}Q_{i}$ with
$Q_{i}=\widehat{g}_{h_{i},\boldsymbol{v}_{i}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})\,\widehat{g}_{h_{i},\boldsymbol{v}_{i}}(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})^{\top}$
can be constructed only with one-pass over the sequential data. In particular,
we could compute $\widehat{Q}_{n}$ sequentially as
$\widehat{Q}_{n}=\frac{1}{n}((n-1)\widehat{Q}_{n-1}+Q_{i})$., the final plug-
in estimator $\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}$ can
also be constructed in an online fashion. Based on Lemma 4.1, we obtain the
following consistency result of the covariance matrix estimator
$\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}$.

######  Theorem 4.3.

Assume Assumptions 1 to 4 hold for $\delta=2$. Set the step size as
$\eta_{n}=\eta_{0}n^{-\alpha}$ for some constant $\eta_{0}>0$ and
$\alpha\in\left(\frac{1}{2},1\right)$, and the spacing parameter as
$h_{n}=h_{0}n^{-\gamma}$ for some constant $h_{0}>0$, and
$\gamma\in\left(\frac{1}{2},1\right)$. We have

$\displaystyle\mathbb{E}\left\|\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}-H^{-1}QH^{-1}\right\|\leq
Cn^{-\alpha/2}.$

We defer the technical proof to Appendix B. Theorem 4.3 establishes the
consistency and the rate of the convergence of our proposed covariance matrix
estimator $\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}$. Given
Theorems 3.3 and 4.3, a confidence interval of the projected true parameter
$\boldsymbol{w}^{\top}\boldsymbol{\theta}^{\star}$ for any
$\boldsymbol{w}\in\mathbb{R}^{d}$ can be constructed via a projection of
$\overline{\boldsymbol{\theta}}_{n}$ and
$\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}$ onto
$\boldsymbol{w}$. Specifically, for a pre-specified confidence level $q$ and
the corresponding $z$-score $z_{q/2}$, we can obtain an asymptotic exact
confidence interval as $n\rightarrow\infty$,

$\displaystyle\mathbb{P}\left\\{\boldsymbol{w}^{\top}\boldsymbol{\theta}^{\star}\in\left[\boldsymbol{w}^{\top}\overline{\boldsymbol{\theta}}_{n}-\frac{z_{q/2}}{\sqrt{n}}\sqrt{\boldsymbol{w}^{\top}\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}\boldsymbol{w}},\;\;\;\boldsymbol{w}^{\top}\overline{\boldsymbol{\theta}}_{n}+\frac{z_{q/2}}{\sqrt{n}}\sqrt{\boldsymbol{w}^{\top}\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}\boldsymbol{w}}\right]\right\\}\rightarrow
1-q.$

###  4.1 Online inference without additional function-value queries

Despite the simplicity of the plug-in approach, the proposed estimator
$\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}$ incurs additional
computational and storage cost as it requires additional function-value
queries for constructing $\widehat{H}_{n}$. It raises a natural question: _is
it possible to conduct inference only based on(KW) iterates
$\\{\boldsymbol{\theta}_{i}\\}_{i=1,2,\dots}$ without additional function-
value queries_?

In this section, we provide an affirmative answer to this question, and
propose an alternative online statistical inference procedure using the
intermediate (KW) iterates only, without requiring any additional function-
value query. Intuitively, the (AKW) estimator in ((AKW)) is constructed as the
average of all intermediate (KW) iterates
$\\{\boldsymbol{\theta}_{i}\\}_{i=1}^{n}$. If all iterates were independent
and identically distributed, the asymptotic covariance could have been
directly estimated by the sample covariance matrix of the iterates
$\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{\theta}_{i}-\overline{\boldsymbol{\theta}})(\boldsymbol{\theta}_{i}-\overline{\boldsymbol{\theta}})^{\top}$.
Unfortunately, the (KW) iterates are far from independent and indeed highly
correlated. Nevertheless, the autocorrelation structure of the iterates can be
carefully analyzed and utilized to construct the estimator of $H^{-1}QH^{-1}$.

In this paper, we adopt an alternative approach to take more advantage of the
autocorrelation structure by leveraging the techniques from robust testing
literature (Abadir and Paruolo, 1997; Kiefer et al., 2000; Lee et al., 2022).
Such an estimator is often referred to as the Fixed Bandwidth
Heteroskedasticity and Autocorrelation Robust estimator (_fixed- $b$_ HAR) in
the econometrics literature. The _fixed- $b$_ HAR estimator is able to
overcome the series correlation and heteroskedasticity in the error terms for
the OLS estimates of the linear regression (e.g. Kiefer et al. (2000)). For
the (RM) scheme, Lee et al. (2022) utilized and generalized this technique to
construct an online statistical inference procedure, and refer to this method
as the _random scaling_ method.

In particular, we present the following theorem based on a functional
extension of the distributional analysis of the intermediate (KW) iterates
$\\{\boldsymbol{\theta}_{t}\\}$ as a stochastic process.

######  Theorem 4.4.

For any $\boldsymbol{w}\in\mathbb{R}^{d}$, under the assumptions in Theorem
3.3, we have

$\displaystyle\sqrt{n}\frac{\boldsymbol{w}^{\top}(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star})}{\sqrt{\boldsymbol{w}^{\top}V_{n}\boldsymbol{w}}}\Longrightarrow\frac{W_{1}}{\sqrt{\int_{0}^{1}\left(W_{r}-rW_{1}\right)^{2}\,\mathrm{d}r}},$
(21)

where
$V_{n}=\frac{1}{n^{2}}\sum_{i=1}^{n}i^{2}(\overline{\boldsymbol{\theta}}_{i}-\overline{\boldsymbol{\theta}}_{n})(\overline{\boldsymbol{\theta}}_{i}-\overline{\boldsymbol{\theta}}_{n})^{\top}$,
and
$\overline{\boldsymbol{\theta}}_{i}=\frac{1}{i}\sum_{\ell=1}^{i}\boldsymbol{\theta}_{\ell}$
is the average of iterates up to the $i$-th iteration, and $\\{W_{t}\\}_{t\geq
0}$ is the standard one-dimensional Brownian motion.

As an important special case, when $\boldsymbol{w}=\boldsymbol{e}_{k}$ for
$k=1,2,\ldots,d$, we have the convergence in each coordinate to the following
pivotal limiting distribution,

$\displaystyle\frac{\sqrt{n}(\overline{\boldsymbol{\theta}}_{n,k}-\boldsymbol{\theta}^{\star}_{k})}{\sqrt{V_{n,kk}}}\Longrightarrow\frac{W_{1}}{\sqrt{\int_{0}^{1}\left(W_{r}-rW_{1}\right)^{2}\,\mathrm{d}r}},$
(22)

For the asymptotic distribution defined on the right hand side in (22), we
repeat the quantiles of the distribution published by Abadir and Paruolo
(1997) in Table 1777Since the distribution on the right hand side of (21) is
symmetric, we only provide one-side quantiles in the table. . Combining the
asymptotic results in (22) and Table 1, we can construct coordinate-wise
confidence intervals for the true parameter $\boldsymbol{\theta}^{\star}$. In
addition, as

$\displaystyle
V_{n}=\frac{1}{n^{2}}\sum_{i=1}^{n}i^{2}(\overline{\boldsymbol{\theta}}_{i}-\overline{\boldsymbol{\theta}}_{n})(\overline{\boldsymbol{\theta}}_{i}-\overline{\boldsymbol{\theta}}_{n})^{\top}=\frac{1}{n^{2}}\sum_{i=1}^{n}i^{2}\overline{\boldsymbol{\theta}}_{i}\overline{\boldsymbol{\theta}}_{i}^{\top}-\frac{2}{n^{2}}\overline{\boldsymbol{\theta}}_{n}\sum_{i=1}^{n}i^{2}\overline{\boldsymbol{\theta}}_{i}^{\top}+\frac{1}{n^{2}}\sum_{i=1}^{n}i^{2}\overline{\boldsymbol{\theta}}_{n}\overline{\boldsymbol{\theta}}_{n}^{\top}$
(23)

can be constructed in an online fashion via the iterative updates of the
matrix
$\sum_{i=1}^{n}i^{2}\overline{\boldsymbol{\theta}}_{i}\overline{\boldsymbol{\theta}}_{i}^{\top}$
and the vector $\sum_{i=1}^{n}i^{2}\overline{\boldsymbol{\theta}}_{i}$, the
proposed online inference procedure only requires one pass over the data.

Quantile | 90% | 95% | 97.5% | 99%
---|---|---|---|---
Abadir and Paruolo (1997) Table 1 | 3.875 | 5.323 | 6.747 | 8.613
Table 1: Cumulative probability table of the limiting distribution.

##  5 Numerical Experiments

In this numerical section, we first investigate the empirical performance of
the proposed inference procedures and their corresponding coverage rates. We
consider linear regression and logistic regression models (Examples 2.1\--2.2)
where $\\{\boldsymbol{x}_{i},y_{i}\\}_{i=1}^{n}$ is an _i.i.d._ sample with
the covariate $\boldsymbol{x}\sim\mathcal{N}(\mathbf{0},\Sigma)$ and the
response $y\in\mathbb{R}$. The true model parameter
$\boldsymbol{\theta}^{\star}\in\mathbb{R}^{d}$ is selected uniformly from the
unit sphere before the experiments. For both models, we consider two different
structures of the covariance matrices $\Sigma$: identity matrix $I_{d}$ and
equicorrelation covariance matrix (Equicorr in the tables), i.e.,
$\Sigma_{k\ell}=0.2$ for all $k\neq\ell$ and $\Sigma_{kk}=1$. The parameter
$\alpha$ in the step size is specified to $\alpha=0.501$. The variance of
noise $\epsilon$ in the linear regression model (Example 2.1) is set to
$\sigma^{2}=0.2$.

| $d$ | $\Sigma$ | Estimation error | Average coverage rate | Average length
---|---|---|---|---|---
|  |  | Parameter | Plug-in Cov. | Plug-in |  Fixed-$b$ | Oracle | Plug-in |  Fixed-$b$ | Oracle
Linear regression |  |  |  |  |  |  |  |
|  |  | 0.0031 | 0.0384 | 0.9448 | 0.9464 | 0.9436 | 1.7555 | 2.1915 | 1.7533
|  | Identity | (0.0010) | (0.0106) | (0.1035) | (0.1174) | (0.1040) | (0.0082) | (0.4184) | -
|  |  | 0.0035 | 0.0342 | 0.9428 | 0.9488 | 0.9412 | 2.0109 | 2.4895 | 2.0078
| 5 | Equicorr | (0.0012) | (0.0092) | (0.1096) | (0.1195) | (0.1102) | (0.0097) | (0.5323) | -
|  |  | 0.0135 | 0.1126 | 0.9319 | 0.9039 | 0.9288 | 3.5337 | 3.7424 | 3.5065
|  | Identity | (0.0023) | (0.0190) | (0.0594) | (0.0657) | (0.0616) | (0.0164) | (0.4292) | -
|  |  | 0.0172 | 0.1124 | 0.9194 | 0.9014 | 0.9170 | 4.3140 | 4.5582 | 4.2753
| 20 | Equicorr | (0.0029) | (0.0199) | (0.0644) | (0.0674) | (0.0656) | (0.0207) | (0.5681) | -
|  |  | 0.0748 | 0.5707 | 0.9309 | 0.7501 | 0.9012 | 8.6675 | 7.6819 | 7.8397
|  | Identity | (0.0062) | (0.0648) | (0.0261) | (0.0397) | (0.0336) | (0.1081) | (0.4924) | -
|  |  | 0.0921 | 0.5615 | 0.9331 | 0.7435 | 0.9044 | 10.7701 | 9.0193 | 9.7508
| 100 | Equicorr | (0.0076) | (0.0647) | (0.0250) | (0.0418) | (0.0320) | (0.1400) | (0.6161) | -
Logistic regression |  |  |  |  |  |  |  |
|  |  | 0.0265 | 0.0587 | 0.9432 | 0.9360 | 0.9440 | 3.1136 | 3.3775 | 3.1078
|  | Identity | (0.0115) | (0.0434) | (0.1219) | (0.1685) | (0.1148) | (0.0936) | (0.2074) | -
|  |  | 0.0299 | 0.0697 | 0.9440 | 0.9364 | 0.9464 | 3.3620 | 3.8028 | 3.2580
| 5 | Equicorr | (0.0131) | (0.0514) | (0.1196) | (0.1566) | (0.1207) | (0.1057) | (0.2558) | -
|  |  | 0.0728 | 0.1030 | 0.9418 | 0.8956 | 0.9403 | 4.8751 | 5.1763 | 4.8374
|  | Identity | (0.0124) | (0.0250) | (0.0532) | (0.1156) | (0.0540) | (0.1973) | (0.4362) | -
|  |  | 0.0799 | 0.1213 | 0.9383 | 0.8949 | 0.9369 | 5.6873 | 5.7532 | 5.6356
| 20 | Equicorr | (0.0146) | (0.0359) | (0.0577) | (0.1106) | (0.0561) | (0.1715) | (0.4064) | -
|  |  | 0.2440 | 0.5236 | 0.9673 | 0.7022 | 0.9082 | 12.0661 | 8.7892 | 10.3175
|  | Identity | (0.0211) | (0.1646) | (0.0193) | (0.0838) | (0.0295) | (0.4642) | (0.6000) | -
|  |  | 0.2867 | 0.7685 | 0.9608 | 0.6950 | 0.9041 | 12.7375 | 9.4884 | 10.4868
| 100 | Equicorr | (0.0253) | (0.2933) | (0.0185) | (0.0728) | (0.0314) | (0.6870) | (0.6170) | -
Table 2: Estimation errors, averaged coverage rates, and average lengths of
the proposed algorithm with search direction (I) and two function queries
($m=1$). Sample size $n=10^{5}$. Corresponding standard errors are reported in
the brackets. We compare the plug-in covariance estimator (plug-in) based
inference (17) and fixed-$b$ HAR (fixed-$b$) based inference (22).

Figure 3: Convergence of the parameter estimation error
$\|\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\|$ and
coverage rates v.s. the sample size $n$ when $d=20$ and the population design
matrix $\Sigma=I$. Plots (a) to (b) show the cases of linear regression and
plots (c) to (d) show the cases of logistic regression. Dashed lines in plots
(b) and (d) correspond to the nominal $95\%$ coverage.

###  5.1 Estimation errors of (AKW) and the performance of inference
procedures

We set the sample size $n=10^{5}$ and the parameter dimension $d=5,20,100$. We
first report the performance of (AKW) with the search direction uniformly
sampled from the natural basis, referred to as (I) in Section 3.1. In Table 2,
we present the estimation error for the parameter
$\boldsymbol{\theta}^{\star}$ in the Euclidean norm and the relative error of
the plug-in covariance estimator in the spectral norm (see the first two
columns), with $100$ Monte-Carlo simulations,

$\displaystyle\frac{\|\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\|}{\|\boldsymbol{\theta}^{\star}\|},\;\;\frac{\|\widehat{H}_{n}^{-1}\widehat{Q}_{n}\widehat{H}_{n}^{-1}-H^{-1}QH^{-1}\|}{\|H^{-1}QH^{-1}\|}.$
(24)

Next, we set the nominal coverage probability as $95\%$ and we project
$\boldsymbol{\theta}\in\mathbb{R}^{d}$ onto
$\boldsymbol{w}=(1,1,\dots,1)^{\top}/\sqrt{d}$ to construct confidence
intervals. In particular, we report the performance of the confidence interval
with the average coverage rate and the average length of the intervals for (1)
the plug-in covariance matrix estimator888Here we use updating probability
$p=1$ for the plug-in estimation. In other words, $d^{2}+1$ queries of
function values are obtained at each iteration. In Section 5.2 below, we
extend the comparison for different $p$. (16) and (2) the fixed-$b$ HAR
procedure in (23). As an oracle benchmark, we also report the length of the
confidence interval with respect to the true covariance matrix $H^{-1}QH^{-1}$
of the plug-in approach and the corresponding coverage rate. As shown from
Table 2, the coverage rate of the plug-in covariance estimator and the oracle
coverage rates are very close to the desired $95\%$ coverage, while the
fixed-$b$ HAR approach is comparable in small dimension $d=5,20$ but has lower
coverage rates for the large dimension $d=100$. The average lengths of both
methods are comparable to the lengths derived by the true limiting covariance.

Then, we fix $d=20$ and the identity design matrix $\Sigma=I$ . We present in
Figure 3 the parameter estimation error
$\|\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\|$ as the
sample size $n$ grows (see subplots (a) and (c) in Figure 3, for linear and
logistic regression, respectively). In subplots (b) and (d) of Figure 3, we
show the coverage rates for the plug-in and fixed-$b$ HAR approaches as the
sample size $n$ increases. As one can see, coverage rates of the plug-in
approach almost match the oracle case using the true asymptotic covariance
matrix $H^{-1}QH^{-1}$. For the linear regression case, the plug-in and
fixed-$b$ HAR approaches are comparable. For the logistic regression case, the
coverage rate of the fixed-$b$ HAR inference procedure is slightly inferior
than that of the plug-in method. On the other hand, the fixed-$b$ HAR approach
does not require additional function queries for the explicit estimation of
the Hessian matrix. Additional simulation results for the equicorrelation
design are relegated to Figure C.1 in the Appendix.

$\Sigma$ | Estimator | Comp. | Estimation error | Average coverage rate | Average length
---|---|---|---|---|---
|  | time | Hessian | Cov. | Estimator | Oracle | Estimator | Oracle
Identity |  |  |  |  |  |  |
|  |  | 0.1780 | 0.3179 | 0.8965 | 0.9288 | 3.6570 | 3.5065
| Plug-in $p=1/400$ | 4.74s | (0.0115) | (0.0423) | (0.0696) | (0.0616) | (0.0195) | -
|  |  | 0.0393 | 0.1503 | 0.9244 | 0.9288 | 3.5511 | 3.5065
| Plug-in $p=1/20$ | 25.42s | (0.0043) | (0.0282) | (0.0665) | (0.0616) | (0.0169) | -
|  |  | 0.0271 | 0.1126 | 0.9319 | 0.9288 | 3.5337 | 3.5065
| Plug-in $p=1$ | 510.53s | (0.0021) | (0.0190) | (0.0594) | (0.0616) | (0.0164) | -
|  |  | - | - | 0.9039 | 0.9288 | 3.7424 | 3.5065
| Fixed-$b$ | 2.82s | - | - | (0.0657) | (0.0616) | (0.4292) | -
Equicorr |  |  |  |  |  |  |
|  |  | 0.0381 | 0.4211 | 0.8815 | 0.9170 | 4.4547 | 4.2753
| Plug-in $p=1/400$ | 4.78s | (0.0043) | (0.0421) | (0.0753) | (0.0656) | (0.0304) | -
|  |  | 0.0117 | 0.1540 | 0.9122 | 0.9170 | 4.3489 | 4.2753
| Plug-in $p=1/20$ | 25.60s | (0.0025) | (0.0271) | (0.0691) | (0.0656) | (0.0293) | -
|  |  | 0.0082 | 0.1124 | 0.9194 | 0.9170 | 4.3140 | 4.2753
| Plug-in $p=1$ | 512.07s | (0.0018) | (0.0199) | (0.0644) | (0.0656) | (0.0207) | -
|  |  | - | - | 0.9014 | 0.9170 | 4.5582 | 4.2753
| Fixed-$b$ | 2.85s | - | - | (0.0674) | (0.0656) | (0.5681) | -
Table 4: Computation time, estimation errors, averaged coverage rates, and average lengths of the proposed algorithm with search direction (I) and two function queries ($m=1$). Sample size $n=10^{5}$, dimension $d=20$ under the linear regression model. Corresponding standard errors are reported in the brackets. We compare the plug-in covariance estimator (plug-in) based inference (19) using $p=1,1/20,1/400$ and fixed-$b$ HAR (fixed-$b$) based inference (23). $\Sigma$ | Estimator | Comp. | Estimation error | Average coverage rate | Average length
---|---|---|---|---|---
|  | time | Hessian | Cov. | Estimator | Oracle | Estimator | Oracle
Identity |  |  |  |  |  |  |
|  |  | 0.1812 | 0.3293 | 0.9039 | 0.9403 | 5.0715 | 4.8374
| Plug-in $p=1/400$ | 5.70s | (0.0281) | (0.0792) | (0.0501) | (0.0540) | (0.2071) | -
|  |  | 0.0737 | 0.1636 | 0.9330 | 0.9403 | 4.9599 | 4.8374
| Plug-in $p=1/20$ | 32.32s | (0.0114) | (0.0393) | (0.0593) | (0.0540) | (0.1833) | -
|  |  | 0.0597 | 0.1030 | 0.9418 | 0.9403 | 4.8751 | 4.8374
| Plug-in $p=1$ | 643.86s | (0.0093) | (0.0250) | (0.0532) | (0.0540) | (0.1973) | -
|  |  | - | - | 0.8956 | 0.9403 | 5.1763 | 4.8374
| Fixed-$b$ | 3.13s | - | - | (0.1156) | (0.0540) | (0.4362) | -
Equicorr |  |  |  |  |  |  |
|  |  | 0.0993 | 0.3620 | 0.8880 | 0.9369 | 5.9456 | 5.6356
| Plug-in $p=1/400$ | 5.75s | (0.0287) | (0.0992) | (0.0540) | (0.0561) | (0.1716) | -
|  |  | 0.0356 | 0.1441 | 0.9288 | 0.9369 | 5.7766 | 5.6356
| Plug-in $p=1/20$ | 32.53s | (0.0120) | (0.0440) | (0.0599) | (0.0561) | (0.1556) | -
|  |  | 0.0240 | 0.1213 | 0.9383 | 0.9369 | 5.6873 | 5.6356
| Plug-in $p=1$ | 645.56s | (0.0101) | (0.0359) | (0.0577) | (0.0561) | (0.1715) | -
|  |  | - | - | 0.8949 | 0.9369 | 5.7532 | 5.6356
| Fixed-$b$ | 3.17s | - | - | (0.1106) | (0.0561) | (0.4064) | -
Table 5: Computation time, estimation errors, averaged coverage rates, and
average lengths of the proposed algorithm with search direction (I) and two
function queries ($m=1$). Sample size $n=10^{5}$, dimension $d=20$ under the
logistic regression model. Corresponding standard errors are reported in the
brackets. We compare the plug-in covariance estimator (plug-in) based
inference (19) using $p=1,1/20,1/400$ and fixed-$b$ HAR (fixed-$b$) based
inference (23).

###  5.2 Comparison of the inference procedures

In this subsection, we provide detailed comparisons of different inference
procedures. Specifically, we fix dimension $d=20$, and compare the performance
for the plug-in and fixed-$b$ HAR schemes. For plug-in estimators, at each
iteration, we update the Hessian estimator $\widehat{H}_{n}$ in (19) using
(16) with probability $p$ chosen from $1,\leavevmode\nobreak\
d^{-1},\leavevmode\nobreak\ d^{-2}$. The fixed-$b$ scheme is updated by (23).
We report the computation time, the estimation error of the Hessian matrix,
and the average coverage rate and length of these candidates based on $100$
replications. The computation time is recorded in a simulation environment
running Python 3.8 with a single 10-core Apple M1 Max chip.

The simulation for linear regression and logistic regression is given below in
Tables 4\--5, respectively. As can be referred from the two tables, the
fixed-$b$ HAR approach gives the fastest execution, due to the fact that no
additional function queries are required for Hessian matrix computation. The
fixed-$b$ HAR method is even faster than the case of $p=d^{-2}$ where we
update only one entry (in expectation) for the Hessian matrix in each (KW)
step. Among plug-in cases, the performance of inference improves as $p$
increases, and it achieves a relatively more reliable coverage for $p\geq
d^{-1}$, (i.e., at least $d$ entries (in expectation) receive updates for the
Hessian estimator per (KW) step), with a significant cost of the computation
time. In practice, we would recommend the fixed-$b$ HAR method for those
computation-sensitive tasks, and the plug-in method with Hessian sampling
probability $p\geq d^{-1}$ in less computation-sensitive tasks.

$d$ | $\mathcal{P}_{\boldsymbol{v}}$ | Estimation error | Average coverage rate | Average length
---|---|---|---|---
|  | Parameter | Plug-in Cov. | Plug-in | Oracle | Plug-in | Oracle
|  | 0.0265 | 0.0587 | 0.9432 | 0.9440 | 3.1136 | 3.1078
| (I) | (0.0115) | (0.0434) | (0.1219) | (0.1148) | (0.8648) | -
|  | 0.0264 | 0.0599 | 0.9396 | 0.9376 | 3.0639 | 3.0625
| (S) | (0.0124) | (0.0453) | (0.1276) | (0.1250) | (0.8211) | -
|  | 0.0312 | 0.0718 | 0.9412 | 0.9420 | 3.6304 | 3.6237
5 | (G) | (0.0139) | (0.0498) | (0.1193) | (0.1176) | (0.9770) | -
|  | 0.0728 | 0.1030 | 0.9418 | 0.9403 | 4.8751 | 4.8374
| (I) | (0.0124) | (0.0250) | (0.0532) | (0.0540) | (0.6441) | -
|  | 0.0711 | 0.1017 | 0.9438 | 0.9419 | 4.8414 | 4.8156
| (S) | (0.0116) | (0.0246) | (0.0523) | (0.0524) | (0.6322) | -
|  | 0.0749 | 0.1054 | 0.9427 | 0.9423 | 5.0873 | 5.0507
20 | (G) | (0.0121) | (0.0248) | (0.0563) | (0.0523) | (0.6654) | -
|  | 0.2440 | 0.5236 | 0.9673 | 0.9082 | 12.0661 | 10.3175
| (I) | (0.0211) | (0.1646) | (0.0193) | (0.0295) | (1.0106) | -
|  | 0.2353 | 0.5122 | 0.9605 | 0.9145 | 13.1366 | 11.1788
| (S) | (0.0205) | (0.1530) | (0.0201) | (0.0358) | (1.0891) | -
|  | 0.2357 | 0.5147 | 0.9614 | 0.9161 | 13.2836 | 11.2901
100 | (G) | (0.0202) | (0.1531) | (0.0205) | (0.0380) | (1.0929) | -
Table 6: Comparison among different direction distributions $\mathcal{P}_{\boldsymbol{v}}$ (Detailed specification of (I),(S),(G) can be referred to Section 3.1). We consider the logistic regression model with design matrix $\Sigma=I$, and the (AKW) estimators are computed based on the case of two function queries ($m=1$). Corresponding standard errors are reported in the brackets. $m;\Sigma$ | $\mathcal{P}_{\boldsymbol{v}}$ | Estimation error | Average coverage rate | Average length
---|---|---|---|---
|  | Parameter | Plug-in Cov. | Plug-in | Oracle | Plug-in | Oracle
|  | 0.0916 | 0.1972 | 0.9547 | 0.9342 | 3.7013 | 3.4794
| (I+WOR) | (0.0103) | (0.1053) | (0.0225) | (0.0330) | (0.2970) | -
|  | 0.0947 | 0.2004 | 0.9551 | 0.9353 | 3.8800 | 3.6383
| (I+WR) | (0.0106) | (0.1025) | (0.0215) | (0.0310) | (0.3053) | -
|  | 0.0958 | 0.2134 | 0.9552 | 0.9320 | 3.8893 | 3.6352
10; Identity | (S) | (0.0118) | (0.1172) | (0.0219) | (0.0368) | (0.3054) | -
|  | 0.1184 | 0.2581 | 0.9404 | 0.9126 | 3.6432 | 3.3700
| (I+WOR) | (0.0122) | (0.1278) | (0.0252) | (0.0382) | (0.2240) | -
|  | 0.1235 | 0.2828 | 0.9431 | 0.9125 | 3.8352 | 3.5234
| (I+WR) | (0.0145) | (0.1573) | (0.0266) | (0.0437) | (0.2498) | -
|  | 0.1224 | 0.2753 | 0.9435 | 0.9135 | 3.8225 | 3.5165
10; Equicorr | (S) | (0.0144) | (0.1501) | (0.0259) | (0.0422) | (0.2614) | -
|  | 0.0261 | 0.0531 | 0.9455 | 0.9438 | 0.8978 | 0.8938
| (I+WOR) | (0.0022) | (0.0135) | (0.0297) | (0.0305) | (0.0290) | -
|  | 0.0333 | 0.0568 | 0.9455 | 0.9441 | 1.4037 | 1.3948
| (I+WR) | (0.0030) | (0.0196) | (0.0253) | (0.0262) | (0.0803) | -
|  | 0.0334 | 0.0556 | 0.9458 | 0.9439 | 1.4034 | 1.3941
100; Identity | (S) | (0.0028) | (0.0199) | (0.0231) | (0.0247) | (0.0816) | -
|  | 0.0328 | 0.0664 | 0.9490 | 0.9441 | 0.9056 | 0.8971
| (I+WOR) | (0.0035) | (0.0199) | (0.0339) | (0.0356) | (0.0493) | -
|  | 0.0453 | 0.0823 | 0.9494 | 0.9444 | 1.4183 | 1.3946
| (I+WR) | (0.0046) | (0.0322) | (0.0240) | (0.0270) | (0.0766) | -
|  | 0.0451 | 0.0821 | 0.9497 | 0.9449 | 1.4157 | 1.3930
100; Equicorr | (S) | (0.0048) | (0.0321) | (0.0249) | (0.0270) | (0.0777) | -
Table 7: Comparison among different sampling schemes for multi-query
algorithms under logistic regression model with dimension $d=100$ and
$m=10,100$, respectively (Detailed specification of (I+WOR),(I+WR),(S) can be
referred to Section 3.1). Corresponding standard errors are reported in the
brackets.

###  5.3 Choices of the search direction distribution

In this subsection, we compare the results for different directions
$\mathcal{P}_{\boldsymbol{v}}$. We report the results for the logistic
regression model with the identity design matrix $\Sigma=I$ in Table 6. Table
6 suggests the (AKW) algorithms with search directions (I), (S), (G) achieve
similar performance for parameter estimation error and average coverage rates,
while the average confidence intervals of (G) are generally larger. The
observations in the numerical experiments match our Proposition 3.5.
Additional simulation results of the linear regression model and the
equicorrelation design are relegated to Tables C.1, C.2, C.3 in the Appendix.

###  5.4 Multi-query (AKW) estimator

We further conduct experiments for the (KW) algorithm with multiple function-
value queries ($m>1$) and compare the performance of $m=10,100$ using
different search directions with sampling schemes (I+WR), (I+WOR), and (S). We
note that (I+WR) and (I+WOR) refer to the uniform sampling from natural basis
with and without replacement, respectively; and (S) refers to the uniform
sampling from the sphere. We report the results of the logistic regression
model in Table 7 and relegate the results of the linear regression to Table
C.4 in the Appendix.

When $m=10$, the (KW) algorithm using all three sampling schemes achieves
similar performance in both estimation and inference. When $m=100$, the
algorithm with (I+WOR) achieves better performance than the other two sampling
schemes by constructing around 30% shorter confidence intervals on average
while achieving comparable coverage rates.

We further present in Figure 8 the estimation error of the parameters and
covariance matrices when we increase the function-query complexity $m$. The
numerical results matches the magnitudes of $Q$ with regard to different $m$
in Theorems 3.6\--3.7, which could help practitioners choose an appropriate
$m$ to balance the accuracy and computational cost. We report the logistic
regression results with the identity design matrix $\Sigma=I$ in Figure 8 and
relegate the equicorrelation design and the linear regression results to
Figures C.2\--C.3 in the Appendix.

Figure 8: The parameter estimation error and the relative covariance
estimation error (see (24)) for multiple function-value evaluations. The
$x$-axis is the number of function evaluations per step (i.e., $m+1$). Here,
we consider the logistic regression model with $n=10^{5}$ and $d=100$.

##  6 Conclusion and Future Work

In this paper, we investigate the statistical inference problem for the
Kiefer-Wolfowitz stochastic optimization algorithm with random search
directions. We show the asymptotic normality for the (KW)-type estimators and
provide consistent estimators of the asymptotic covariance matrix to
facilitate the inference. Our theoretical analysis provides a comprehensive
comparison on the impact of different random search directions, the number of
multi-query evaluations, and sampling schemes. Our findings are validated by
numerical experiments.

For future works, our results and estimation methods may be potentially useful
to understand asymptotic behaviors of other gradient-free variants of
stochastic optimization algorithms, e.g. moment-adjusted stochastic gradients
(Liang and Su, 2019), stochastic optimization under constraints (Duchi and
Ruan, 2021), high dimensional stochastic algorithms (Chao and Cheng, 2019; Shi
et al., 2021), and SGD in contextual bandit settings Chen et al. (2021).

## Appendix

##  Appendix A Supplementary Material for Section 3

Throughout the Appendix, we will assume, without loss of generality,
$F(\cdot)$ achieves its minimum at $\boldsymbol{\theta}^{\star}=\mathbf{0}$
and $F(\mathbf{0})=0$. We now introduce some notations as follows,

$\displaystyle\boldsymbol{\xi}_{n}$ $\displaystyle=\nabla
F(\boldsymbol{\theta}_{n-1})-\mathbb{E}_{n-1}(\frac{1}{h_{n}}[F(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n})-F(\boldsymbol{\theta}_{n-1})]\boldsymbol{v}_{n}),$
$\displaystyle\boldsymbol{\gamma}_{n}$
$\displaystyle=\mathbb{E}_{n-1}\frac{1}{h_{n}}[F(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n})-F(\boldsymbol{\theta}_{n-1})]\boldsymbol{v}_{n}-\frac{1}{h_{n}}[F(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n})-F(\boldsymbol{\theta}_{n-1})]\boldsymbol{v}_{n},$
$\displaystyle\boldsymbol{\varepsilon}_{n}$
$\displaystyle=\frac{1}{h_{n}}[F(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n})-F(\boldsymbol{\theta}_{n-1})]\boldsymbol{v}_{n}-\frac{1}{h_{n}}[f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n};\boldsymbol{\zeta}_{n})-f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})]\boldsymbol{v}_{n}.$

###  A.1 Two-query approximation

#### Proof of Lemma 2.4

###### Proof.

By definition,
$\mathbb{E}_{\boldsymbol{\zeta}}\;\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})=\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v}=\frac{1}{h}\left[F(\boldsymbol{\theta}+h\boldsymbol{v})-F(\boldsymbol{\theta})\right]\boldsymbol{v}$.
For the first inequality, we have

$\displaystyle\left\|\mathbb{E}\;\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta})\right\|$
$\displaystyle=\left\|\mathbb{E}\;\frac{1}{h}\left[F(\boldsymbol{\theta}+h\boldsymbol{v})-F(\boldsymbol{\theta})\right]\boldsymbol{v}-\nabla
F(\boldsymbol{\theta})\right\|$
$\displaystyle=\left\|\mathbb{E}\;\boldsymbol{v}\boldsymbol{v}^{\top}\nabla
F(\boldsymbol{\theta})+\frac{1}{2}h\mathbb{E}\;\boldsymbol{v}\boldsymbol{v}^{\top}\nabla^{2}F(\boldsymbol{\theta}_{h,\boldsymbol{v}})\boldsymbol{v}-\nabla
F(\boldsymbol{\theta})\right\|$
$\displaystyle=\frac{1}{2}h\left\|\mathbb{E}\;\boldsymbol{v}\boldsymbol{v}^{\top}\nabla^{2}F(\boldsymbol{\theta}_{h,\boldsymbol{v}})\boldsymbol{v}\right\|$
$\displaystyle\leq\frac{1}{2}hL_{f}\mathbb{E}\|\boldsymbol{v}\|^{3},$ (25)

where in the third equality we use the Taylor expansion of
$F(\boldsymbol{\theta})$, and $\boldsymbol{\theta}_{h,\boldsymbol{v}}$ comes
from the remainder term of the Taylor expansion. ∎

#### Proof of Proposition 3.1

######  Proposition 3.1.

Assume Assumptions 1, 2, and 4 hold. Set the step size as
$\eta_{n}=\eta_{0}n^{-\alpha}$ for some constant $\eta_{0}>0$ and
$\alpha\in\left(\frac{1}{2},1\right)$ and the spacing parameter as
$h_{n}=h_{0}n^{-\gamma}$ for some constant $h_{0}>0$, and
$\gamma\in\left(\frac{1}{2},1\right)$. The (KW) iterate
$\boldsymbol{\theta}_{n}$ converges to $\boldsymbol{\theta}^{\star}$ almost
surely. Furthermore, for sufficiently large $n$, we have for $0<\delta\leq 2$,

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}-\boldsymbol{\theta}^{\star}\|^{2+\delta}$
$\displaystyle\leq Cn^{-\alpha(2+\delta)/2}.$

where the constant $C$ depends on $d,\lambda,L_{f},\alpha,\gamma,$
$\eta_{0},h_{0}$.

######  Remark A.1.

The parameter dependency in Proposition 3.1 could be given explicitly as
follows,

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}-\boldsymbol{\theta}^{*}\|^{2}$
$\displaystyle\leq\exp\left(CM_{1}\eta_{0}/(2\alpha-1)+CM_{2}/(2\beta-1)-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right)\|\boldsymbol{\theta}_{0}\|^{2}$
$\displaystyle\quad+M_{3}\left(\exp\left(-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right)+\frac{\eta_{0}n^{-\alpha}}{\lambda}\right)$
$\displaystyle\quad+\frac{M_{3}}{M_{1}}\exp\left(CM_{1}\eta_{0}/(2\alpha-1)+CM_{2}/(2\beta-1)-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right),$

where the constant $C$ above is a universal constant that does not depend on
any constant/parameters in the assumptions. The other terms
$M_{1},M_{2},M_{3}$ above are given below,

$\displaystyle M_{1}$
$\displaystyle=C\left(L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{4}+M^{\frac{2}{2+\delta}}\mathbb{E}\|\boldsymbol{v}\|^{4}+L_{f}^{2}\right),$
$\displaystyle M_{2}$
$\displaystyle=CL_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{3},$ $\displaystyle
M_{3}$
$\displaystyle=C\left(\mathbb{E}\|\boldsymbol{v}\|^{3}+\left(h_{n}^{2}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{6}+M^{\frac{2}{2+\delta}}\mathbb{E}\|\boldsymbol{v}\|^{4}(h_{n}^{2}\|\boldsymbol{v}\|^{2}+1)\right)\right).$

We will prove both Proposition 3.1 and Remark A.1 below.

###### Proof.

We first give some bounds on
$\boldsymbol{\xi}_{n},\boldsymbol{\gamma}_{n},\boldsymbol{\varepsilon}_{n}$.
By definition,
$\mathbb{E}_{n-1}\boldsymbol{\gamma}_{n}=\mathbb{E}_{n-1}\boldsymbol{\varepsilon}_{n}=0$.
From (A.1),

$\displaystyle\|\boldsymbol{\xi}_{n}\|\leq\frac{1}{2}h_{n}L_{f}\mathbb{E}\|\boldsymbol{v}\|^{3}.$
(26)

We can bound $\boldsymbol{\gamma}_{n}$ by the following

$\displaystyle\mathbb{E}\|\boldsymbol{\gamma}_{n}\|^{2}$
$\displaystyle\leq\mathbb{E}\left\|\frac{1}{h_{n}}[F(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{v}_{n})-F(\boldsymbol{\theta}_{n-1})]\boldsymbol{v}_{n}\right\|^{2}$
$\displaystyle\leq\mathbb{E}\|\langle\nabla
F(\boldsymbol{\theta}_{n-1}),\boldsymbol{v}_{n}\rangle\boldsymbol{v}_{n}\|^{2}+\frac{1}{4}h_{n}^{2}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{6}$
$\displaystyle\leq
L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{4}\mathbb{E}\|\boldsymbol{\theta}_{n-1}\|^{2}+\frac{1}{4}h_{n}^{2}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{6}.$
(27)

We also have the following fact for $\boldsymbol{\varepsilon}$.

$\displaystyle\mathbb{E}_{n-1}\left[\|\boldsymbol{\varepsilon}_{n}\|^{2}|\boldsymbol{v}_{n}\right]$
$\displaystyle=\;\;$
$\displaystyle\mathbb{E}_{n-1}\left[\left\|\frac{1}{h_{n}}\int_{0}^{h_{n}}\langle\nabla
F(\boldsymbol{\theta}_{n-1}+s\boldsymbol{v}_{n})-\nabla
f(\boldsymbol{\theta}_{n-1}+s\boldsymbol{v}_{n};\boldsymbol{\zeta}_{n}),\boldsymbol{v}_{n}\rangle\boldsymbol{v}_{n}\,\mathrm{d}s\right\|^{2}\bigg{|}\boldsymbol{v}_{n}\right]$
$\displaystyle\leq\;\;$
$\displaystyle\|\boldsymbol{v}_{n}\|^{4}\mathbb{E}_{n-1}\left[\frac{1}{h_{n}}\int_{0}^{h_{n}}\left\|\nabla
F(\boldsymbol{\theta}_{n-1}+s\boldsymbol{v}_{n})-\nabla
f(\boldsymbol{\theta}_{n-1}+s\boldsymbol{v}_{n};\boldsymbol{\zeta}_{n})\right\|^{2}\,\mathrm{d}s\big{|}\boldsymbol{v}_{n}\right]$
$\displaystyle\leq\;\;$ $\displaystyle
M^{\frac{2}{2+\delta}}\|\boldsymbol{v}_{n}\|^{4}\frac{1}{h_{n}}\int_{0}^{h_{n}}(\|\boldsymbol{\theta}_{n-1}+s\boldsymbol{v}_{n}\|^{2}+1)\,\mathrm{d}s$
$\displaystyle\leq\;\;$ $\displaystyle
M^{\frac{2}{2+\delta}}\|\boldsymbol{v}_{n}\|^{4}(\|\boldsymbol{\theta}_{n-1}\|^{2}+h_{n}^{2}\|\boldsymbol{v}_{n}\|^{2}+1),$
(28)

where in the second inequality, we use Assumption 2.

Now decompose the update step as follows,

$\displaystyle\boldsymbol{\theta}_{n}$
$\displaystyle=\boldsymbol{\theta}_{n-1}-\eta_{n}\frac{1}{h_{n}}[f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{u}_{n};\boldsymbol{\zeta}_{n})-f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})]$
$\displaystyle=\boldsymbol{\theta}_{n-1}-\eta_{n}\nabla
F(\boldsymbol{\theta}_{n-1})+\eta_{n}(\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}).$

Therefore, we can derive that,

$\displaystyle\|\boldsymbol{\theta}_{n}\|^{2}\leq$
$\displaystyle\|\boldsymbol{\theta}_{n-1}\|^{2}-2\eta_{n}\langle\nabla
F(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n-1}\rangle+2\eta_{n}\langle\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n},\boldsymbol{\theta}_{n-1}\rangle$
$\displaystyle+\eta_{n}^{2}\|\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}-\nabla
F(\boldsymbol{\theta}_{n-1})\|^{2}.$ (29)

For the first part in the RHS of the above inequality, we have,

$\displaystyle\langle\nabla
F(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n-1}\rangle\geq
F(\boldsymbol{\theta}_{n-1})+\frac{\lambda}{2}\|\boldsymbol{\theta}_{n-1}\|^{2}\geq\lambda\|\boldsymbol{\theta}_{n-1}\|^{2},$

with strong convexity property. Moreover,

$\displaystyle\left|\eta_{n}\mathbb{E}_{n-1}\langle\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n},\boldsymbol{\theta}_{n-1}\rangle\right|$
$\displaystyle=\eta_{n}\left|\mathbb{E}_{n-1}\langle\boldsymbol{\xi}_{n},\boldsymbol{\theta}_{n-1}\rangle\right|$
$\displaystyle\leq\frac{1}{2}\eta_{n}h_{n}L_{f}\|\boldsymbol{\theta}_{n-1}\|\mathbb{E}\|\boldsymbol{v}\|^{3}$
$\displaystyle\leq
CL_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{3}h_{n}^{2}\|\boldsymbol{\theta}_{n-1}\|^{2}+C\mathbb{E}\|\boldsymbol{v}\|^{3}\eta_{n}^{2},$
(30)
$\displaystyle\mathbb{E}_{n-1}\|\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}-\nabla
F(\boldsymbol{\theta}_{n-1})\|^{2}$ $\displaystyle\leq
4\|\boldsymbol{\xi}_{n}\|^{2}+4\|\boldsymbol{\gamma}_{n}\|^{2}+4\|\boldsymbol{\varepsilon}_{n}\|^{2}+4\|\nabla
F(\boldsymbol{\theta}_{n-1})\|^{2}$ $\displaystyle\leq
h_{n}^{2}L_{f}^{2}\mathbb{E}(\|\boldsymbol{v}\|^{3})^{2}+4L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{4}\|\boldsymbol{\theta}_{n-1}\|^{2}+h_{n}^{2}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{6}$
$\displaystyle\quad+4M^{\frac{2}{2+\delta}}\mathbb{E}\|\boldsymbol{v}_{n}\|^{4}(\|\boldsymbol{\theta}_{n-1}\|^{2}+h_{n}^{2}\|\boldsymbol{v}_{n}\|^{2}+1)+4L_{f}^{2}\|\boldsymbol{\theta}_{n-1}\|^{2}$
$\displaystyle:=M_{1}\|\boldsymbol{\theta}_{n-1}\|^{2}+M_{2}$ (31)

where we use Cauchy-Schwarz inequality in (A.1), (A.1) and
$M_{1}=C\bigl{(}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{4}+M^{\frac{2}{2+\delta}}\mathbb{E}\|\boldsymbol{v}\|^{4}+L_{f}^{2}\bigr{)},M_{2}=C\bigl{(}h_{n}^{2}L_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{6}+M^{\frac{2}{2+\delta}}\mathbb{E}\|\boldsymbol{v}\|^{4}(h_{n}^{2}\|\boldsymbol{v}\|^{2}+1)\bigr{)}$.
So combining all inequalities, we have

$\displaystyle\mathbb{E}_{n-1}\|\boldsymbol{\theta}_{n}\|^{2}\leq\left[1-2\lambda\eta_{n}+M_{1}\eta_{n}^{2}+M_{3}h_{n}^{2}\right]\|\boldsymbol{\theta}_{n-1}\|^{2}+M_{4}\eta_{n}^{2},$
(32)

where $M_{3},M_{4}$ is defined by
$M_{3}=CL_{f}^{2}\mathbb{E}\|\boldsymbol{v}\|^{3}$,
$M_{4}=C(\mathbb{E}\|\boldsymbol{v}\|^{3}+M_{2})$. Following the proof of
Theorem 1 of Moulines and Bach (2011), we can apply the recursion and get

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2}\leq\prod_{k=1}^{n}\left[1-2\lambda\eta_{k}+M_{1}\eta_{k}^{2}+CM_{3}h_{k}^{2}\right]\|\boldsymbol{\theta}_{0}\|^{2}+M_{4}\sum_{k=1}^{n}\prod_{i=k+1}^{n}\left[1-2\lambda\eta_{i}+M_{1}\eta_{k}^{2}+M_{3}h_{k}^{2}\right]\eta_{k}^{2}.$

We can then bound the first term on the RHS,

$\displaystyle\prod_{k=1}^{n}\left[1-2\lambda\eta_{k}+M_{1}\eta_{k}^{2}+M_{3}h_{k}^{2}\right]\leq\exp\left(-2\lambda\sum_{k=1}^{n}\eta_{k}\right)\exp\left(M_{1}\sum_{k=1}^{n}\eta_{k}^{2}\right)\exp\left(M_{3}\sum_{k=1}^{n}h_{k}^{2}\right),$

as well as the second term on the RHS

$\displaystyle\sum_{k=1}^{n}\prod_{i=k+1}^{n}\left[1-2\lambda\eta_{i}+M_{1}\eta_{k}^{2}+M_{3}h_{k}^{2}\right]\eta_{k}^{2}$
$\displaystyle\leq\;\;$
$\displaystyle\exp\left(-\lambda\sum_{k=m+1}^{n}\eta_{k}\right)\sum_{k=1}^{n}\eta_{k}^{2}+\frac{\eta_{m}}{\lambda}+\frac{1}{M_{1}}\exp\left(M_{1}\sum_{k=1}^{n_{0}}\eta_{k}^{2}\right)\exp\left(M_{3}\sum_{k=1}^{n_{0}}h_{k}^{2}\right)\exp\left(-\lambda\sum_{k=1}^{n}\eta_{k}\right),$

where we denote by
$n_{0}=\inf\\{k\in\mathbb{N},1-2\lambda\eta_{k}+M_{1}\eta_{k}^{2}+M_{3}h_{k}^{2}\leq
1-\lambda\eta_{k}\\}$ and $m$ is any integer in $\\{1,\dots,n\\}$. Choose
$m=n/2$ and bound $n_{0}$ by $n$. Notice that $\sum_{k=1}^{n}\eta_{k}^{2}$
converge. So we can get

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2}$
$\displaystyle\leq\exp\left(CM_{1}\eta_{0}/(2\alpha-1)+CM_{3}/(2\beta-1)-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right)\|\boldsymbol{\theta}_{0}\|^{2}$
$\displaystyle\quad+M_{4}\left(\exp\left(-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right)+\frac{\eta_{0}n^{-\alpha}}{\lambda}\right)$
$\displaystyle\quad+\frac{M_{4}}{M_{1}}\exp\left(CM_{1}\eta_{0}/(2\alpha-1)+CM_{3}/(2\beta-1)-C\lambda\eta_{0}n^{1-\alpha}/(1-\alpha)\right).$

Only the term $M_{4}\eta_{0}n^{-\alpha}/\lambda$ decreases at the order of
$O(n^{-\alpha})$ while all the other terms decrease much faster.

Notice that up to this point, all $C$’s are universal constants which do not
depend on any parameters in the assumptions. From now on, we will absorb all
parameters (other than $n$) into $C$ to make the asymptotic analysis more
clear.

By martingale convergence theorem, $\|\boldsymbol{\theta}_{n}\|$ converges
almost surely. Because its second moment converges to $\mathbf{0}$, it must
converge to $\mathbf{0}$ almost surely.

We now show that,

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}-\boldsymbol{\theta}^{\star}\|^{2+\delta}\leq
Cn^{-\alpha(2+\delta)/2}.$

By same arguments as in (26), (A.1), (A.1), we can get
$\|\boldsymbol{\xi}_{n}\|^{2+\delta}\leq Ch_{n}^{2+\delta}$,
$\mathbb{E}_{n-1}\|\boldsymbol{\gamma}_{n}\|^{2+\delta}\leq\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+Ch_{n}^{2+\delta}$,
$\mathbb{E}_{n-1}\left[\|\boldsymbol{\varepsilon}_{n}\|^{2+\delta}\right]\leq
C(\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+1)$.

By similar arguments as in Lemma B.3, there exists constants $C$ such that for
any $\boldsymbol{a},\boldsymbol{b}$,

$\displaystyle\|\boldsymbol{a}+\boldsymbol{b}\|^{2+\delta}\leq\|\boldsymbol{a}\|^{2+\delta}+(2+\delta)\langle\boldsymbol{a},\boldsymbol{b}\rangle\|\boldsymbol{a}\|^{\delta}+C\|\boldsymbol{a}\|^{\delta}\|\boldsymbol{b}\|^{2}+C\|\boldsymbol{b}\|^{2+\delta}.$

So we have the bound

$\displaystyle\mathbb{E}_{n-1}\|\boldsymbol{\theta}_{n}\|^{2+\delta}$
$\displaystyle\leq\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+\eta_{n}(2+\delta)\mathbb{E}_{n-1}\langle\boldsymbol{\theta}_{n-1},-\nabla
F(\boldsymbol{\theta}_{n-1})+\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}\rangle\|\boldsymbol{\theta}_{n-1}\|^{\delta}$
$\displaystyle\quad+C\eta_{n}^{2}\|\boldsymbol{\theta}_{n-1}\|^{\delta}\mathbb{E}_{n-1}\|-\nabla
F(\boldsymbol{\theta}_{n-1})+\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}\|^{2}$
$\displaystyle\quad+C\eta_{n}^{2+\delta}\mathbb{E}_{n-1}\|-\nabla
F(\boldsymbol{\theta}_{n-1})+\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}\|^{2+\delta}$
$\displaystyle\leq(1-(2+\delta)\lambda\eta_{n})\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+C\eta_{n}h_{n}\|\boldsymbol{\theta}_{n-1}\|^{1+\delta}$
$\displaystyle\quad+C\eta_{n}^{2}(\|\boldsymbol{\theta}_{n-1}\|^{2}+1)\|\boldsymbol{\theta}_{n-1}\|^{\delta}+C\eta_{n}^{2+\delta}(\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+1).$

If $0<\delta\leq 1$, by previous bound
$\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2}\leq Cn^{-\alpha}$, we can get
$\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{1+\delta}\leq
Cn^{-\alpha(1+\delta)/2}$ and
$\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{\delta}\leq Cn^{-\alpha\delta/2}$ by
Hölder’s inequality. So we can further get

$\displaystyle\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2+\delta}\leq(1-Cn^{-\alpha}+Cn^{-2\alpha})\mathbb{E}\|\boldsymbol{\theta}_{n-1}\|^{2+\delta}+Cn^{-(2+\delta)\alpha/2},$

which implies $\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2+\delta}\leq
Cn^{-(2+\delta)\alpha/2}$ as in the above proof after (32).

Now the case for $0<\delta\leq 1$ is proved. We can then use induction. If
$\mathbb{E}\|\boldsymbol{\theta}_{n}\|^{2+\delta}\leq
Cn^{-(2+\delta)\alpha/2}$ for all $\delta\leq n$, then we can use the same
method to prove the same inequality holds for $\delta\in(n,n+1]$. Thus the
inequality holds for all $\delta$. ∎

#### Proof of Lemma 3.2

###### Proof.

By Assumption 2, we know that

$\displaystyle\mathbb{E}\|\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta})\|^{2+\delta}\leq
M(\|\boldsymbol{\theta}\|^{2+\delta}+d^{2+\delta}).$

Therefore, the following holds some constant $C>0$,

$\displaystyle\mathbb{E}\|\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta})\|^{2}\leq C(\|\boldsymbol{\theta}\|^{2}+d^{2}).$ (33)

In particular,

$\displaystyle\mathbb{E}\|\nabla f(\mathbf{0};\boldsymbol{\zeta})-\nabla
F(\mathbf{0})\|^{2}\leq C.$ (34)

From Assumption 3, we can get the following estimate for the Hessian matrix
$\nabla^{2}f(\boldsymbol{\theta};\boldsymbol{\zeta})$,

$\displaystyle\mathbb{E}\|\nabla^{2}f(\boldsymbol{\theta};\boldsymbol{\zeta})\|^{2}$
$\displaystyle\leq
2\mathbb{E}\|\nabla^{2}f(\mathbf{0};\boldsymbol{\zeta})\|^{2}+2\mathbb{E}\left\|\nabla^{2}f(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla^{2}f(\mathbf{0};\boldsymbol{\zeta})\right\|^{2}$
$\displaystyle\leq C(1+\|\boldsymbol{\theta}\|^{2}).$

Using the above observation, we find that

$\displaystyle\quad\mathbb{E}\|\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla F(\boldsymbol{\theta})-\nabla
f(\mathbf{0};\boldsymbol{\zeta})+\nabla F(\mathbf{0})\|^{2}$
$\displaystyle\leq C\|\boldsymbol{\theta}\|^{2}+2\mathbb{E}\|\nabla
f(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
f(\mathbf{0};\boldsymbol{\zeta}))\|^{2}$
$\displaystyle=C\|\boldsymbol{\theta}\|^{2}+2\mathbb{E}\left\|\int_{0}^{1}\nabla^{2}f(s\boldsymbol{\theta};\boldsymbol{\zeta})\boldsymbol{\theta}\,\mathrm{d}s\right\|^{2}$
$\displaystyle\leq
C\|\boldsymbol{\theta}\|^{2}+2\mathbb{E}\int_{0}^{1}\|\nabla^{2}f(s\boldsymbol{\theta};\boldsymbol{\zeta})\boldsymbol{\theta}\|^{2}\,\mathrm{d}s$
$\displaystyle\leq
C\|\boldsymbol{\theta}\|^{2}(1+\int_{0}^{1}\mathbb{E}\|\nabla^{2}f(s\boldsymbol{\theta};\boldsymbol{\zeta})\|^{2}\,\mathrm{d}s)$
$\displaystyle\leq
C\|\boldsymbol{\theta}\|^{2}(1+\|\boldsymbol{\theta}\|^{2}).$ (35)

Define the function $\Sigma(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2})$
by

$\displaystyle\Sigma(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}):=\mathbb{E}(\nabla
f(\boldsymbol{\theta}_{1};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{1}))(\nabla
f(\boldsymbol{\theta}_{2};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{2}))^{\top}.$

Then combining inequalities (33), (34), (A.1), we have

$\displaystyle\|\Sigma(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2})-S\|$
$\displaystyle\leq\mathbb{E}\|(\nabla
f(\boldsymbol{\theta}_{1};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{1}))(\nabla
f(\boldsymbol{\theta}_{2};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{2}))^{\top}$ $\displaystyle\quad-(\nabla
f(\mathbf{0};\boldsymbol{\zeta})-\nabla F(\mathbf{0}))(\nabla
f(\mathbf{0};\boldsymbol{\zeta})-\nabla F(\mathbf{0}))^{\top}\|$
$\displaystyle\leq\mathbb{E}\|\nabla
f(\boldsymbol{\theta}_{1};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{1})\|\|\nabla
f(\boldsymbol{\theta}_{2};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{2})-\nabla f(\mathbf{0};\boldsymbol{\zeta})+\nabla
F(\mathbf{0})\|$ $\displaystyle\quad+\mathbb{E}\|\nabla
f(\boldsymbol{\theta}_{1};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta}_{2})-\nabla f(\mathbf{0};\boldsymbol{\zeta})+\nabla
F(\mathbf{0})\|\|\nabla f(\mathbf{0};\boldsymbol{\zeta})-\nabla
F(\mathbf{0})\|$ $\displaystyle\leq
C(d+\|\boldsymbol{\theta}_{1}\|)\|\boldsymbol{\theta}_{2}\|(1+\|\boldsymbol{\theta}_{2}\|)+C\|\boldsymbol{\theta}_{1}\|(1+\|\boldsymbol{\theta}_{1}\|).$
(36)

Notice that

$\displaystyle\quad\mathbb{E}_{\boldsymbol{\zeta}}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})^{\top}-(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})^{\top}$
$\displaystyle=\mathbb{E}_{\boldsymbol{\zeta}}(\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})-\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})(\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})-\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})^{\top}$
$\displaystyle=\frac{1}{h^{2}}\mathbb{E}_{\boldsymbol{\zeta}}\boldsymbol{v}(f(\boldsymbol{\theta}+h\boldsymbol{v};\boldsymbol{\zeta})-f(\boldsymbol{\theta};\boldsymbol{\zeta})-F(\boldsymbol{\theta}+h\boldsymbol{v})+F(\boldsymbol{\theta}))^{2}\boldsymbol{v}^{\top}$
$\displaystyle=\frac{1}{h^{2}}\mathbb{E}_{\boldsymbol{\zeta}}\boldsymbol{v}\boldsymbol{v}^{\top}\bigg{[}\int_{0}^{h}\int_{0}^{h}\left(\nabla
F(\boldsymbol{\theta}+s_{1}\boldsymbol{v})-\nabla
f(\boldsymbol{\theta}+s_{1}\boldsymbol{v};\boldsymbol{\zeta})\right)$
$\displaystyle\quad\left(\nabla
F(\boldsymbol{\theta}+s_{2}\boldsymbol{v})-\nabla
f(\boldsymbol{\theta}+s_{2}\boldsymbol{v};\boldsymbol{\zeta})\right)^{\top}\,\mathrm{d}s_{1}\,\mathrm{d}s_{2}\bigg{|}\bigg{]}\boldsymbol{v}\boldsymbol{v}^{\top}$
$\displaystyle=\frac{1}{h^{2}}\mathbb{E}_{\boldsymbol{\zeta}}\boldsymbol{v}\boldsymbol{v}^{\top}\int_{0}^{h}\int_{0}^{h}\Sigma(\boldsymbol{\theta}+s_{1}\boldsymbol{v},\boldsymbol{\theta}+s_{2}\boldsymbol{v})\,\mathrm{d}s_{1}\,\mathrm{d}s_{2}\boldsymbol{v}\boldsymbol{v}^{\top}.$

We can use (A.1) and derive that

$\displaystyle\quad\|\mathbb{E}_{\boldsymbol{\zeta}}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})^{\top}-(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})^{\top}-\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\|$
$\displaystyle\leq
C\|\boldsymbol{v}\|^{4}(\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(1+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(d+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|).$

Now we have

$\displaystyle\quad\|\mathbb{E}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})^{\top}-\mathbb{E}(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})^{\top}-\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\|$
$\displaystyle\leq
C\mathbb{E}\|\boldsymbol{v}\|^{4}(\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(1+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(d+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|).$
(37)

By the same argument,

$\displaystyle\quad\|\mathbb{E}(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})(\frac{1}{h}\Delta_{h,\boldsymbol{v}}F(\boldsymbol{\theta})\boldsymbol{v})^{\top}\|$
$\displaystyle\leq\frac{1}{h^{2}}\mathbb{E}\bigg{\|}\boldsymbol{v}\boldsymbol{v}^{\top}\bigg{[}\int_{0}^{h}\int_{0}^{h}\left(\nabla
F(\boldsymbol{\theta}+s_{1}\boldsymbol{v})\right)\left(\nabla
F(\boldsymbol{\theta}+s_{2}\boldsymbol{v})\right)^{\top}\,\mathrm{d}s_{1}\,\mathrm{d}s_{2}\bigg{|}\bigg{]}\boldsymbol{v}\boldsymbol{v}^{\top}\bigg{\|}$
$\displaystyle\leq
C\mathbb{E}\|\boldsymbol{v}\|^{4}(\|\boldsymbol{\theta}\|^{2}+h^{2}\|\boldsymbol{v}\|^{2}).$

So we finally get

$\displaystyle\|\mathbb{E}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})^{\top}-\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\|\leq
C\mathbb{E}\|\boldsymbol{v}\|^{4}(\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(1+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|)(d+\|\boldsymbol{\theta}\|+h\|\boldsymbol{v}\|).$

for some constant $C>0$. ∎

#### Proof of Theorem 3.3

###### Proof.

We follow the proof in Polyak and Juditsky (1992). The update step is

$\displaystyle\boldsymbol{\theta}_{n}$
$\displaystyle=\boldsymbol{\theta}_{n-1}-\eta_{n}\nabla
F(\boldsymbol{\theta}_{n-1})+\eta_{n}(\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n})$
$\displaystyle=(I_{d}-\eta_{n}H)\boldsymbol{\theta}_{n-1}+\eta_{n}(H\boldsymbol{\theta}_{n-1}-\nabla
F(\boldsymbol{\theta}_{n-1})+\boldsymbol{\xi}_{n}+\boldsymbol{\gamma}_{n}+\boldsymbol{\varepsilon}_{n}).$

By the argument in Polyak and Juditsky (1992), we only need to prove the
following three conditions. First,

$\displaystyle\sum_{i=1}^{\infty}\frac{1}{\sqrt{i}}\mathbb{E}\|H\boldsymbol{\theta}_{i-1}-\nabla
F(\boldsymbol{\theta}_{i-1})+\boldsymbol{\xi}_{i}\|,$ (38)

is bounded almost surely. Furthermore, we have

$\displaystyle\mathbb{E}\|\boldsymbol{\gamma}_{i}+\boldsymbol{\varepsilon}_{i}\|^{2},$
(39)

is bounded almost surely. Moreover, when $t\rightarrow\infty$, we have the
following convergence in probability,

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}(\boldsymbol{\gamma}_{i}+\boldsymbol{\varepsilon}_{i})\Longrightarrow\mathcal{N}(\mathbf{0},Q).$
(40)

By Assumption 3, we know that

$\displaystyle\|\nabla^{2}F(\boldsymbol{\theta})-\nabla^{2}F(\boldsymbol{\theta^{\prime}})\|^{2}\leq
L_{g}\|\boldsymbol{\theta}-\boldsymbol{\theta^{\prime}}\|^{2}.$

By Taylor expansion,

$\displaystyle\|H\boldsymbol{\theta}_{i-1}-\nabla
F(\boldsymbol{\theta}_{i-1})\|\leq C\|\boldsymbol{\theta}_{i-1}\|^{2}.$

By facts (26) to (A.1), we know that

$\displaystyle\mathbb{E}\|H\boldsymbol{\theta}_{i-1}-\nabla
F(\boldsymbol{\theta}_{i-1})+\boldsymbol{\xi}_{i}\|\leq Ci^{-\alpha},$

which indicates that condition (38) holds.

Because $\boldsymbol{\gamma}_{i}$ converges to $\mathbf{0}$ almost surely and
$\boldsymbol{\varepsilon}_{i}$ has bounded variance. So condition (39) holds.
To prove condition (40), it suffices to verify that,

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\boldsymbol{\varepsilon}_{i}\Longrightarrow\mathcal{N}(\mathbf{0},Q).$

By martingale central limit theorem (Durrett, 2019, Theorem 8.2.8), we only
need to verify two conditions,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{i-1}[\boldsymbol{\varepsilon}_{i}\boldsymbol{\varepsilon}_{i}^{\top}]\rightarrow
Q,$ (41)
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\|\boldsymbol{\varepsilon}_{i}\|^{2}\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]\rightarrow
0,$ (42)

in probability for all $a>0$.

Notice that (A.1) is equivalent to the following inequality,

$\displaystyle\|\mathbb{E}_{n-1}\boldsymbol{\varepsilon}_{n}\boldsymbol{\varepsilon}_{n}^{\top}-\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}\|\leq
C(\|\boldsymbol{\theta}_{n-1}\|+h_{n})(1+\|\boldsymbol{\theta}_{n-1}\|^{3}+h_{n}^{3}).$
(43)

Thus
$\mathbb{E}_{n-1}[\boldsymbol{\varepsilon}_{n}\boldsymbol{\varepsilon}_{n}^{\top}]$
converges almost surely to $Q$ and condition (41) holds.

Now consider the quantity in (42), by Proposition 3.1,

$\displaystyle\mathbb{E}_{i-1}\left[\|\boldsymbol{\varepsilon}_{i}\|^{2}\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]\leq\left[\mathbb{E}_{i-1}\left[\|\boldsymbol{\varepsilon}_{i}\|^{2+\delta}\right]\right]^{\frac{2}{2+\delta}}\left[\mathbb{E}_{i-1}\left[\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]\right]^{\frac{\delta}{2+\delta}}.$

Note that

$\displaystyle\mathbb{E}_{i-1}\left[\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]=\mathbb{P}_{i-1}\left(\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}|\boldsymbol{\theta}_{i-1}\right)\leq\frac{1}{a\sqrt{n}}\mathbb{E}_{i-1}\|\boldsymbol{\varepsilon}_{i}\|.$

Therefore, it can be bounded by

$\displaystyle\mathbb{E}_{i-1}\left[\|\boldsymbol{\varepsilon}_{i}\|^{2}\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]\leq
C\left(\frac{1}{a\sqrt{n}}\right)^{\frac{\delta}{2+\delta}}\left(1+\|\boldsymbol{\theta}_{i-1}\|^{2+\delta}\right)^{\frac{2}{2+\delta}}\left(1+\|\boldsymbol{\theta}_{i-1}\|\right)^{\frac{\delta}{2+\delta}},$

from which we can obtain that

$\displaystyle\mathbb{E}[\|\boldsymbol{\varepsilon}_{i}\|^{2}\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}]\leq
C\left(\frac{1}{a\sqrt{n}}\right)^{\frac{\delta}{2+\delta}}.$ (44)

We find that condition (42) holds when $n$ goes to infinity:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}\left[\|\boldsymbol{\varepsilon}_{i}\|^{2}\mathbf{1}_{\|\boldsymbol{\varepsilon}_{i}\|>a\sqrt{n}}\right]\leq
C\left(\frac{1}{a\sqrt{n}}\right)^{\frac{\delta}{2+\delta}}\rightarrow 0.$

Therefore, we conclude the result. ∎

#### Proof of Proposition 3.5

###### Proof.

For $Q^{\tt{(G)}}$, let $\boldsymbol{z}\sim\mathcal{N}(\mathbf{0},I_{d})$, and
we now calculate
$\mathbb{E}\boldsymbol{z}\boldsymbol{z}^{\top}S\boldsymbol{z}\boldsymbol{z}^{\top}$.
The $(i,i)$-th entry is

$\displaystyle\mathbb{E}\sum_{j,k}z_{i}z_{j}S_{jk}z_{k}z_{i}=\sum_{j\neq
i}S_{jj}+3S_{ii}=2S_{ii}+\mathrm{tr}(S).$

For $i\neq j$, the $(i,j)$th entry is

$\displaystyle\mathbb{E}\sum_{k,l}z_{i}z_{k}S_{kl}z_{l}z_{j}=2S_{ij}.$

So
$\mathbb{E}\boldsymbol{z}\boldsymbol{z}^{\top}S\boldsymbol{z}\boldsymbol{z}^{\top}=2S+\mathrm{tr}(S)I_{d}$.

For $Q^{\tt{(S)}}$, let $\boldsymbol{v}$ be sampled from the uniform
distribution on the sphere $\|\boldsymbol{v}\|=d$. The Gaussian vector
$\boldsymbol{z}$ can be decomposed into independent radius part and spherical
part,

$\displaystyle\mathbb{E}[\boldsymbol{z}\boldsymbol{z}^{\top}]=\mathbb{E}\left[\|\boldsymbol{z}\|^{2}\frac{\boldsymbol{z}}{\|\boldsymbol{z}\|}\frac{\boldsymbol{z}^{\top}}{\|\boldsymbol{z}\|}\right]$
$\displaystyle=\ \mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top},$
$\displaystyle\mathbb{E}[\boldsymbol{z}\boldsymbol{z}^{\top}S\boldsymbol{z}\boldsymbol{z}^{\top}]=\mathbb{E}\left[\|\boldsymbol{z}\|^{4}\frac{\boldsymbol{z}}{\|\boldsymbol{z}\|}\frac{\boldsymbol{z}^{\top}}{\|\boldsymbol{z}\|}S\frac{\boldsymbol{z}}{\|\boldsymbol{z}\|}\frac{\boldsymbol{z}^{\top}}{\|\boldsymbol{z}\|}\right]$
$\displaystyle=\frac{d+2}{d}\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}.$

Now we have

$\displaystyle\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}=I_{d},\;\;\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}=\frac{d}{d+2}(2S+\mathrm{tr}(S)I_{d}).$

For $Q^{\tt{(U)}}$, let $\boldsymbol{u}$ obey the uniform distribution on
$\\{\sqrt{d}e_{1},\dots,\sqrt{d}e_{d}\\}$. By direct calculation, we have

$\displaystyle\mathbb{E}\boldsymbol{u}\boldsymbol{u}^{\top}S\boldsymbol{u}\boldsymbol{u}^{\top}=\sum_{j=1}^{d}\frac{1}{d}\cdot
d^{2}S_{jj}=d\ \mathrm{diag}(S).$

The final two cases for $Q^{\tt{(U)}},Q^{\tt{(P)}}$ can also be verified by
direct calculation. ∎

###  A.2 Extensions to local strong convexity and nonsmoothness

#### Asymptotic behavior for locally strongly convex loss function

To comply with the settings of the logistic regression, we need to consider a
relaxed version of Assumption 1 as follows,

######  Assumption 1′.

The population loss function $F(\boldsymbol{\theta})$ is twice continuously
differentiable, convex and $L_{f}$-smooth. In addition, there exists
$\delta_{1}>0$ such that for all $\boldsymbol{\theta}$ in the
$\delta_{1}$-ball centered at $\boldsymbol{\theta}^{\star}$, the Hessian
matrix $\nabla^{2}F(\boldsymbol{\theta})$ is positive-definite.

Assumption 1′ considers local strong convexity of the population objective
$F(\cdot)$ at the minimizer $\boldsymbol{\theta}^{\star}$. Intuitively, after
a number of steps in the (KW) SGD update, the estimated parameter
$\boldsymbol{\theta}_{n}$ would be sufficiently close to
$\boldsymbol{\theta}^{\star}$ and we have the strong convexity from there.
This assumption naturally suites the settings of the logistic regression.

######  Theorem A.2.

Let Assumption 1′, and 2 to 4 hold. Set the step size as
$\eta_{n}=\eta_{0}n^{-\alpha}$ for some constant $\eta_{0}>0$ and
$\alpha\in\left(\frac{1}{2},1\right)$, and the spacing parameter as
$h_{n}=h_{0}n^{-\gamma}$ for some constant $h_{0}>0$, and
$\gamma\in\left(\frac{1}{2},1\right)$. The averaged estimator
$\overline{\boldsymbol{\theta}}_{n}$ satisfies,

$\displaystyle\sqrt{n}\;\left(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}QH^{-1}\right),\qquad\text{as}\quad
n\rightarrow\infty.$

###### Proof.

Under assumption 1′, using Lemma B.1 in Su and Zhu (2018), for all
$\boldsymbol{\theta}$ in the $\delta_{1}$-ball centered at $\mathbf{0}$, we
have

$\displaystyle\langle\boldsymbol{\theta},\nabla
F(\boldsymbol{\theta})\rangle\geq\rho\|\boldsymbol{\theta}\|\min\left\\{\|\boldsymbol{\theta}\|,\delta_{1}\right\\}.$
(45)

for some $\rho>0$. For the first part in the RHS of the previous inequality
(A.1), using inequality (45), we have

$\displaystyle\langle\nabla
F(\boldsymbol{\theta}_{n-1}),\boldsymbol{\theta}_{n-1}\rangle\geq\rho\|\boldsymbol{\theta}_{n-1}\|\min\left\\{\|\boldsymbol{\theta}_{n-1}\|,\delta_{1}\right\\},$

Furthermore, by (A.1) and (A.1), we have

$\displaystyle\mathbb{E}_{n-1}\|\boldsymbol{\theta}_{n}\|^{2}\leq\left(1+C_{2}^{2}\eta_{n}^{2}\right)\|\boldsymbol{\theta}_{n-1}\|^{2}-2\eta_{n}\rho\|\boldsymbol{\theta}_{n-1}\|\min\left\\{\|\boldsymbol{\theta}_{n-1}\|,\delta_{1}\right\\}+C_{2}\eta_{n}^{2}.$

By Robbins-Siegmund theorem, $\|\boldsymbol{\theta}_{n}\|^{2}$ converges to
some random variable almost surely and

$\displaystyle\sum_{n=1}^{\infty}2\eta_{n}\rho\|\boldsymbol{\theta}_{n}\|\min\left\\{\|\boldsymbol{\theta}_{n}\|,\delta_{1}\right\\}<\infty.$

Combining the fact that $\sum_{n=1}^{\infty}\eta_{n}=\infty$ we can yield that
$\boldsymbol{\theta}_{n}$ converges almost surely to $\mathbf{0}$. The rest
part follows from the proof of Theorem 3.3. ∎

#### Asymptotic behavior of (AKW) estimator for nonsmooth loss functions

We present Theorem 3.3 in the main paper with strengthened Assumptions 2 and
3, where we assume the existence of the gradient of the inaccessible (RM)
stochastic gradient $g(\boldsymbol{\theta};\boldsymbol{\zeta})$ and its Gram
matrix
$S=\mathbb{E}[g(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})g(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}]$.
The theoretical analysis of the asymptotic distribution of the (AKW) estimator
remains working with a weakened assumption, which is a natural fit to some
nonsmooth loss functions $F(\boldsymbol{\theta})$ including the quantile
regression in Example 2.3.

######  Assumption 5.

Assume there exists $C_{1}>0$ such that
$\big{\|}\mathbb{E}\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta})-\nabla
F(\boldsymbol{\theta})\big{\|}\leq C_{1}h$ for any $h>0$ and
$\boldsymbol{\theta}\in\mathbb{R}^{d}$. Further assume there exists $C_{2}>0$
such that
$\mathbb{E}\left[\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta}^{\star};\boldsymbol{\zeta})^{\top}\right]=Q+\Delta_{h}$
for some matrix $Q\in\mathbb{R}^{d\times d}$ and $\|\Delta_{h}\|\leq
C_{2}h^{\iota}$, for some $\iota>0$. Moreover, for some $0<\delta\leq 2$,
there exists $M>0$,

$\displaystyle\mathbb{E}\|\widehat{g}_{h,\boldsymbol{v}}(\boldsymbol{\theta};\boldsymbol{\zeta}_{n})-\nabla
F(\boldsymbol{\theta})\|^{2+\delta}\leq
M(\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\star}\|^{2+\delta}+h^{2+\delta}).$

######  Theorem A.3.

Let Assumption 1 and 5 hold. Under the step size and spacing parameter
conditions specified in Theorem 3.3, the averaged estimator
$\overline{\boldsymbol{\theta}}_{n}$ satisfies,

$\displaystyle\sqrt{n}\;\left(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star}\right)\Longrightarrow\mathcal{N}\left(\mathbf{0},H^{-1}QH^{-1}\right),\qquad\text{as}\quad
n\rightarrow\infty.$ (46)

###### Proof.

Under Assumption 1 and 5, the conclusions in Lemma 2.4 and Lemma 3.2 naturally
hold. The rest of the proof follows from the proof in Proposition 3.1 and
Theorem 3.3. ∎

###  A.3 Illustration of choices of directions $\mathcal{P}_{\boldsymbol{v}}$

We first note that $Q^{\tt{(G)}}\succ Q^{\tt{(S)}}$ regardless of the
dimension $d$ and Gram matrix $S$. Intuitively, when the direction
$\boldsymbol{v}$ is generated by Gaussian (G), it can be decomposed into two
independent random variables: the radical part $\|\boldsymbol{v}\|$ and the
spherical part $\boldsymbol{v}/\|\boldsymbol{v}\|$. The spherical part
$\boldsymbol{v}/\|\boldsymbol{v}\|$ follows the same distribution as the
uniform distribution on the sphere with radius $d$ (which is identical to
(S)). The extra randomness in the radical part
$\|\boldsymbol{v}\|^{2}\sim\chi^{2}(d)$ leads to a larger magnitude of $Q$
compared to that of (S). Therefore the (AKW) estimator with Gaussian
directions (G) is always inferior to that with spherical directions (S),
asymptotically. However, for the other candidates, they are not directly
comparable, and the optimal choice of $\mathcal{P}_{\boldsymbol{v}}$ depends
on the optimality criterion, and Gram matrix $S$.

As a simple illustration, we consider $S=\mathrm{diag}(1,r_{0})$ for some
$r_{0}>0$. We have

  1. (S)

Spherical:
$Q^{\tt{(S)}}=\mathrm{diag}\left(\frac{r_{0}+3}{2},\frac{3r_{0}+1}{2}\right)$.

  2. (I)

Uniform in a natural coordinate basis: $Q^{\tt{(I)}}=\mathrm{diag}(2,2r_{0})$.

  3. (U)

Uniform in an arbitrary orthonormal basis $U$: when
$U=(\cos\omega,\sin\omega;-\sin\omega,\cos\omega)$ and $\omega=0$, we have
$Q^{\tt{(U)}}=Q^{\tt{(I)}}=\mathrm{diag}(2,2r_{0})$; when $\omega=\pi/4$, we
have $Q^{\tt{(U)}}=\mathrm{diag}(1+r_{0},1+r_{0})$.

  4. (P)

Non-uniform in a natural coordinate basis:
$\mathrm{diag}\left(\frac{1}{p_{1}},\frac{r_{0}}{1-p_{1}}\right),p_{1}\in(0,1)$.

From the above we can see that, the choices of the distribution of direction
vectors $\mathcal{P}_{\boldsymbol{v}}$ depends on the optimality-criteria on
comparing the covariance matrices. Specifically in the above example, if one
seeks to minimize

  * •

the trace of covariance matrix, we have

$\displaystyle\mathrm{tr}(Q^{\tt{(S)}})=\mathrm{tr}(Q^{\tt{(I)}})=\mathrm{tr}(Q^{\tt{(U)}})=2+2r_{0},\;\;\mathrm{tr}(Q^{\tt{(P)}})=\frac{1}{p_{1}}+\frac{r_{0}}{1-p_{1}},$

and the optimal distribution that minimizes the trace depends on the value of
$p_{1}$.

  * •

the determinant of covariance matrix, we have

$\displaystyle\det(Q^{\tt{(S)}})=\frac{3r_{0}^{2}+10r_{0}+3}{4},$
$\displaystyle\det(Q^{\tt{(I)}})=4r_{0},$
$\displaystyle\det(Q^{\tt{(U)}})=\frac{-\cos(4\omega)(r_{0}-1)^{2}+r_{0}^{2}+6r_{0}+1}{2},\;\;$
$\displaystyle\det(Q^{\tt{(P)}})=\frac{r_{0}}{p_{1}(1-p_{1})}.$

By a simple derivation, we have
$\det(Q^{\tt{(S)}})\geq\det(Q^{\tt{(U)}})\geq\det(Q^{\tt{(I)}})$ and
$\det(Q^{\tt{(P)}})\geq\det(Q^{\tt{(I)}})$.

  * •

the operator norm of covariance matrix, i.e., the largest eigenvalue, we have

$\displaystyle\lambda_{\max}(Q^{\tt{(S)}})=\frac{r_{0}+3}{2},$
$\displaystyle\lambda_{\max}(Q^{\tt{(I)}})=2,$
$\displaystyle\lambda_{\max}(Q^{\tt{(P)}})=\max\left\\{\frac{1}{p_{1}},\frac{r_{0}}{1-p_{1}}\right\\},\;\;$
$\displaystyle\lambda_{\max}(Q^{\tt{(U)}})=r_{0}+1+(1-r_{0})\left|\cos(2\omega)\right|.$

The smallest operator norm for $Q^{\tt{(P)}}$ is given by
$p_{1}=\frac{1}{1+r_{0}}$. When $r_{0}\leq 1$, and $0\leq\omega\leq\pi/6$, we
have
$\lambda_{\max}(Q^{\tt{(I)}})\geq\lambda_{\max}(Q^{\tt{(U)}})\geq\lambda_{\max}(Q^{\tt{(S)}})\geq\lambda_{\max}(Q^{\tt{(P)}})$.
When $r_{0}\geq 1$, and $0\leq\omega\leq\pi/6$, we have
$\lambda_{\max}(Q^{\tt{(P)}})\geq\lambda_{\max}(Q^{\tt{(S)}})\geq\lambda_{\max}(Q^{\tt{(U)}})\geq\lambda_{\max}(Q^{\tt{(I)}})$.
For other choices of $\omega$, we can obtain a comparison analogously.

Figure A.1: Comparison of $Q$ matrices under different direction distributions
$\mathcal{P}_{\boldsymbol{v}}$ when $S=\mathrm{diag}(1,1/2)$.

In general, it is natural to use Loewner order to compare two positive semi-
definite matrix $A,B\in\mathbb{R}^{d\times d}$, i.e., $A\succeq B$ if
$\boldsymbol{x}^{\top}A\boldsymbol{x}\geq\boldsymbol{x}^{\top}B\boldsymbol{x}$
for any $\boldsymbol{x}\in\mathbb{R}^{d}$. It is equivalent to say, for any
positive constant $c>0$, the ellipsoid
$\\{\boldsymbol{x}\in\mathbb{R}^{d}:\boldsymbol{x}^{\top}A\boldsymbol{x}\leq
c\\}$ contains the ellipsoid
$\\{\boldsymbol{x}\in\mathbb{R}^{d}:\boldsymbol{x}^{\top}B\boldsymbol{x}\leq
c\\}$. To better illustrate the result, we consider the 2-dimensional case
where $S=\mathrm{diag}(1,1/2)$ and plot the ellipse
$\\{\boldsymbol{x}\in\mathbb{R}^{2}:\boldsymbol{x}^{\top}Q^{(\cdot)}\boldsymbol{x}=2\\}$.
In Figure A.1, we compare $Q^{\tt{(S)}}$, $Q^{\tt{(I)}}$ (as a special case of
$Q^{\tt{(U)}}$ with $\theta=0$), $Q^{\tt{(U)}}$ with $\theta=\frac{\pi}{6}$,
and $Q^{\tt{(P)}}$ with $p_{1}=\frac{1}{1+r_{0}}=\frac{2}{3}$. As can be
inferred from the plot, none of the ellipsoids contain any other ellipsoids.

As shown in this illustrative example, there is no unique optimal direction
distribution, and a practitioner might choose a search direction based on her
favorable optimality criterion.

Lastly, in the following Remark A.4, we show that, if the optimality criterion
degenerates to one dimension, one may utilize the non-uniform distribution (P)
to obtain a smaller limiting variance. In particular, consider the application
where we are only interested in the first coordinate of
$\boldsymbol{\theta}^{\star}$, in which cases the optimality criterion of the
limiting variance is on $\theta^{\star}_{1}$. We will show that the (AKW)
estimator with the non-uniform distribution $\tt{(P)}$ achieves the Cramér-Rao
lower bound.

######  Remark A.4.

Assume the population loss function $F(\cdot)$ has Hessian $H=I_{d}$.
Considering a non-uniform sampling $\tt{(P)}$ from
$\\{\boldsymbol{e}_{k}\\}_{k=1}^{d}$ for the direction distribution
$\mathcal{P}_{\boldsymbol{v}}$. We choose $\boldsymbol{v}=\boldsymbol{e}_{k}$
with probability $p_{k}$ for $k=1,2,\dots,d$, where $p_{1}=1-p$ for some
constant $p\in(0,1]$ and $p_{k}=p/(d-1)$ for $k\neq 1$. Define i.i.d. random
variables $k_{n}$ where $k_{n}=1$ with probability $1-p$ and
$k_{n}=2,\ldots,d$ uniformly with probability $p/(d-1)$. The gradient
estimator is defined by,

$\displaystyle\widehat{g}(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})=\frac{f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{e}_{k_{n}};\boldsymbol{\zeta}_{n})-f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})}{h_{n}p_{n}}\boldsymbol{e}_{k_{n}},$

where $p_{n}=1-p$ if $k_{n}=1$, $p_{n}=p/(d-1)$ for $k_{n}>1$. By the same
argument as Proposition 3.5, the variance for
$\overline{\boldsymbol{\theta}}_{n}$ in the direction $\boldsymbol{e}_{1}$ is,

$\displaystyle
n\mathrm{Var}\left(\boldsymbol{e}_{1}^{\top}(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star})\right)=\frac{S_{11}}{1-p}.$

As $p\rightarrow 0$, we approximately obtain the optimal variance given by
Cramér-Rao lower bound in the direction $\boldsymbol{e}_{1}$. However, in
order to approach the optimal variance in the direction $\boldsymbol{e}_{1}$,
we increase the magnitude of variance in all other directions, where the
variance in other directions is given by
$n\mathrm{Var}\left(\boldsymbol{e}_{k}^{\top}(\overline{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}^{\star})\right)=(d-1)S_{kk}/p$
for $k=2,\ldots,d$.

###  A.4 Multi-query approximation

#### Proof of Theorem 3.6

###### Proof.

The convergence result can be obtained as in the two function evaluation case.
The only difference is the following calculation:

$\displaystyle\mathbb{E}\left(\frac{1}{m}\sum_{i=1}^{m}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}\right)S\left(\frac{1}{m}\sum_{i=1}^{m}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}\right)=\frac{1}{m}\mathbb{E}\boldsymbol{v}\boldsymbol{v}^{\top}S\boldsymbol{v}\boldsymbol{v}^{\top}+\frac{m-1}{m}S,$

which implies the desired result. ∎

#### Proof of Theorem 3.7

###### Proof.

It is clear that $Q_{m}=S$ for $m=d$. We need to compute the quantity

$\displaystyle
Q_{m}=\frac{d^{2}}{m^{2}}\mathbb{E}\bigg{(}\sum_{i=1}^{m}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}\bigg{)}S\bigg{(}\sum_{i=1}^{m}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}\bigg{)},$

which can be simplifies to

$Q_{m}=\frac{d^{2}}{m^{2}}\mathbb{E}\bigg{(}\sum_{i=1}^{m}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}S\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}\bigg{)}+\frac{d^{2}}{m^{2}}\mathbb{E}\bigg{(}\sum_{i\neq
j}\boldsymbol{v}_{i}\boldsymbol{v}_{i}^{\top}S\boldsymbol{v}_{j}\boldsymbol{v}_{j}^{\top}\bigg{)}.$

By symmetry, it equals to

$Q_{m}=\frac{d^{2}}{m}\mathbb{E}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}S\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}+\frac{d^{2}(m-1)}{m}\mathbb{E}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}S\boldsymbol{v}_{2}\boldsymbol{v}_{2}^{\top}.$

We know
$\mathbb{E}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}S\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}=\frac{1}{d^{2}}Q$
and $Q_{d}=S$. So we can solve for
$\mathbb{E}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}S\boldsymbol{v}_{2}\boldsymbol{v}_{2}^{\top}$
and get

$\mathbb{E}\boldsymbol{v}_{1}\boldsymbol{v}_{1}^{\top}S\boldsymbol{v}_{2}\boldsymbol{v}_{2}^{\top}=\frac{1}{d(d-1)}(\frac{1}{d}Q-\mathrm{diag}(S)).$

Therefore,

$\displaystyle Q_{m}$
$\displaystyle=\frac{1}{m}Q+\frac{d(m-1)}{m(d-1)}(\frac{1}{d}Q-\mathrm{diag}S)$
$\displaystyle=\frac{d-m}{m(d-1)}Q+\frac{d(m-1)}{m(d-1)}S.\qed$

##  Appendix B Proofs of Results in Section 4

###  B.1 Proof of Lemma 4.1

Before we come to the proof of the Hessian estimator (17) in Lemma 4.1, we
first introduce a naive method to estimate Hessian matrix $H$ which we omit in
the main text.

Inspired by the previous gradient estimator, we can estimate the Hessian
matrix $H$ by the following

$\displaystyle\widehat{G}_{n}=\frac{1}{mh_{n}^{2}}\sum_{j=1}^{m}\left[\Delta_{h_{n}\boldsymbol{v}_{n}^{(j)}}f(\boldsymbol{\theta}_{n-1}+h_{n}\boldsymbol{u}_{n}^{(j)};\boldsymbol{\zeta}_{n})-\Delta_{h_{n}\boldsymbol{v}_{n}^{(j)}}f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})\right]\boldsymbol{u}_{n}^{(j)}\boldsymbol{v}_{n}^{(j)\top},$

where $\\{\boldsymbol{u}_{n}^{(j)}\\}_{j=1}^{m}$ and
$\\{\boldsymbol{v}_{n}^{(j)}\\}_{j=1}^{m}$ are _i.i.d._ random vectors and
$m>0$ is a parameter (which might be different from $m$ in the previous
section). Therefore, our naive Hessian estimator is,

$\displaystyle\widetilde{H}_{n}=\frac{1}{n}\sum_{i=1}^{n}\frac{\widehat{G}_{i}+\widehat{G}_{i}^{\top}}{2}.$
(47)

where the $(\widehat{G}_{i}+\widehat{G}_{i}^{\top})/2\,$ term ensures the
symmetry of $\widetilde{H}_{n}$. The function query complexity is
$\mathcal{O}(m)$ per step for this Hessian estimation.

Now we restate our Lemma 4.1 for the both estimators (47) and (17).

######  Lemma B.1.

Under the assumptions in Theorem 3.3, we have the following result for the
Hessian estimator (47),

$\displaystyle\mathbb{E}\|\widetilde{H}_{n}-H\|^{2}\leq
C_{1}n^{-\alpha}+C_{2}\left(1+\frac{1}{m}\right)n^{-1}.$ (48)

The Hessian estimator (17) satisfies,

$\displaystyle\mathbb{E}\|\widetilde{H}_{n}-H\|^{2}\leq
C_{1}n^{-\alpha}+C_{2}p^{-1}n^{-1}.$ (49)

###### Proof.

In the case of naive Hessian estimator (47), we decompose
$\widetilde{H}_{n}-H$ as follows,

$\displaystyle\widetilde{H}_{n}-H$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\widehat{G}_{i}+\widehat{G}_{n}^{\top}}{2}-H$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{\widehat{G}_{i}+\widehat{G}_{i}^{\top}}{2}-\left(\frac{1}{m}\sum_{j=1}^{m}\boldsymbol{u}_{i}^{(j)}\boldsymbol{u}_{i}^{(j)\top}\right)\nabla^{2}f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})\left(\frac{1}{m}\sum_{j=1}^{m}\boldsymbol{v}_{i}^{(j)}\boldsymbol{v}_{i}^{(j)\top}\right)\right)$
$\displaystyle\quad+\frac{1}{n}\sum_{i=1}^{n}\left(\left(\frac{1}{m}\sum_{j=1}^{m}\boldsymbol{u}_{i}^{(j)}\boldsymbol{u}_{i}^{(j)\top}\right)\nabla^{2}f(\boldsymbol{\theta}_{n-1};\boldsymbol{\zeta}_{n})\left(\frac{1}{m}\sum_{j=1}^{m}\boldsymbol{v}_{i}^{(j)}\boldsymbol{v}_{i}^{(j)\top}\right)-\nabla^{2}f(\boldsymbol{\theta}_{i-1};\boldsymbol{\zeta}_{i})\right)$