An Asymptotic Analysis of Minibatch-Based Momentum Methods for Linear
Regression Models
Yuan Gaoa,b, Xuening Zhuc, Haobo Qi${}^{d},$ Guodong Lie, Riquan Zhanga,b,
Hansheng Wangd

a School of Statistics, East China Normal University, Shanghai, China
b Key Laboratory of Advanced Theory and Application in Statistics and Data
Science - MOE, East China Normal University, Shanghai, China
cSchool of Data Science, Fudan University, Shanghai, China
dGuanghua School of Management, Peking University, Beijing, China
eDepartment of Statistics and Actuarial Science, University of Hong Kong, Hong
Kong, China

###### Abstract

Momentum methods have been shown to accelerate the convergence of the standard
gradient descent algorithm in practice and theory. In particular, the
minibatch-based gradient descent methods with momentum (MGDM) are widely used
to solve large-scale optimization problems with massive datasets. Despite the
success of the MGDM methods in practice, their theoretical properties are
still underexplored. To this end, we investigate the theoretical properties of
MGDM methods based on the linear regression models. We first study the
numerical convergence properties of the MGDM algorithm and further provide the
theoretically optimal tuning parameters specification to achieve faster
convergence rate. In addition, we explore the relationship between the
statistical properties of the resulting MGDM estimator and the tuning
parameters. Based on these theoretical findings, we give the conditions for
the resulting estimator to achieve the optimal statistical efficiency.
Finally, extensive numerical experiments are conducted to verify our
theoretical results.

KEYWORDS: Gradient Descent, Momentum Method, Fixed Minibatch, Shuffled
Minibatch

1\. INTRODUCTION

Modern statistical analysis often involves massive datasets and sophisticated
models, which makes the corresponding statistical computation extremely
challenging. For example, the well-known ImageNet dataset, contains more than
one million photos and takes up approximately 150GB space on a hard drive
(Deng et al.,, 2009). Obviously, most computers available to researchers
cannot easily load such an enormous dataset as a whole into its memory. Then,
how to process these massive datasets efficiently with statistical guarantees
becomes a problem of great interest. To address this issue, it is in urgent
need to develop novel scalable algorithms for statistical analysis.

In recent yeas, a great number of researches have been devoted to dealing with
large-scale statistical optimization problems. For example, the distributed
computing methods, based on the divide-and-conquer strategy from computer
science literature, become very popular in the statistics community. Many
statistical learning problems have been investigated on a distributed system,
such as generalized linear models (Chen and Xie,, 2014; Battey et al.,, 2018;
Tang et al.,, 2020), quantile regression (Chen et al.,, 2019; Chen et al.,
2020b, ), kernel ridge regression (Chang et al.,, 2017; Lin et al.,, 2020),
support vector machine (Wang et al.,, 2019), principal component analysis (Fan
et al.,, 2019; Chen et al., 2021a, ), and references therein. Clearly, the
distributed computing frameworks reduce the computational burden of one single
machine by employing a computer cluster. However, it still replies on a
distributed architecture for implementation, which can be expensive for many
practical problems. With limited computational resources, it is desirable to
conduct computation with a single computer and feasible optimization
algorithms. In this regard, stochastic optimization methods, especially
stochastic gradient descent (SGD) and its derivatives, play an important role.

Before introducing SGD and its derivatives, we first discuss their classical
versions. As a first-order optimization method, gradient descent (GD) is
arguably the simplest one. Unlike second-order methods, such as Newton-Raphson
algorithm, a first-order method requires no computation of the Hessian
matrices. In fact, a standard GD algorithm only requires the gradients of the
objective function and a learning rate. Therefore, it is computationally more
feasible for sophisticated models with massive datasets. Furthermore, the cost
in computer memory for storing Hessian matrices can be significantly reduced.
Unfortunately, there is no free lunch. The tradeoff for GD methods’ excellent
practical feasibility is slow convergence. In other words, the convergence
rate of a GD algorithm is usually much slower than that of a Newton-Raphson
algorithm, if the Newton-Raphson algorithm can be implemented. This is
particularly true for ill-conditioned problems (Polyak,, 1987). In this case,
feature directions with tiny variability are likely to be ignored by a
standard GD algorithm. To partially fix this problem, Polyak, (1964) proposed
a heavy-ball method to accelerate the standard GD algorithm. Specifically, it
uses the information from not only the current step but also the previous step
(i.e., the momentum). As a consequence, the search direction can be slightly
modified. With this slight but novel modification, the convergence rate of
gradient descent with momentum (GDM) could be much improved compared to the
standard GD algorithm.

Despite its theoretical popularity, the practical implementation of GDM could
be nontrivial on a massive dataset. Under a traditional setup with small
sample size and simple model structure, a GDM algorithm can be easily
implemented by loading the entire dataset into the computer memory as a whole.
Subsequently, various gradients and momentums can be efficiently computed and
processed. Unfortunately, such a common practice could be seriously challenged
if the target model structure is sufficiently complex and the dataset is
massive in size. For example, the ResNet50 model (He et al.,, 2016) is a
classical deep learning model with more than 25 million parameters. In this
case, it would take a long time to compute the gradient and momentum needed
for performing one single step of iteration. Thus, implementing GDM algorithms
on such a complex model (i.e., the ResNet50 model) with such a massive dataset
(i.e., the ImageNet dataset) is very challenging.

Since a massive dataset cannot be easily loaded into computer memory as a
whole, it has to be processed in a batch-by-batch manner. Here, “batch” refers
to a subsample of the original whole dataset. Subsequently, the important
quantities (e.g., gradient and momentum) can be computed based on the batch.
This approach leads to various stochastic gradient descent (SGD) methods.
Here, “stochastic” refers to the batch being randomly subsampled from the
original whole dataset. Then, according to how this subsample (i.e., batch) is
generated and utilized, different SGD algorithms can be developed. The history
of SGD can date back to the pioneering work of Robbins and Monro, (1951). The
computational efficiency of SGD finds wide application in the large-scale
optimization problems, and attracts a large number of follow-up studies (see,
e.g., Wang et al.,, 2013, 2017; Yang et al.,, 2019; Chen et al., 2021b, , and
references therein). Although many existing works focus on the numerical
convergence rate of the stochastic optimization algorithms, there are also a
few of them that investigate the resulting estimator from a statistical
viewpoint. For example, Toulis et al., (2017) proposed an implicit SGD
procedure, and analyzed the asymptotic behavior of the estimators resulting
from the standard and implicit SGD procedures. Chen et al., 2020a  and Zhu et
al., (2021) investigated the statistical inference problem based on the
averaged SGD estimator, and provided several methods to estimate the
asymptotic covariance. Luo and Song, (2020) developed a renewable estimation
procedure that extends the standard SGD and averaged implicit SGD algorithm,
and further studied the asymptotic properties of the resulting estimator. In
general, the algorithms that they discussed are one-pass in nature and should
be very useful for streaming data analysis. However, for non-streaming
situations where the massive dataset is already on the hard drive, the one-
pass type of algorithm is less frequently used in practice. Instead, multiple
passes (or epochs) are more commonly used for better statistical efficiency.

In fact, the many popularly used SGD algorithms in practice are based on
shuffled minibatches. In other words, the entire SGD algorithm is divided into
different epochs. For each epoch, the original whole dataset is randomly
shuffled and then sequentially divided into different batches. Subsequently,
different batches are sequentially loaded into the computer memory and then
processed. Once all the batches have been successfully executed, the whole
dataset might need to be reshuffled. This process accomplishes one epoch of
the minibatch-based SGD algorithm. Accordingly, each sample in the whole
dataset is ensured to be used exactly once for every epoch iteration. Our
limited numerical experiments suggest that this approach often leads to
estimators with satisfactory statistical efficiency. As a consequence,
shuffled minibatch-based GD algorithms are indeed the most widely used SGD
algorithms in practice. In fact, it has been implemented by many important
software libraries, such as TensorFlow (Abadi et al.,, 2016) and PyTorch
(Paszke et al.,, 2017).

However, our theoretical understanding of this most popularly used SGD
algorithm (i.e., minibatch-based gradient descent with momentum, MGDM) is
extremely limited. A number of important theoretical questions remain unknown.
For example, under what conditions should an MGDM algorithm numerically
converge? Past literature suggests that for a standard GD algorithm, the
learning rate should not be too large. Otherwise, the algorithm might not
numerically converge. Consequently, it is natural to question whether this
remains true for MGDM algorithms. We are particularly interested in
understanding the important role played by the momentum parameter. Next,
assuming that the MGDM algorithm converges numerically and leads to the MGDM
estimator, what are the asymptotic properties of the resulting MGDM estimator?
This would facilitate the understanding of the roles played by both the
learning rate and momentum parameter.

To address these important questions, we start with a relatively simple model
(i.e., a linear regression model) and one particular type of MGDM algorithm
(i.e., the fixed minibatch-based GDM algorithm, FMGDM). By doing so, we are
able to analyze the theoretical properties of the MGDM method both numerically
and asymptotically. This leads to fruitful theoretical findings. In
particular, we identify the nearly sufficient and necessary conditions for the
tuning parameters (i.e., the learning rate and the momentum parameter) to
ensure numerical convergence of the algorithm. In addition, we further
investigate the numerical convergence rate and provide the optimal tuning
parameters specification in theory. Subsequently, we study the asymptotic
properties of the resulting MGDM estimator. We find that a diminishing
learning rate is needed for the resulting estimator to achieve global
statistical efficiency. Finally, we conduct extensive numerical experiments to
validate our theoretical findings.

In summary, we attempt to make the following important contributions to the
existing literature. First, we systematically study the MGDM algorithms for
both shuffled and fixed minibatch. A relatively complete theory is provided
for both numerical and asymptotic properties of the MGDM algorithm. Second,
our theory leads to nontrivial findings on the important relationship between
statistical efficiency of the resulting estimator and the tuning parameters.
In particular, we find that as the learning rate decreases, both MGDM
estimators (i.e., fixed and shuffled) can be statistically as efficient as the
global OLS estimator. Last, extensive numerical studies are presented to
verify our theoretical findings and demonstrate the outstanding performance of
our proposed FMGDM method.

The rest of the article is organized as follows. Section 2 first introduces
the MGDM method under the linear regression setup. Then, the numerical
convergence properties and the statistical efficiency of the resulting
estimator are investigated. Section 3 includes extensive numerical experiments
to corroborate our theoretical findings. We conclude the article in Section 4
and discuss a number of interesting topics for future study.

2\. MOMENTUM GRADIENT DESCENT WITH FIXED MINI-BATCH

2.1. Classical Gradient Descent with Momentum

We first introduce some frequently used notations for this paper. Let $A$ be a
matrix in $\mathbb{R}^{p_{1}\times p_{2}}$. We use $\lambda_{\max}(A)$ and
$\lambda_{\min}(A)$ to denote the largest and smallest eigenvalue ,
respectively, of $A$, if $A$ is symmetric. We define
$\|A\|=\sqrt{\lambda_{\max}(A^{\top}A)}$ and
$\|A\|_{F}=\sqrt{\mbox{tr}(A^{\top}A)}$ as the operator norm and the Frobenius
norm of $A$, respectively. In addition,
$\rho(A)=\max\\{|\lambda|:\lambda\text{ is the eigenvalue of }A\\}$, which
stands for the spectral radius of $A$. As usual, $\|v\|=\sqrt{v^{\top}v}$ is
the standard $l_{2}$-norm of an arbitrary vector $v\in\mathbb{R}^{p}$. We
start with a classical linear regression model. Let $(Y_{i},X_{i})$ be the
observation collected from $i$-th subject with $1\leq i\leq N$. Assume

$Y_{i}=X_{i}^{\top}\theta+\varepsilon_{i},\ \text{with}\ 1\leq i\leq N,$ (2.1)

where $X_{i}\in\mathbb{R}^{p}$ is the covariate, $Y_{i}\in\mathbb{R}$ is the
response, $\varepsilon_{i}$ is the random error with mean $0$ and variance
$\sigma^{2}$, and $\theta\in\mathbb{R}^{p}$ is the unknown parameter vector
with true value given by $\theta_{0}=(\theta_{01},\dots,\theta_{0p})^{\top}$.
Suppose the whole sample is indexed by $\mathbb{S}=\\{1,\dots,N\\}$. To
estimate the unknown parameter vector $\theta$, we usually minimize the
ordinary least squares (OLS) loss function as

$\mathcal{L}(\theta)=N^{-1}\sum_{i=1}^{N}\ell_{i}(\theta)=(2N)^{-1}\sum_{i=1}^{N}\Big{(}Y_{i}-X_{i}^{\top}\theta\Big{)}^{2},$

where $\ell_{i}(\theta)=2^{-1}(Y_{i}-X_{i}^{\top}\theta)^{2}$ is the loss
evaluated on the $i$-th observation. We know that the solution to this problem
is the classical ordinary least squares (OLS) estimator that has an explicit
form as
$\widehat{\theta}_{\text{ols}}=\arg\min_{\theta}\mathcal{L}(\theta)=\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xy}$,
where
$\widehat{\Sigma}_{xx}=N^{-1}\sum_{i=1}^{N}X_{i}X_{i}^{\top}\in\mathbb{R}^{p\times
p}$ and
$\widehat{\Sigma}_{xy}=N^{-1}\sum_{i=1}^{N}X_{i}Y_{i}\in\mathbb{R}^{p}$. Under
some mild conditions, one can show that
$\sqrt{N}(\widehat{\theta}_{\text{ols}}-\theta_{0})\to_{d}N(\boldsymbol{0},\sigma^{2}\Sigma_{xx}^{-1})$,
where $\Sigma_{xx}$ is the covariance matrix of $X_{i}$ (Rao et al.,, 1973).
To practically compute $\widehat{\theta}_{\text{ols}}$, the sample covariance
matrix $\widehat{\Sigma}_{xx}$ needs to be calculated and inverted. This leads
to a computational complexity of order $O(Np^{2}+p^{3})$ in general.

Another popular way to compute $\widehat{\theta}_{\text{ols}}$ is to use the
first-order optimization methods in an iterative way. Arguably, the most basic
first-order optimization method is the gradient descent (GD) algorithm. For
the OLS problem, a standard GD algorithm should update the estimates
iteratively as
$\widehat{\theta}^{(t)}=\widehat{\theta}^{(t-1)}-\alpha\dot{\mathcal{L}}(\widehat{\theta}^{(t-1)})$,
where $\dot{\mathcal{L}}(\theta)$ denotes the gradient of $\mathcal{L}$ at
$\theta$ and $\alpha>0$ is the learning rate controlling the step size. As
shown by Nesterov et al., (2018), with an appropriately selected $\alpha$, the
estimates $\\{\widehat{\theta}^{(t)}\\}$ generated by the GD algorithm should
satisfy
$\|\widehat{\theta}^{(t)}-\widehat{\theta}_{\text{ols}}\|\leq\rho^{t}\|\widehat{\theta}^{(0)}-\widehat{\theta}_{\text{ols}}\|$
for some convergence factor $\rho\in[0,1)$ and an arbitrary initial value
$\widehat{\theta}^{(0)}$. The optimal convergence factor can be achieved if we
choose
$\alpha=2/\\{\lambda_{\max}(\widehat{\Sigma}_{xx})+\lambda_{\min}(\widehat{\Sigma}_{xx})\\}$.
In this case, the optimal convergence factor should be
$\rho=(\kappa-1)/(\kappa+1)$ , where
$\kappa=\lambda_{\max}(\widehat{\Sigma}_{xx})/\lambda_{\min}(\widehat{\Sigma}_{xx})\geq
1$ is the condition number of $\widehat{\Sigma}_{xx}$.

However, when the condition number $\kappa$ is large, the convergence factor
$\rho$ could be very close to $1$. Consequently, the convergence rate could be
painfully slow. To accelerate it, Polyak, (1964) proposed the heavy-ball
method, which is also called the momentum method (Sutskever et al.,, 2013).
Specifically, in this work, we investigate the gradient descent with the
classical momentum, which updates the estimates as

$\displaystyle\widehat{\theta}^{(t)}=\widehat{\theta}^{(t-1)}-v^{(t)},$ (2.2)
$\displaystyle v^{(t)}=\gamma
v^{(t-1)}+\alpha\dot{\mathcal{L}}(\widehat{\theta}^{(t-1)}),$ (2.3)

where $\alpha>0$ is the learning rate and $\gamma>0$ is the momentum
parameter. Polyak, (1964) showed that, by choosing appropriate tuning
parameters $\alpha$ and $\gamma$, the optimal convergence factor of the
momentum method could be achieved as
$\rho=(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$. In contrast, that of a standard GD
algorithm is $\rho=(\kappa-1)/(\kappa+1)$, which is strictly larger than that
of the GD with momentum (GDM) method as long as $\kappa>1$. This suggests that
the optimal convergence rate of a GDM algorithm should be faster than that of
a standard GD algorithm.

2.2. A Fixed Mini-batch Based Gradient Descent with Momentum

Despite its theoretical attractiveness, a standard GDM method cannot be
immediately used to handle massive datasets. In many cases, the dataset could
be too large to be comfortably processed as a whole by computer memory. As a
consequence, it has to be placed on the hard drive and then processed in a
batch-by-batch manner. In practice, this approach leads to the most popularly
used minibatch-based gradient descent methods with momentum (MGDM).
Specifically, to practically implement an MGDM algorithm, we first need to
randomly partition the whole sample $\mathbb{S}$ into $M$ disjoint minibatches
as $\mathbb{S}=\cup_{m=1}^{M}\mathcal{S}_{(m)}$ with
$\mathcal{S}_{(m_{1})}\cap\mathcal{S}_{(m_{2})}=\emptyset$ for any $m_{1}\neq
m_{2}$. Throughout the article, we assume that $M$ is fixed. Without loss of
generality, we assume that $|\mathcal{S}_{(m)}|=n$ for each $1\leq m\leq M$.
Obviously, we should have $N=Mn$. As mentioned before, we start with the fixed
minibatch-based GDM (FMGDM) method due to its analytical simplicity. In other
words, we assume that the random partition
$\mathbb{S}=\cup_{m=1}^{M}\mathcal{S}_{(m)}$, once given, should be fixed
throughout the rest of the algorithm.

In the fixed minibatch setting, the OLS loss function based on
$\mathcal{S}_{(m)}$ can be defined as
$\mathcal{L}_{(m)}(\theta)=(2n)^{-1}\sum_{i\in\mathcal{S}_{(m)}}(Y_{i}-X_{i}^{\top}\theta)^{2}$.
The corresponding gradient is
$\dot{\mathcal{L}}_{(m)}(\theta)=\widehat{\Sigma}_{xx}^{(m)}\theta-\widehat{\Sigma}_{xy}^{(m)}$,
where
$\widehat{\Sigma}_{xx}^{(m)}=n^{-1}\sum_{i\in\mathcal{S}_{(m)}}X_{i}X_{i}^{\top}\in\mathbb{R}^{p\times
p}$ and
$\widehat{\Sigma}_{xy}^{(m)}=n^{-1}\sum_{i\in\mathcal{S}_{(m)}}X_{i}Y_{i}\in\mathbb{R}^{p}$.
Subsequently, a standard GDM algorithm should be executed in a batch-by-batch
manner. Specifically, let $\widehat{\theta}^{(t,m)}$ be the estimate obtained
from the $m$-th minibatch in the $t$-th epoch. We then should have

$\displaystyle\widehat{\theta}^{(t,m)}=\widehat{\theta}^{(t,m-1)}-v^{(t,m)},$
(2.4) $\displaystyle v^{(t,m)}=\gamma
v^{(t,m-1)}+\alpha\dot{\mathcal{L}}_{(m)}(\widehat{\theta}^{(t,m-1)}),$ (2.5)

where $\widehat{\theta}^{(t,0)}=\widehat{\theta}^{(t-1,M)}$,
$v^{(t,0)}=v^{(t-1,M)}$, and $\widehat{\theta}^{(0,m)}$ ($1\leq m\leq M$) is
the initial estimate specified for the $m$-th minibatch. This immediately
leads to three important questions. First, theoretically, if we can treat
(2.4) and (2.5) as a linear dynamical system, does this dynamical system have
a stable solution? Second, does the FMGDM algorithm converge to this stable
solution? Last, if the FMGDM algorithm does converge to the stable solution
under appropriate conditions, what is the statistical properties of the
resulting estimator?

2.3. The Linear Dynamical System and Stable Solution

We investigate the first question in this subsection: does there exist a
stable solution to the FMGDM algorithm? To address this important problem, we
temporarily assume that there does exist a numerical limit
$\widehat{\theta}^{(m)}$ such that
$\widehat{\theta}^{(t,m)}\to\widehat{\theta}^{(m)}$ as $t\to\infty$ for each
$1\leq m\leq M$. This leads to an analytical solution as
$\widehat{\theta}^{(m)}$ for each $1\leq m\leq M$. This analytical solution
further helps us to understand the conditions needed for its unique existence.
Specifically, we write the stable solution vector as
$\widehat{\theta}^{*}=(\widehat{\theta}^{(1)\top},\dots,\widehat{\theta}^{(M)\top})^{\top}\in\mathbb{R}^{q}$
with $q=Mp$. Then, by (2.4) and (2.5), we know that $\widehat{\theta}^{*}$
should be the solution to the following linear dynamical system,

$\displaystyle\begin{cases}\widehat{\theta}^{(t,1)}=\widehat{\Delta}^{(1)}\widehat{\theta}^{(t,M)}-\gamma\widehat{\theta}^{(t,M-1)}+\alpha\widehat{\Sigma}_{xy}^{(1)}\\\
\widehat{\theta}^{(t,2)}=\widehat{\Delta}^{(2)}\widehat{\theta}^{(t,1)}-\gamma\widehat{\theta}^{(t,M)}+\alpha\widehat{\Sigma}_{xy}^{(2)}\\\
\cdots\\\
\widehat{\theta}^{(t,M)}=\widehat{\Delta}^{(M)}\widehat{\theta}^{(t,M-1)}-\gamma\widehat{\theta}^{(t,M-2)}+\alpha\widehat{\Sigma}_{xy}^{(t,M)},\\\
\end{cases}$ (2.6)

where
$\widehat{\Delta}^{(m)}=(1+\gamma)I_{p}-\alpha\widehat{\Sigma}_{xx}^{(m)}\in\mathbb{R}^{p\times
p}$ for $1\leq m\leq M$. We further write
$\widehat{\Sigma}_{xy}^{*}=\big{(}\widehat{\Sigma}_{xy}^{(1)\top},\dots,\widehat{\Sigma}_{xy}^{(M)\top}\big{)}^{\top}\in\mathbb{R}^{q}$,
and

$\displaystyle\widehat{\Omega}=\begin{bmatrix}I_{p}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\gamma
I_{p}&-\widehat{\Delta}^{(1)}\\\
-\widehat{\Delta}^{(2)}&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\gamma
I_{p}\\\ \gamma
I_{p}&-\widehat{\Delta}^{(3)}&I_{p}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\gamma
I_{p}&-\widehat{\Delta}^{(M)}&I_{p}\end{bmatrix}\in\mathbb{R}^{q\times q}.$
(2.7)

Then, we should have
$\widehat{\Omega}\widehat{\theta}^{*}=\alpha\widehat{\Sigma}_{xy}^{*}$ or
equivalently
$\widehat{\theta}^{*}=\alpha\widehat{\Omega}^{-1}\widehat{\Sigma}_{xy}^{*}$,
provided $\widehat{\Omega}$ is invertible. Consequently, whether
$\widehat{\Omega}$ is invertible or not determines the existence of the stable
solution. Then, in the following theorem, we elaborate the conditions to
ensure its invertibility in probability.

######  Theorem 1.

(Stable Solution) Assume $\Sigma_{xx}\in\mathbb{R}^{p\times p}$ is a positive
definite matrix with eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{p}>0$. Further assume that
$\gamma\neq 1$, $\alpha\neq 0$, and $\alpha\neq 2(1+\gamma)/\lambda_{j}$ for
each $1\leq j\leq p$. Then, $\widehat{\Omega}$ is invertible with probability
tending to one as $n\to\infty$.

The proof of Theorem 1 can be found in Appendix A.1. By Theorem 1, we know
that the stable solution does exist for a wide range of different $\alpha$ and
$\gamma$ values. In fact, as long as $\gamma\neq 1$, $\alpha\neq 0$, and
$\gamma\neq 2(1+\gamma)/\lambda_{j}$ for every $1\leq j\leq p$, we should have
$\widehat{\Omega}$ invertible with high probability, and thus, the stable
solution should exist. In particular, we find that the stable solution exists
as $\gamma\to\infty$. Nevertheless, we should remark that the existence of the
stable solution and the FMGDM estimator might be totally different. By the
FMGDM estimator, we mean that the FMGDM algorithm should converge to a
numerical limit $\widehat{\theta}^{(m)}$, and this numerical limit
$\widehat{\theta}^{(m)}$ is referred to as the FMGDM estimator. Obviously, as
long as such a numerical limit does exist, it must be a stable solution to
(2.6). Nevertheless, this does not imply that every possible stable solution
can be numerically approached by the FMGDM algorithm. In fact, there do exist
stable solutions (e.g., those stable solutions associated with very large
$\gamma$ values) that cannot be achieved using the FMGDM algorithm.
Consequently, it is not a FMGDM estimator. Thus, we know that the existence of
the stable solution should be a necessary condition for the existence of the
FMGDM estimator but might not be sufficient. Accordingly, we investigate the
conditions to assure the numerical convergence of the FMGDM algorithm in the
next subsection.

2.4. Numerical Convergence Properties

We devote this subsection to studying the second important problem, that is,
the numerical convergence properties of the FMGDM algorithm. Specifically, we
are interested in studying how the $\alpha$ and $\gamma$ would affect the
numerical convergence of the FMGDM algorithm. Recall that $\lambda_{j}\ (1\leq
j\leq p)$ is the $j$-th largest eigenvalue of $\Sigma_{xx}$. The main results
are summarized in the following theorem.

######  Theorem 2.

(Numerical Convergence) (a) If $0\leq\gamma<1$ and
$0<\alpha<2(1+\gamma)/\lambda_{1}$ and assuming $n\to\infty$, then, we have
$\widehat{\theta}^{(t,m)}\to\widehat{\theta}^{(m)}$ as $t\to\infty$ holds with
probability tending to one. (b) Assume
$\widehat{\theta}^{(0,m)}\neq\widehat{\theta}^{(m)}$ and $n\to\infty$. For
$\gamma>1$ or $\alpha>2(1+\gamma)/\lambda_{1}$, we have
$\widehat{\theta}^{(t,m)}\not\to\widehat{\theta}^{(m)}$ as $t\to\infty$ holds
with probability tending to one.

The proof of Theorem 2 can be found in Appendix A.2. Theorem 2 provides nearly
sufficient and necessary conditions for numerical convergence. Specifically,
from Theorem 2(a), we know that to ensure the numerical convergence, we must
have $0\leq\gamma<1$. This corroborates the common practice in the existing
literature very well (Sutskever et al.,, 2013). Otherwise, by Theorem 2(b), we
know that the numerical convergence can never be achieved if $\gamma>1$.
Consequently, the FMGDM estimator does not exist in this case. Furthermore,
with a given momentum parameter $0\leq\gamma<1$, the learning rate $\alpha$
should be sufficiently small. Otherwise, the algorithm should not converge.
Remarkably, the upper bound of the learning rate $\alpha$ is related to the
largest eigenvalue of $\Sigma_{xx}$ (i.e., $\lambda_{1}$) and the momentum
parameter $\gamma$. This upper bound should be larger if the momentum
parameter $\gamma$ is larger. As a result, this enables a larger learning rate
(thus faster convergence) to be used for a FMGDM algorithm compared to a
standard GD algorithm ($\gamma=0$). Theorem 2 examines the conditions for the
numerical convergence for the FMGDM algorithm. However, it does not provide
the exact convergence rate. In particular, we want to understand how the
$\alpha$ and $\gamma$ specification would affect the numerical convergence
rate. We next develop the following theorem.

######  Theorem 3.

(Numerical Convergence Speed) If $0\leq\gamma<1$ and
$0<\alpha<2(1+\gamma)/\lambda_{1}$, then there exists a constant
$\rho\in[\gamma,1)$, such that
$\|\widehat{\theta}^{(t,m)}-\widehat{\theta}^{(m)}\|\leq(\rho^{M}+\zeta_{n}+\eta_{t})^{t}\big{(}\|\widehat{\theta}^{(0,m)}-\widehat{\theta}^{(m)}\|+\|\widehat{\theta}^{(0,m-1)}-\widehat{\theta}^{(m-1)}\|\big{)}$
for each $1\leq m\leq M$, where $\zeta_{n}\to_{p}0$ as $n\to\infty$ and
$\eta_{t}\to 0$ as $t\to\infty$. Furthermore, we have

  1. (a)

The minimal (i.e., optimal) $\rho$ is given by
$\rho=(\sqrt{\lambda_{1}/\lambda_{p}}-1)/(\sqrt{\lambda_{1}/\lambda_{p}}+1)$,
when $\alpha=4/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$ and
$\gamma=(\sqrt{\lambda_{1}}-\sqrt{\lambda_{p}})^{2}/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$.

  2. (b)

If $\gamma=0$, we should have
$\rho=\max\\{|1-\alpha\lambda_{1}|,|1-\alpha\lambda_{p}|\\}$. In particular,
the minimal (i.e., optimal) $\rho$ is given by
$\rho=(\lambda_{1}/\lambda_{p}-1)/(\lambda_{1}/\lambda_{p}+1)$, when
$\alpha=2/(\lambda_{1}+\lambda_{p})$.

  3. (c)

Assume $0<\alpha<1/\lambda_{1}$. (c.1) If $\gamma=0$, then the $\rho$ is given
by $\rho=1-\alpha\lambda_{p}$. (c.2) If
$0<\gamma<(1-{\alpha\lambda_{p}})^{2}$, then $\rho<1-\alpha\lambda_{p}$. In
particular, the minimal (i.e., optimal) $\rho$ is given by
$\rho=1-\sqrt{\alpha\lambda_{p}}$, when
$\gamma=(1-\sqrt{\alpha\lambda_{p}})^{2}$.

The proof of Theorem 3 is shown in Appendix A.3. From Theorem 3, we know that
the FMGDM algorithm should enjoy a linear convergence rate by choosing
appropriate $\alpha$ and $\gamma$, as long as both $n$ and $t$ are large
enough. By Theorem 3 (a), we know that the optimal convergence factor
$\rho=(\sqrt{\lambda_{1}/\lambda_{p}}-1)/(\sqrt{\lambda_{1}/\lambda_{p}}+1)$
can be achieved by carefully specified $\alpha$ and $\gamma$. This
specification is closely related to the eigenvalue structure of the covariance
matrix $\Sigma_{xx}$. Our findings are in line with theoretical findings in
the classical GDM method without minibatch (Polyak,, 1964). It is known that
this convergence factor also attains the achievable lower bound of various
first-order methods (Nesterov et al.,, 2018). Moreover, Theorem 3 (b) claims
that, for a standard GD algorithm with $\gamma=0$, the best convergence factor
should be $\rho=(\lambda_{1}/\lambda_{p}-1)/(\lambda_{1}/\lambda_{p}+1)$. This
represents a convergence rate that is strictly slower than that of the FMGDM
algorithm as long as $\lambda_{1}>\lambda_{p}$.

Although the optimal numerical convergence can be guaranteed by (a) and (b)
under different circumstances, the stable solution is not necessarily optimal
in statistical efficiency. For most SGD-related methods, the decaying learning
rate $\alpha$ is typically needed to improve statistical efficiency (Robbins
and Monro,, 1951; Polyak and Juditsky,, 1992). In this case, from Theorem 3
(c.1) we know that the convergence factor in this case tends to $1$ as
$\alpha\to 0$, which represents a very slow convergence rate. Moreover, from
Theorem 3 (c.2) we find that by choosing an appropriate $\gamma$, the
convergence factor of the FMGDM could be strictly smaller than that of a
standard GD algorithm ($\gamma=0$). This finding suggests that the momentum
term can indeed accelerate the numerical convergence. In addition, Theorem 3
(c.2) discovers a very interesting and novel interactive relationship between
the learning rate $\alpha$ and the momentum parameter $\gamma$. Specifically,
it suggests that $\gamma$ should approach $1$ as $\alpha\to 0$. This result
theoretically explains why practitioners often like to specify $\gamma$ to be
a value very close to $1$ to strike a balance between both numerical and
statistical convergence (Sutskever et al.,, 2013).

2.5. The Statistical Properties of the Stable Solution

In this subsection, we investigate how the tuning parameters $\alpha$ and
$\gamma$ would affect the statistical properties of the stable solution. Note
that the OLS estimator is the global optimal solution to our problem. Ideally,
we should have the resulting FMGDM estimator stay with the OLS estimator as
closely as possible. Consequently, we should examine the relationship between
the stable solution and the OLS estimator. To this end, we need to analyze
$\widehat{\Omega}$ and its inverse (assume it exists) carefully. Accordingly,
we decompose it into $\widehat{\Omega}=A+\alpha\widehat{B}$, where

$\displaystyle
A=\begin{bmatrix}I_{p}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\gamma
I_{p}&-(1+\gamma)I_{p}\\\
-(1+\gamma)I_{p}&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\gamma
I_{p}\\\ \gamma
I_{p}&-(1+\gamma)I_{p}&I_{p}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\gamma
I_{p}&-(1+\gamma)I_{p}&I_{p}\end{bmatrix}\in\mathbb{R}^{q\times q},$ (2.8)

and

$\displaystyle\widehat{B}=\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(1)}\\\
\widehat{\Sigma}_{xx}^{(2)}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(M)}&\boldsymbol{0}\end{bmatrix}\in\mathbb{R}^{q\times
q}.$ (2.9)

Note that $A$ and $\widehat{B}$ do not involve the learning rate $\alpha$.
Furthermore, it can be easily verified that
$A^{\top}I^{*}=AI^{*}=\boldsymbol{0}$, where $I^{*}=\boldsymbol{1}_{M}\otimes
I_{p}$. Thus, the columns of $I^{*}$ are the eigenvectors of $A$ and
$A^{\top}$ corresponding to the eigenvalue $0$. Next, let
$P_{1}=I^{*}/\sqrt{M}\in\mathbb{R}^{q\times p}$. Then, we can construct
$P_{2}\in\mathbb{R}^{q\times(q-p)}$ that has orthonormal columns and is
orthogonal to $P_{1}$ (i.e., $P_{2}^{\top}P_{2}=I_{p},\
P_{2}^{\top}P_{1}=\boldsymbol{0}$). We then have
$P=[P_{1},P_{2}]\in\mathbb{R}^{q\times q}$ as an orthogonal matrix. Simple
algebra gives
$P^{\top}AP=[\boldsymbol{0},\boldsymbol{0};\boldsymbol{0},P_{2}^{\top}AP_{2}]$.
Consequently, we obtain that

$\displaystyle P^{\top}\widehat{\Omega}P=P^{\top}AP+\alpha
P^{\top}\widehat{B}P=\begin{bmatrix}\alpha P_{1}^{\top}\widehat{B}P_{1}&\alpha
P_{1}^{\top}\widehat{B}P_{2}\\\ \alpha
P_{2}^{\top}\widehat{B}P_{1}&P_{2}^{\top}AP_{2}+\alpha
P_{2}^{\top}\widehat{B}P_{2}\end{bmatrix}.$ (2.10)

Denote $\widehat{B}_{11}^{P}=P_{1}^{\top}\widehat{B}P_{1}$,
$\widehat{B}_{12}^{P}=P_{1}^{\top}\widehat{B}P_{2}$,
$\widehat{B}_{21}^{P}=P_{2}^{\top}\widehat{B}P_{1}$,
$\widehat{B}_{22}^{P}=P_{2}^{\top}\widehat{B}P_{2}$, and
$A_{22}^{P}=P_{2}^{\top}AP_{2}$. Then, (2.10) is simplified to
$P^{\top}\widehat{\Omega}P=[\alpha\widehat{B}_{11}^{P},\alpha\widehat{B}_{12}^{P};\,\alpha\widehat{B}_{21}^{P},A_{22}^{P}+\alpha\widehat{B}_{22}^{P}]$.
With the help of this transformation, we can obtain the following theorem,
which elaborates the difference between the stable solution and the OLS
estimator.

######  Theorem 4.

(Relationship with the OLS Estimator) Suppose the conditions in Theorem 1
hold. Further denote
$\widehat{\theta}_{\text{ols}}^{*}=\mbox{\boldmath$1$}_{M}\otimes\widehat{\theta}_{\text{ols}}$.
Then we have

  1. (a)

If $\gamma\geq 0$ is fixed and $\alpha\to 0$, then
$\widehat{\theta}^{*}=\widehat{\theta}_{\text{ols}}^{*}+\alpha\widehat{E}\\{1+O_{p}(\alpha)\\}$,
where
$\widehat{E}=\big{(}P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}-I_{q}\big{)}\big{\\{}P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}\big{\\}}\big{(}\widehat{B}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{*}\big{)}$.

  2. (b)

If $\alpha>0$ is fixed and $\gamma\to\infty$ , then
$\widehat{\theta}^{*}=\widehat{\theta}_{\text{ols}}^{*}+\alpha\widehat{E}\\{1+O_{p}(\|(A_{22}^{P})^{-1}\|)\\}$,
where $\widehat{E}$ is the same as that in (a) and $\|(A_{22}^{P})^{-1}\|\to
0$ as $\gamma\to\infty$.

The proof of Theorem 4 can be found in Appendix A.4. With a fixed $\gamma$, it
can be verified that $\sqrt{n}\widehat{E}=O_{p}(1)$ as $n\to\infty$. Then by
Theorem 4, we should have $\alpha\sqrt{n}\widehat{E}\to_{p}0$ if $\alpha\to
0$. Hence it can be ignored since we have
$\sqrt{n}(\widehat{\theta}_{\text{ols}}-\theta_{0})=O_{p}(1)$. With a fixed
$\alpha$, it can be verified that $(A_{22}^{P})^{-1}\to 0$ as
$\gamma\to\infty$. This further implies that $\sqrt{n}\widehat{E}\to_{p}0$ as
$\gamma\to\infty$. Consequently, we should have
$\widehat{\theta}^{(m)}\to\widehat{\theta}_{\text{ols}}$ for $1\leq m\leq M$
as long as $\alpha\to 0$ or $\gamma\to\infty$. This immediately suggests two
different ways to force the stable solution to stay close to the global OLS
estimator. They are, respectively, $\alpha\to 0$ and $\gamma\to\infty$.
However, from Theorem 2 (b), we know that the FMGDM algorithm does not
converge if $\gamma>1$. Consequently, this leaves us only one way to obtain
the best statistical efficiency, which is forcing $\alpha\to 0$. On the other
hand, if $\alpha$ is set to be too small, then the FMGDM algorithm would
converge painfully slowly. Thus, practitioners have to consider relatively
large $\alpha$ values at the early stage of the algorithm for faster numerical
convergence. Then we turn to smaller values for improved statistical
efficiency. Moreover, as suggested by Theorem 3 (c.2), we can choose a
momentum parameter close to $1$ to accelerate the numerical convergence.

To further investigate the statistical property of the stable solution, we
establish its asymptotic normality for fixed $\gamma$ in the following
theorem.

######  Theorem 5.

(Asymptotic Normality) Suppose the conditions in Theorem 1 hold. Further
assume that $\gamma\geq 0$ is fixed and $\alpha\to 0$. This yields

$\sqrt{N}\big{(}\widehat{\theta}^{(m)}-\theta_{0}+\widehat{b}_{m}(\alpha)\big{)}\to_{d}N_{p}\Big{(}\boldsymbol{0},\sigma^{2}V_{m}(\alpha)\Big{)}$

as $N\to\infty$, where $\widehat{b}_{m}(\alpha)=O_{p}(\alpha^{2}n^{-3/2})$,
$V_{m}(\alpha)=\Sigma_{xx}^{-1}+\alpha^{2}Md_{\gamma}^{2}\Sigma_{xx}+O(\alpha^{4})$,
and $d_{\gamma}$ is the $l_{2}$-norm of the first row of the matrix
$P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}$.

The proof of Theorem 5 is shown in Appendix A.5 and the explicit formulas for
calculating $d_{\gamma}$ can be found in Appendix A.6. Theorem 5 gives the
asymptotic behavior of the stable solution. Specifically, for fixed $\gamma$
and $M$, the stable solution is still $\sqrt{N}$-consistent as long as
$\sqrt{N}\widehat{b}_{m}(\alpha)=o_{p}(1)$. In addition, from the bias and
variance term, we can see that as $\alpha\to 0$, the stable solution should be
asymptotically equivalent to the OLS estimator.

2.6. Shuffled Minibatch-Based Gradient Descent with Momentum

As mentioned before, the shuffled minibatch-based GDM (SMDGM) algorithm is in
fact the most popularly used MGDM algorithm in practice, but little is
understood regarding its statistical properties. However, Theorem 4 developed
in the previous subsection provides us an opportunity to obtain some
theoretical insights. Specifically, let $\mathcal{S}^{(k,m)}$ be the $m$-th
minibatch used in the $k$-th shuffled partition ($1\leq k\leq K$). We should
have $\mathbb{S}=\cup_{m=1}^{M}\mathcal{S}^{(k,m)}$ for every $k$ and
$\mathcal{S}^{(k,m_{1})}\cap\mathcal{S}^{(k,m_{2})}=\emptyset$ for any
$m_{1}\neq m_{2}$ and every $k$. Different shuffled partitions subsequently
determine different dynamical systems and thus different stable solutions.
More specifically, consider one particular shuffled partition
$\\{\mathcal{S}^{(k,m)}:1\leq m\leq M\\}$ for one particular $1\leq k\leq K$.
Then, the corresponding dynamical system becomes

$\displaystyle\begin{cases}\widehat{\theta}^{(t,1)}=\widehat{\Delta}^{(k,1)}\widehat{\theta}^{(t,M)}-\gamma\widehat{\theta}^{(t,M-1)}+\alpha\widehat{\Sigma}_{xy}^{(k,1)}\\\
\widehat{\theta}^{(t,2)}=\widehat{\Delta}^{(k,2)}\widehat{\theta}^{(t,1)}-\gamma\widehat{\theta}^{(t,M)}+\alpha\widehat{\Sigma}_{xy}^{(k,2)}\\\
\cdots\\\
\widehat{\theta}^{(t,M)}=\widehat{\Delta}^{(k,M)}\widehat{\theta}^{(t,M-1)}-\gamma\widehat{\theta}^{(t,M-2)}+\alpha\widehat{\Sigma}_{xy}^{(k,M)},\\\
\end{cases}$ (2.11)

where
$\widehat{\Delta}^{(k,m)}=(1+\gamma)I_{p}-\alpha\widehat{\Sigma}_{xx}^{(k,m)}\in\mathbb{R}^{p\times
p}$,
$\widehat{\Sigma}_{xx}^{(k,m)}=n^{-1}\sum_{i\in\mathcal{S}^{(k,m)}}X_{i}X_{i}^{\top}$,
and
$\widehat{\Sigma}_{xy}^{(k,m)}=n^{-1}\sum_{i\in\mathcal{S}^{(k,m)}}X_{i}Y_{i}$
for each $1\leq m\leq M$. We denote the corresponding stable solution by
$\widehat{\theta}^{(k)*}$. Then, from Theorem 4, we know that
$\widehat{\theta}^{(k)*}-\widehat{\theta}_{\text{ols}}^{*}=\alpha\widehat{E}^{(k)}\\{1+O_{p}(\alpha)\\}$,
where
$\widehat{E}^{(k)}=\big{(}P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}-I_{q}\big{)}\big{\\{}P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}\big{\\}}\big{(}\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\big{)}$,
$\widehat{B}^{(k)}$ has the same form as (2.9) but with
$\widehat{\Sigma}_{xx}^{(m)}$ replaced by $\widehat{\Sigma}_{xx}^{(k,m)}$ for
each $1\leq m\leq M$, and
$\widehat{\Sigma}_{xy}^{(k)*}=\big{(}\widehat{\Sigma}_{xy}^{(k,1)\top},\dots,\widehat{\Sigma}_{xy}^{(k,M)\top}\big{)}^{\top}$.
Consequently, as long as we can show that $\max_{1\leq k\leq
K}\|\widehat{E}^{(k)}\|=o_{p}(1/\sqrt{N})$, then the difference between the
shuffled estimator $\widehat{\theta}^{(k)*}$ and the global OLS estimator is
ignorable uniformly over all shuffled partitions.

Specifically, we first note that $\max_{1\leq k\leq
K}\|\widehat{E}^{(k)}\|\leq\big{(}\max_{1\leq k\leq
K}\big{\|}P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}\big{\|}+1\big{)}\
\big{\|}P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}\big{\|}\max_{1\leq k\leq
K}\big{\|}\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\big{\|}$,
where

$\displaystyle
P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}=M^{-1}\begin{bmatrix}\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,2)}&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,3)}&\cdots&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,1)}\\\
\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,2)}&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,3)}&\cdots&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,1)}\\\
\vdots&\vdots&\vdots&\vdots\\\
\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,2)}&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,3)}&\cdots&\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,1)}\\\
\end{bmatrix},$

and
$\|\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\|=\sum_{m=1}^{M}\|\widehat{\Sigma}_{xx}^{(k,m)}\widehat{\theta}_{\text{ols}}-\widehat{\Sigma}_{xy}^{(k,m)}\|$.
For simplicity, we next assume that both the covariate $X_{i}$ and the random
error $\varepsilon_{i}$ are normally distributed. We are then able to study
the asymptotic behavior of $\max_{1\leq k\leq K}\|\widehat{E}^{(k)}\|$. To
this end, we need to study $\max_{1\leq k\leq
K}\|P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}\|$ and
$\max_{1\leq k\leq
K}\big{\|}\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\|$
separately. We first consider $\max_{1\leq k\leq
K}\|P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}\|$. By Lemma
A.3 in Bickel et al., (2008), there should exist constants $C_{1}>0$ and
$C_{2}>0$ such that $P(\max_{1\leq m\leq
M}\|\widehat{\Sigma}_{xx}^{(k,m)}-\Sigma_{xx}\|_{F}\geq\nu)\leq
C_{1}\exp(-C_{2}N\nu^{2})$ for any sufficiently small $\nu>0$, as long as $M$
and $p$ are fixed. Consequently, we should have
$P(\max_{m,k}\|\widehat{\Sigma}_{xx}^{(k,m)}-\Sigma_{xx}\|_{F}\geq\nu)\leq
C_{1}\exp(\log K-C_{2}N\nu^{2})$, which tends to $0$ as long as $\log K=o(N)$.
Therefore, we can derive that
$\max_{m,k}\|\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xx}^{(k,m)}\|\leq\|\widehat{\Sigma}_{xx}^{-1}\|\max_{m,k}\|\widehat{\Sigma}_{xx}^{(k,m)}\|_{F}=O_{p}(1)$.
As a consequence, we should have $\max_{1\leq k\leq
K}\|P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}^{(k)}\|=O_{p}(1)$
as long as $\log K=o(N)$.

Next, we study $\max_{1\leq k\leq
K}\big{\|}\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\|$.
Let
$\widehat{\Sigma}_{x\varepsilon}=N^{-1}\sum_{i=1}^{N}X_{i}\varepsilon_{i}$,
and
$\widehat{\Sigma}_{x\varepsilon}^{(k,m)}=n^{-1}\sum_{i\in\mathcal{S}^{(k,m)}}X_{i}\varepsilon_{i}$.
Then, we should have
$\widehat{\Sigma}_{xy}^{(k,m)}=\widehat{\Sigma}_{xx}^{(k,m)}\theta_{0}+\Sigma_{x\varepsilon}^{(k,m)}$,
and
$\widehat{\theta}_{\text{ols}}=\theta_{0}+\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{x\varepsilon}$.
Note that
$\|\widehat{\Sigma}_{xx}^{(k,m)}\widehat{\theta}_{\text{ols}}-\widehat{\Sigma}_{xy}^{(k,m)}\|=\|\widehat{\Sigma}_{xx}^{(k,m)}\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{x\varepsilon}-\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|\leq\|\widehat{\Sigma}_{xx}^{(k,m)}\|\|\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{x\varepsilon}\|+\|\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|$.
Then by Bernstein’s inequality (Vershynin,, 2018, Theorem 2.8.2), we have
$P(\sqrt{N/\log K}\max_{1\leq m\leq
M}\|\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|\geq\nu)\leq
C_{1}\exp\big{[}-C_{2}\min\\{(\log K)\nu^{2},\sqrt{N\log K}\nu\\}\big{]}$ for
any $\nu>0$ as long as $p$ and $M$ are fixed, where $C_{1},C_{2}$ are two
positive constants. Hence, it follows that $P(\sqrt{N/\log
K}\max_{m,k}\|\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|\geq\nu)\leq
C_{1}\exp\big{[}\log K-C_{2}\min\\{(\log K)\nu^{2},\sqrt{N\log
K}\nu\\}\big{]}$ for any $\nu>0$. As a consequence, for any $\delta>0$, there
should be a sufficiently large constant $\nu_{\delta}>0$ such that
$P(\sqrt{N/\log
K}\max_{m,k}\|\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|\geq\nu_{\delta})\leq\delta$
for each $N$, as long as $\log K=o(N)$. Thus, we should have
$\max_{m,k}\|\widehat{\Sigma}_{x\varepsilon}^{(k,m)}\|=O_{p}(\sqrt{(\log
K)/N})$. Since $\max_{k,m}\|\widehat{\Sigma}_{xx}^{(k,m)}\|=O_{p}(1)$ and
$\|\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{x\varepsilon}\|=O_{p}(1/\sqrt{N})$,
we can conclude that $\max_{1\leq k\leq
K}\big{\|}\widehat{B}^{(k)}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{(k)*}\|=O_{p}(\sqrt{(\log
K)/N})$ as long as $\log K=o(N)$.

The results above demonstrate that $\max_{1\leq k\leq
K}\|\widehat{E}^{(k)}\|=O_{p}(\sqrt{(\log K)/N})$ provided $\log K=o(N)$.
These results further suggest that differences between the SMGDM estimators
and the global OLS estimator are uniformly bounded by $o_{p}(1/\sqrt{N})$, as
long as $\alpha\sqrt{\log K}\to 0$ and $\log K=o(N)$. Consequently, they
should share the same statistical efficiency as the global OLS estimator.

3\. NUMERICAL STUDIES

In this section, we conduct extensive numerical experiments to corroborate our
theoretical findings. First, we study the numerical convergence rate of the
FMGDM method. Second, we explore the statistical properties of the stable
solutions corresponding to the different tuning parameters specifications.
Last, we compare the performances of the different types of MGDM methods.

3.1. Numerical Convergence

In this subsection, we demonstrate how the tuning parameters $\alpha$ and
$\gamma$ affect the numerical convergence speed of the FMGDM method.
Specifically, we consider a standard linear regression model (2.1) with
$N=5000$ and $p=50$. The random noise $\varepsilon_{i}$’s are simulated from
the standard normal distribution with mean $0$ and variance $\sigma^{2}=1$.
The true regression coefficient $\theta_{0}$ is given by
$\theta_{0}=(\theta_{01},\dots,\theta_{0p})^{\top}$ with
$\theta_{0j}=10\exp(-0.5j)$ for $1\leq j\leq p$. Furthermore, the covariate
$X_{i}$’s are simulated from a multivariate normal distribution
$N_{p}(\boldsymbol{0},\Sigma_{xx})$, where
$\Sigma_{xx}=I_{p}+\kappa\mbox{\boldmath$1$}_{p}\mbox{\boldmath$1$}_{p}^{\top}$.
One can easily verify that $\lambda_{1}=\lambda_{\max}(\Sigma_{xx})=\kappa
p+1$ and $\lambda_{p}=\lambda_{\min}(\Sigma_{xx})=1$. Hence the condition
number of $\Sigma_{xx}$ is given by $\kappa p+1$. To compare the convergence
speed of the fixed minibatch-based GD (FMGD, i.e., $\gamma=0$) and FMGDM
methods, we consider two different specifications: (i) (FMGD)
$\alpha=2/(\lambda_{1}+\lambda_{p})$ and $\gamma=0$ and (ii) (FMGDM)
$\alpha=4/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$ and
$\gamma=(\sqrt{\lambda_{1}}-\sqrt{\lambda_{p}})^{2}/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$.
As indicated by Theorem 3, these two specifications correspond to the optimal
tuning parameters configurations of the FMGD and FMGDM methods, respectively.
Different $\kappa$ values are studied (i.e., $\kappa=1,2,3$), and the
minibatch size is fixed to $n=500$. Accordingly, the number of minibatches is
given by $M=N/n=10$. For each setting, we first compute the corresponding
stable solution on the $M$-th minibatch $\widehat{\theta}^{(M)}$ by solving
the linear dynamical system (2.6) directly. Next, we study how the algorithm
converges to it. To this end, we run the algorithm for a total of $T=30$
epochs. Recall that $\widehat{\theta}^{(t,M)}$ is the estimate obtained from
the $M$-th minibatch in the $t$-th epoch. We then record the difference
$\delta^{(t)}=\|\widehat{\theta}^{(t,M)}-\widehat{\theta}^{(M)}\|$ for each
epoch $1\leq t\leq T$. We randomly replicate the experiments a total number of
$100$ times. This leads to a total of $100$ $\delta^{(t)}$ values for each
epoch $1\leq t\leq T$. The medians of these values are then plotted in Figure
1.

Figure 1: The numerical convergence rate of different algorithms (i.e., the
FMGD and FMGDM methods). The vertical axis corresponds to the $\delta^{(t)}$
value, while the horizontal axis corresponds to the epoch. The red solid line
stands for the FMGDM method and the blue dashed line stands for the FMGD
method.

From Figure 1, we can obtain the following interesting observations. First, we
find that the curve of the FMGDM method converges towards $0$ much faster than
that of the FMGD method. This result confirms that our theoretical claims made
in Theorem 3. Second, we can see that the convergence rates of both methods
become slower as the conditional number increases. Moreover, the FMGDM method
exhibits a more substantial advantage over the FMGD method. These findings
also corroborate our theoretical results in Theorem 3 (a) and (b) very well.

3.2. Statistical Efficiency of the Stable Solution

This subsection is devoted to verifying the theoretical claims made in Section
2.5, that is, how the two tuning parameters (i.e., $\alpha$ and $\gamma$)
affect the statistical efficiency of the stable solution. Specifically, we
consider two different simulation experiments. One experiment studies the
effects of the learning rate $\alpha$ with a fixed momentum parameter
$\gamma$. The other experiment examines the effect of the momentum parameter
$\gamma$ while holding the learning rate $\alpha$ fixed. In the first case, we
fix $\gamma=0.5$ and let $\alpha=0.15,0.1,0.05,0.01,0.001$, while in the
second case, we fix $\alpha=0.1$ and let $\gamma=0,0.5,0.8,1,2,5,10$.
Subsequently, we generate the simulated dataset as in Section 3.1 with
$\kappa=1$. For each given $(\alpha,\gamma)$-specification, the stable
solution is computed for the last minibatch as $\widehat{\theta}^{(M)}$ by
solving the linear dynamical system (2.6) directly. Then we evaluate the
efficiency of different estimators by computing the estimation error (EE) as
$\|\widehat{\theta}^{(M)}-\theta_{0}\|^{2}$, where $\theta_{0}$ is the true
regression coefficient. For comparison, we also compute the EE values of the
OLS estimator. This experiment is randomly replicated a total of $100$ times,
leading to a total of $100$ EE values for each estimator, which are then log-
transformed and shown as boxplots in Figure 2.

Figure 2: Boxplots of $\log(\operatorname{EE})$ of the FMGDM method under
different $(\alpha,\gamma)$-specifications. In the left panel, we fix the
momentum parameter as $\gamma=0.5$ and let learning rate $\alpha$ range from
$0.15$ to $0.001$. In the right panel, we fix the learning rate as
$\alpha=0.1$ and let the momentum parameter $\gamma$ range from $0$ to $10$.

From Figure 2, we can draw the following conclusions. First, from the left
panel, we can see that, when the learning rate is large (e.g., $0.2,0.1$), the
EE values of corresponding stable solutions are much larger than that of the
global OLS estimator. As the learning rate $\alpha$ decays, the EE values
steadily drop and then converge to that of the OLS estimator. These findings
verify the results in Theorem 4 (a) and Theorem 5. That is, by letting
$\alpha\to 0$, the stable solution should have the same statistical efficiency
as the global OLS estimator. Second, from the right panel of Figure 2 we can
see that when $\gamma$ is small (e.g., $\gamma\leq 2$), the EE values of the
corresponding stable solutions are much larger than that of the OLS estimator.
As $\gamma$ increases, the EE values finally converge to that of the OLS
estimator. This finding is in agreement with the claim in Theorem 4 (b).
Finally, we find that the EE values of the stable solution seem to increase
from $\gamma=0$ to $\gamma=0.8$. To investigate this phenomenon, recall that
by Theorem 5, the extra variance term $\alpha^{2}Md_{\gamma}^{2}\Sigma_{xx}$
should be the dominant term of the EE value. Here, $d_{\gamma}$ is the only
factor involving $\gamma$. Therefore, we calculate values of $d_{\gamma}$ for
$\gamma$ ranging from $0$ to $10$ by the explicit formulas in Appendix A.6.
These values are then plotted, as shown in Figure 3. It can be easily seen
that the curve of $d_{\gamma}$ has a similar pattern as that of the EE values
in the right panel of Figure 2. That is, the value of $d_{\gamma}$ first
increases with increasing $\gamma$ from $0$ to some value slightly smaller
than $1$, and then, $d_{\gamma}$ decreases as $\gamma$ increases. This result
further validates the asymptotic variance formula given in Theorem 5.

Figure 3: The curve of $d_{\gamma}$ with $\gamma$ ranging from $0$ to $10$.
The values of $d_{\gamma}$ are calculated by the explicit formulas (A.10) and
(A.12) in Appendix A.6.

3.3. Different Types of MGDM Methods

Next, we compare the performance of the different types of MGDM methods.
Specifically, we consider three different types of minibatch methods, i.e.,
fixed, shuffled and random minibatches. The corresponding MGDM methods are
referred to as FMGDM, SMGDM and RMGDM methods, respectively. For the RMGDM
method, we use simple random sampling without replacement to generate each
minibatch. In this experiment, we consider three different learning rates as
$\alpha=0.04,0.02,0.01$ and fix the momentum parameter as $\gamma=0.9$ by
convention. Subsequently, we simulate the dataset as in Section 3.1 with
$\kappa=1$. Once the data are generated, we run each MGDM algorithm for a
total of $T=50$ epochs and record the resulting estimator
$\widehat{\theta}^{(T,M)}$ from the $M$-th minibatch. For a fair comparison,
passing over every $M=N/n=10$ minibatches is considered as one epoch for each
MGDM method. Then we compute the corresponding estimation error (EE) as
$\|\widehat{\theta}^{(T,M)}-\theta_{0}\|^{2}$ for each method under each
$(\alpha,\gamma)$-specification. For comparison, we also compute the EE values
of the OLS estimator. We randomly replicate the experiment a total of $100$
times, leading to a total of $100$ EE values for each estimator. These EE
values are then log-transformed and shown as boxplots in Figure 4.

Figure 4: Boxplots of $\log(\operatorname{EE})$ of the FMGDM, SMGDM, RMGDM and
OLS estimators. For the three MGDM methods, we consider three different
learning rates $\alpha=0.04,0.02,0.01$ and fix the momentum parameter as
$\gamma=0.9$. Red boxes represent the FMGDM estimator, blue boxes represent
the SMGDM estimator, green boxes represent the RMGDM estimator, and the black
box represents the OLS estimator.

From Figure 4, we immediately observe that the FMGDM estimator has the
smallest EE values among the three MGDM methods. For the SMGDM method, we can
find that the corresponding EE values converge to that of the OLS estimator as
the learning rate $\alpha$ decreases. This finding corroborates our
justifications in Section 2.6 very well. That is, the SMGDM estimator should
enjoy the same statistical efficiency as the OLS estimator with a learning
rate of $\alpha\to 0$. The EE values of the RMGDM estimator also decrease with
decreasing learning rate. However, there is still a clear gap between the
RMGDM estimator and the OLS estimator in terms of EE values, even though the
learning rate is sufficiently small (i.e., $\alpha=0.01$). This result may be
due to the extra random errors created by simple random sampling. These
results demonstrate that the FMGDM and SMGDM methods are more preferable in
terms of statistical efficiency compared to the RMGDM method.

3.4. Application to Airline Data

To demonstrate the performance of the MGDM methods on a real-world dataset, we
consider here the U.S. Airline Dataset. The whole dataset is available for
download at http://stat-computing.org/dataexpo/2009. It contains the flight
arrival and departure details for all commercial flights within the USA from
1987 to 2008. We use the data from 2008, and aim to predict the delay time in
the arrival of a flight given other flight information. Each record of the
data contains the arrival delay, departure time, scheduled arrival and
departure time, flight date, distance of the flight, information of the
carrier, origin and destination. The detailed variable information is
described in Table 1. The six continuous variables are standardized to be mean
0 and variance 1, and the five categorical variables are converted to dummy
variables with appropriate dimensions. Finally, a total of 126 predictors and
one response variable are used in the linear regression model, and the total
sample size is over 2.5 million observations.

In this experiment, we set the minibatch size as $n=1000$, momentum parameter
as $\gamma=0.9$ and consider three different learning rates
$\alpha=0.4,0.2,0.1$. For comparison purpose, we first randomly choose a
subsample of size $N=10^{5}$ from the whole sample, and then run each MGDM
algorithm with $M=N/n=100$ for a total of $T=100$ epochs. Next, we compute the
estimation error (EE) as
$\|\widehat{\theta}^{(T,M)}-\widehat{\theta}_{\text{ols}}\|^{2}$ for each
method under each $(\alpha,\gamma)$-specification, where
$\widehat{\theta}_{\text{ols}}$ is the OLS estimator of the regression
coefficient based on the whole sample. We replicate the experiment through
$50$ times random subsampling, leading to a total of $50$ EE values for each
estimator. These EE values are then log-transformed and shown as boxplots in
Figure 5. From Figure 5, we can see a similar pattern as shown in Figure 4.
That is, as the learning rate $\alpha$ decreases, the EE values of all three
methods decrease. In addition, the FMGDM and SMGDM methods show a competitive
performance. The RMGDM method has the worst performance under all three
$\alpha$ specifications.

Figure 5: Boxplots of $\log(\operatorname{EE})$ of the FMGDM, SMGDM, and RMGDM
estimators based on the airline data. For the three MGDM methods, we consider
three different learning rates $\alpha=0.4,0.2,0.1$ and fix the momentum
parameter as $\gamma=0.9$. Red boxes represent the FMGDM estimator, blue boxes
represent the SMGDM estimator, green boxes represent the RMGDM estimator.

Table 1: Variable description for the U.S. airline data from 2008. Numerical variables are standardized to have mean 0 and variance 1. Variable | Description | Variable used in the model
---|---|---
`ArrDelay` | How long was the delay in the arrival of the flight | Used as the numerical response variable
`Month` | Which month of the year | Converted to 11 dummies
`DayofMonth` | Which day of the month | Used as numerical variable
`DayOfWeek` | Which day of the week | Converted to 6 dummies
`DepTime` | Actual departure time | Used as numerical variable
`CRSDepTime` | Scheduled departure time | Used as numerical variable
`CRSArrTime` | Scheduled arrival time | Used as numerical variable
`Distance` | Distance between the origin and destination in miles | Used as numerical variable
`UniqueCarrier` | Flight carrier code for 29 carriers | Top 7 carries converted to 6 dummies
`Origin` | Departing origin | Top 50 origin cities converted to 49 dummies
`Dest` | Destination of the flight | Top 50 origin cities converted to 49 dummies

4\. CONCLUDING REMARKS

In this paper, we study the theoretical properties of MGDM algorithms for
linear regression models. We focus mainly on two different types of MGDM
methods, i.e., fixed MGDM (FMGDM) and shuffled MGDM (SMGDM) methods. We first
investigate the FMGDM method by introducing a linear dynamical system. Then,
we provide nearly sufficient and necessary conditions for numerical
convergence. Furthermore, the numerical convergence rate of FMGDM algorithm is
investigated. This analysis leads to the theoretically optimal specification
of the tuning parameters. Theoretical findings reveal that the extra momentum
term could greatly speed up numerical convergence compared to a standard
gradient descent method without momentum. This is particularly true for
problems with large condition numbers. Subsequently, the statistical
properties of the resulting FMGDM estimator are carefully investigated and the
asymptotic normality is further established. We then find that with a
sufficiently small learning rate, the MGDM estimators (both fixed and
shuffled) can share the same statistical efficiency as the global estimator.
Finally, these theoretical results are fully validated through extensive
numerical studies. To conclude this article, we discuss here several
interesting topics for future study. First, we examine the FMGDM algorithm for
the quadratic objective function, due to its analytical simplicity. Of note,
many learning problems have non-quadratic objection functions, such as
logistic regression and quantile regression models. However, it seems to us
that the theoretical properties of either the fixed or shuffled MGDM
estimators for these problems remain unknown. This is a research topic of
great interest. Second, modern statistical inference often involves ultrahigh-
dimensional problems. It is valuable to extend our results on the MGDM methods
to the models with diverging $p$. Last, there are many momentum-related
methods (e.g., Nesterov’s momentum method, RMSprop and ADAM) that are widely
used in practice. Therefore, it is also very important to investigate their
theoretical properties under a fixed or shuffled minibatch setup.

APPENDIX A

Appendix A.1 Proof of Theorem 1

Construct a block matrix analogous to $\widehat{\Omega}$ in (2.7) as

$\displaystyle\Omega=\begin{bmatrix}I_{p}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\gamma
I_{p}&-\Delta\\\
-\Delta&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\gamma
I_{p}\\\ \gamma
I_{p}&-\Delta&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\gamma
I_{p}&-\Delta&I_{p}\end{bmatrix}\in\mathbb{R}^{q\times q},$ (A.1)

where $\Delta=(1+\gamma)I_{p}-\alpha\Sigma_{xx}\in\mathbb{R}^{p\times p}$. By
Lemma 1 in the Appendix B, we know that $\Omega$ is invertible. Note that

$\displaystyle\widehat{\Omega}-\Omega=\alpha\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(1)}-\Sigma_{xx}\\\
\widehat{\Sigma}_{xx}^{(2)}-\Sigma_{xx}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(3)}-\Sigma_{xx}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(M)}-\Sigma_{xx}&\boldsymbol{0}\\\
\end{bmatrix}.$

By law of large number, we should have
$\|\widehat{\Sigma}_{xx}^{(m)}-\Sigma_{xx}\|\to_{p}0$ as $n\to\infty$ for each
$1\leq m\leq M$. Hence, $\|\widehat{\Omega}-\Omega\|\to_{p}0$ as $n\to\infty$.
Since determinant is a continuous function, we should have
$|\widehat{\Omega}|\to_{p}|\Omega|\neq 0$ as $n\to\infty$. Thus, the claim of
the theorem follows immediately.

Appendix A.2 Proof of Theorem 2

By (2.4) and (2.5), we have

$\displaystyle\widehat{\theta}^{(t,m)}$
$\displaystyle=\widehat{\theta}^{(t,m-1)}-\gamma(\widehat{\theta}^{(t,m-2)}-\widehat{\theta}^{(t,m)})-\alpha\Big{\\{}\widehat{\Sigma}_{xx}^{(m)}\widehat{\theta}^{(t,m-1)}-\widehat{\Sigma}_{xy}^{(m)}\Big{\\}}$
$\displaystyle=\Big{\\{}(1+\gamma)I_{p}-\alpha\widehat{\Sigma}_{xx}^{(m)}\Big{\\}}\widehat{\theta}^{(t,m-1)}-\gamma\widehat{\theta}^{(t,m-2)}+\alpha\widehat{\Sigma}_{xy}^{(m)}.$

Note the stable solution also satisfies
$\widehat{\theta}^{(m)}=\Big{\\{}(1+\gamma)I_{p}-\alpha\widehat{\Sigma}_{xx}^{(m)}\Big{\\}}\widehat{\theta}^{(m-1)}-\gamma\widehat{\theta}^{(m-2)}+\alpha\widehat{\Sigma}_{xy}^{(m)}$.
Hence, taking a difference of above two equations gives

$\displaystyle\widehat{\theta}^{(t,m)}-\widehat{\theta}^{(m)}$
$\displaystyle=\Big{\\{}(1+\gamma)I_{p}-\alpha\widehat{\Sigma}_{xx}^{(m)}\Big{\\}}\Big{(}\widehat{\theta}^{(t,m-1)}-\widehat{\theta}^{(m-1)}\Big{)}-\gamma\Big{(}\widehat{\theta}^{(t,m-2)}-\widehat{\theta}^{(m-2)}\Big{)},$

where $\widehat{\theta}^{(t,0)}=\widehat{\theta}^{(t-1,M)}$ and
$\widehat{\theta}^{(t,-1)}=\widehat{\theta}^{(t-1,M-1)}$. If we denote
$\delta^{(t,m)}=\widehat{\theta}^{(t,m)}-\widehat{\theta}^{(m)}$, then we have

$\displaystyle r^{(t,m)}=\begin{bmatrix}\delta^{(t,m)}\\\
\delta^{(t,m-1)}\end{bmatrix}=\begin{bmatrix}\widehat{\Delta}^{(m)}&-\gamma
I_{p}\\\ I_{p}&\boldsymbol{0}\end{bmatrix}\begin{bmatrix}\delta^{(t,m-1)}\\\
\delta^{(t,m-2)}\end{bmatrix}=D^{(m)}r^{(t,m-1)}.$

Consequently, we obtain the recursion formula as
$r^{(t,m)}=C^{(m)}r^{(t-1,m)}$, where

$\displaystyle C^{(1)}=D^{(1)}D^{(M)}D^{(M-1)}\cdots D^{(2)},$ $\displaystyle
C^{(m)}=D^{(m)}\cdots D^{(1)}D^{(M)}\cdots D^{(m+1)},\ 2\leq m\leq M-1,$
$\displaystyle C^{(M)}=D^{(M)}D^{(M-1)}\cdots\cdots D^{(1)}.$

Proof of (a). Consider the block matrix $D=\big{[}\Delta,-\gamma
I_{p};I_{p},\boldsymbol{0}\big{]}\in\mathbb{R}^{(2p)\times(2p)}$, where
$\Delta=(1+\gamma)I_{p}-\alpha\Sigma_{xx}\in\mathbb{R}^{p\times p}$. By
assumed conditions in (a) and Lemma 2(a), we know that $\rho(D)<1$. Note that
$\|D^{(m)}-D\|=\alpha\|\widehat{\Sigma}_{xx}^{(m)}-\Sigma_{xx}\|\to_{p}0$ as
$n\to\infty$. Since $M$ is fixed, it follows immediately that
$\|C^{(m)}-D^{M}\|\to_{p}0$ as $n\to\infty$. Because the spectral radius is a
continuous function of the entries of a matrix, we obtain that
$\rho(C^{(m)})\to_{p}\rho(D^{M})=\rho^{M}(D)<1$ as $n\to\infty$. Consequently,
we should have $\rho\\{(C^{(m)})^{t}\\}=\rho^{t}(C^{(m)})\to_{p}0$ as
$n\to\infty$. Then by recursion formula, we have
$r^{(t,m)}=\\{C^{(m)}\\}^{t}r^{(0,m)}\to_{p}\boldsymbol{0}$ as $n\to\infty$.
Thus, (a) follows immediately.

Proof of (b). By assumed conditions in (b) and Lemma 2(e), we have
$\rho(D)>1$. Then using the similar procedures in the proof of (a), we can
obtain that $\rho(C^{(m)})\to_{p}\rho^{M}(D)>1$ as $n\to\infty$. Consequently,
$\\{C^{(m)}\\}^{t}\not\to\boldsymbol{0}$ as $t\to\infty$ holds with
probability tending to one as $n\to\infty$ (Horn and Johnson,, 2012, Theorem
5.6.12). Thus, as long as $\delta^{(0,m)}\neq\boldsymbol{0}$, $\delta^{(t,m)}$
will not converge to zero with probability tending to one as $n\to\infty$. We
have proved (b) and complete the proof.

Appendix A.3 Proof of Theorem 3

From the proof of Theorem 2, we know that $\gamma\leq\rho(D)<1$ if
$0\leq\gamma<1$ and $0<\alpha<2(1+\gamma)/\lambda_{1}$, and
$\rho(C^{(m)})=\rho^{M}(D)+\zeta_{n}$, where $\zeta_{n}\to_{p}0$ as
$n\to\infty$. Gelfand formula (Horn and Johnson,, 2012, Corollary 5.6.14)
guarantees that $\lim_{t\to\infty}\|\\{C^{(m)}\\}^{t}\|^{1/t}=\rho(C^{(m)})$.
Then we have $\|\\{C^{(m)}\\}^{t}\|=(\rho(C^{(m)})+\eta_{t})^{t}$, where
$\eta_{t}\to 0$ as $t\to\infty$. Recall that
$r^{(t,m)}=\\{C^{(m)}\\}^{t}r^{(0,m)}$. Consequently,
$\|\delta^{(t,m)}\|\leq\|r^{(t,m)}\|\leq\|\\{C^{(m)}\\}^{t}\|\
\|r^{(0,m)}\|\leq\\{\rho^{M}(D)+\zeta_{n}+\eta_{t}\\}^{t}\big{(}\|\delta^{(0,m)}\|+\|\delta^{(0,m-1)}\|\big{)}$.
Then the conclusions (a), (b) and (c) follow directly from Lemma 2(c), (b) and
(d) in Appendix B, respectively.

Appendix A.4 Proof of Theorem 4

Recall that $P=[P_{1},P_{2}]$ is a orthogonal matrix, where
$P_{1}=I^{*}/\sqrt{M}\in\mathbb{R}^{q\times p}$. In fact, we can let
$P_{2}=Z\otimes I_{p}\in\mathbb{R}^{q\times(q-p)}$, where

$\displaystyle
Z=\begin{bmatrix}1/\sqrt{2}&1/\sqrt{6}&\cdots&1/\sqrt{M(M-1)}\\\
-1/\sqrt{2}&1/\sqrt{6}&\cdots&1/\sqrt{M(M-1)}\\\
0,&-2/\sqrt{6}&\cdots&1/\sqrt{M(M-1)}\\\ \vdots&\vdots&\vdots&\vdots\\\
0&0&\cdots&-(M-1)/\sqrt{M(M-1)}\end{bmatrix}\in\mathbb{R}^{M\times(M-1)}.$
(A.2)

Proof of (a). We will prove (a) in two steps. In the first step, we obtain the
expression of $\widehat{\Omega}^{*}=(P^{\top}\widehat{\Omega}P)^{-1}$. In the
second step, we investigate stable solution by $\widehat{\theta}^{*}=\alpha
P(P^{\top}\widehat{\Omega}P)^{-1}P^{\top}\widehat{\Sigma}_{xy}^{*}$.

Step 1. For the simplicity, we further denote
$\widehat{\Omega}^{*}=[\widehat{\Omega}_{11}^{*},\widehat{\Omega}_{12}^{*};\widehat{\Omega}_{21}^{*},\widehat{\Omega}_{22}^{*}]$,
where $\widehat{\Omega}_{11}^{*}\in\mathbb{R}^{p\times p}$,
$\widehat{\Omega}_{12}^{*}\in\mathbb{R}^{p\times(q-p)}$,
$\widehat{\Omega}_{21}^{*}\in\mathbb{R}^{(q-p)\times p}$, and
$\widehat{\Omega}_{22}^{*}\in\mathbb{R}^{(q-p)\times(q-p)}$. We then compute
each block of $\widehat{\Omega}^{*}$ subsequently. By techniques used in the
proof of Lemma 1, one can verify that $A$ is diagonalizable with $q-p$ non-
zero eigenvalues. Hence,
$\mbox{rank}(A)=\mbox{rank}(P^{\top}AP)=\mbox{rank}(A_{22}^{p})=q-p$. i.e.,
$A_{22}^{p}$ is invertible. Note that
$\widehat{B}_{11}^{P}=P_{1}^{\top}\widehat{B}P_{1}=M^{-1}\sum_{m=1}^{M}\widehat{\Sigma}_{xx}^{(m)}=\widehat{\Sigma}_{xx}$
is also invertible almost surely. Consequently, it follows from the formula
for the inverse of a block matrix (Seber,, 2008) and Taylor expansion that

$\displaystyle\widehat{\Omega}_{11}^{*}$
$\displaystyle=\Big{[}\alpha\widehat{B}_{11}^{P}-(\alpha\widehat{B}_{12}^{P})(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}(\alpha\widehat{B}_{21}^{P})\Big{]}^{-1}$
$\displaystyle=\alpha^{-1}\Big{[}\widehat{\Sigma}_{xx}-\alpha\widehat{B}_{12}^{P}\Big{\\{}(A_{22}^{P})^{-1}-(A_{22}^{P})^{-1}\widehat{B}_{22}^{P}(A_{22}^{P})^{-1}+O_{p}(\alpha^{2})\Big{\\}}\widehat{B}_{21}^{P}\Big{]}^{-1}$
$\displaystyle=\alpha^{-1}\widehat{\Sigma}_{xx}^{-1}+\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}-\alpha\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}\widehat{B}_{22}^{P}(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}+O_{p}(\alpha^{2}).$

Similarly, it can be derived that

$\displaystyle\widehat{\Omega}_{12}^{*}=$
$\displaystyle-\widehat{\Omega}_{11}^{*}(\alpha\widehat{B}_{12}^{P})(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}$
$\displaystyle=$
$\displaystyle-\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}+\alpha\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}\widehat{B}_{22}^{P}(A_{22}^{P})^{-1}$
$\displaystyle-\alpha\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}+O_{p}(\alpha^{2}),$
$\displaystyle\widehat{\Omega}_{21}^{*}=$
$\displaystyle-(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}(\alpha\widehat{B}_{21}^{P})\widehat{\Omega}_{11}^{*}$
$\displaystyle=$
$\displaystyle-(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}+\alpha(A_{22}^{P})^{-1}\widehat{B}_{22}^{P}(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}$
$\displaystyle-\alpha(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}+O_{p}(\alpha^{2}),$
$\displaystyle\widehat{\Omega}_{22}^{*}=$
$\displaystyle(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}+(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}(\alpha\widehat{B}_{21}^{P})\widehat{\Omega}_{11}^{*}(\alpha\widehat{B}_{12}^{P})(A_{22}^{P}+\alpha\widehat{B}_{22}^{P})^{-1}$
$\displaystyle=$
$\displaystyle(A_{22}^{P})^{-1}-\alpha(A_{22}^{P})^{-1}\widehat{B}_{22}^{P}(A_{22}^{P})^{-1}+\alpha(A_{22}^{P})^{-1}\widehat{B}_{21}^{P}\widehat{\Sigma}_{xx}^{-1}\widehat{B}_{12}^{P}(A_{22}^{P})^{-1}+O_{p}(\alpha^{2}).$

Denote $\widehat{H}_{1}=P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}$ and
$H_{2}=P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}$. Then we should have

$\displaystyle\widehat{\Omega}^{-1}=$ $\displaystyle
P\widehat{\Omega}^{*}P^{\top}=P_{1}\widehat{\Omega}_{11}^{*}P_{1}^{\top}+P_{1}\widehat{\Omega}_{12}^{*}P_{2}^{\top}+P_{2}\widehat{\Omega}_{21}^{*}P_{1}^{\top}+P_{2}\widehat{\Omega}_{22}^{*}P_{2}^{\top}$
$\displaystyle=$
$\displaystyle\alpha^{-1}\widehat{H}_{1}+(\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}-\widehat{H}_{1}\widehat{B}H_{2}-H_{2}\widehat{B}\widehat{H}_{1}+H_{2})$
$\displaystyle-\alpha(\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}-\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}H_{2}-H_{2}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}+H_{2}\widehat{B}H_{2})$
$\displaystyle+\alpha(H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}-\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}-H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1})+O_{p}(\alpha^{2})$
$\displaystyle=$
$\displaystyle\alpha^{-1}\widehat{H}_{1}+(\widehat{H}_{1}\widehat{B}-I_{q})H_{2}(\widehat{B}\widehat{H}_{1}-I_{q})-\alpha(\widehat{H}_{1}\widehat{B}-I_{q})(H_{2}BH_{2})(\widehat{B}\widehat{H}_{1}-I_{q})$
$\displaystyle+\alpha(\widehat{H}_{1}\widehat{B}-I_{q})(H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2})(\widehat{B}\widehat{H}_{1}-I_{q})-\alpha\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}+O_{p}(\alpha^{2}).$
(A.3)

Step 2. Recall that $P_{1}=I^{*}/\sqrt{M}$. Hence,
$P_{1}^{\top}\widehat{\Sigma}_{xy}^{*}=M^{-1/2}\sum_{m=1}^{M}\widehat{\Sigma}_{xy}^{(m)}=\sqrt{M}\widehat{\Sigma}_{xy}$.
Thus,
$\widehat{H}_{1}\widehat{\Sigma}_{xy}^{*}=\sqrt{M}P_{1}\widehat{\Sigma}_{xx}^{-1}\widehat{\Sigma}_{xy}=I^{*}\widehat{\theta}_{\text{ols}}=\widehat{\theta}_{\text{ols}}^{*}$.
Combine the results above, it can be concluded that

$\widehat{\theta}^{*}=\alpha\widehat{\Omega}^{-1}\widehat{\Sigma}_{xy}^{*}=\widehat{\theta}_{\text{ols}}^{*}+\alpha\widehat{E}\big{\\{}1+O_{p}(\alpha)\big{\\}},$
(A.4)

where

$\displaystyle\widehat{E}=$
$\displaystyle(\widehat{H}_{1}\widehat{B}-I_{q})H_{2}(\widehat{B}\widehat{H}_{1}-I_{q})\widehat{\Sigma}_{xy}^{*}$
$\displaystyle=$
$\displaystyle\big{(}P_{1}\widehat{\Sigma}_{xx}^{-1}P_{1}^{\top}\widehat{B}-I_{q}\big{)}\big{\\{}P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}\big{\\}}\big{(}\widehat{B}\widehat{\theta}_{\text{ols}}^{*}-\widehat{\Sigma}_{xy}^{*}\big{)}.$

We finish the proof of (a).

Proof of (b). From the proof of (a), we can find that $\gamma$ affects
$\widehat{\Omega}^{*}$ only through $A_{22}^{P}$. Note that
$I_{q}=PP^{\top}=P_{1}P_{1}^{\top}+P_{2}P_{2}^{\top}$. Then we have

$\displaystyle\|A_{22}^{P}\|_{F}^{2}=\mbox{tr}\Big{\\{}A_{22}^{P}(A_{22}^{P})^{\top}\Big{\\}}$
$\displaystyle=\mbox{tr}(P_{2}^{\top}AP_{2}P_{2}^{\top}A^{\top}P_{2})=\mbox{tr}\Big{\\{}A(I_{q}-P_{1}P_{1}^{\top})A^{\top}(I_{q}-P_{1}P_{1}^{\top})\Big{\\}}$
$\displaystyle=\mbox{tr}(AA^{\top})-2\mbox{tr}(AP_{1}P_{1}^{\top}A^{\top})+\mbox{tr}(AP_{1}P_{1}^{\top}AP_{1}P_{1}^{\top})$
$\displaystyle=\mbox{tr}(AA^{\top})=\big{\\{}(1+\gamma)^{2}+\gamma^{2}+1\big{\\}}q,$

where we have used the fact $AP_{1}=\boldsymbol{0}$. Hence, we conclude that
$\|(A_{22}^{P})^{-1}\|\to 0$ as $\gamma\to\infty$. Then by analogous
procedures to derive (A.4), we can obtain that
$\widehat{\theta}^{*}=\widehat{\theta}_{\text{ols}}^{*}+\alpha\widehat{E}\Big{\\{}1+O_{p}(\|(A_{22}^{P})^{-1}\|)\Big{\\}}$,
where $\widehat{E}$ is the same as that in (A.4). The only difference is that
the remainder term is $O_{p}(\|(A_{22}^{P})^{-1}\|)$ instead of
$O_{p}(\alpha)$. We have proved (b) and complete the proof.

Appendix A.5 Proof of Theorem 5

Note that
$\widehat{\Sigma}_{xx}^{(m)}=n^{-1}\sum_{i\in\mathcal{S}_{(m)}}X_{i}X_{i}^{\top}$
and $\widehat{\Sigma}_{xy}^{(m)}=n^{-1}\sum_{i\in\mathcal{S}_{(m)}}X_{i}Y_{i}$
are sample mean of $n$ independent random variables. Then by central limit
theorem, we should have
$\sqrt{n}(\widehat{\Sigma}_{xy}^{(m)}-\widehat{\Sigma}_{xx}^{(m)}\theta_{0})\to_{d}N(\boldsymbol{0},\sigma^{2}\Sigma_{xx})$
as $N\to\infty$ for each $1\leq m\leq M$. Write the block diagonal matrix as
$\widehat{\Sigma}_{xx}^{*}=\mbox{diag}\\{\widehat{\Sigma}_{xx}^{(1)},\dots,\widehat{\Sigma}_{xx}^{(M)}\\}\in\mathbb{R}^{q\times
q}$. Then it follows that
$\sqrt{n}\big{\\{}\widehat{\Sigma}_{xy}^{*}-\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})\big{\\}}\to_{d}N_{q}(\boldsymbol{0},\sigma^{2}I_{M}\otimes\Sigma_{xx})$,
where
$\widehat{\Sigma}_{xy}^{*}=\big{(}\widehat{\Sigma}_{xy}^{(1)\top},\dots,\widehat{\Sigma}_{xy}^{(M)\top}\big{)}^{\top}\in\mathbb{R}^{q}$.
Recall that the stable solution has the form
$\widehat{\theta}^{*}=\alpha\widehat{\Omega}^{-1}\widehat{\Sigma}_{xy}^{*}$.
Consequently,
$\sqrt{n}\big{\\{}\widehat{\theta}^{*}-\widehat{\mu}(\alpha)\big{\\}}\to_{d}N_{q}\big{(}\boldsymbol{0},\sigma^{2}V(\alpha)\big{)}$,
where
$\widehat{\mu}(\alpha)=\alpha\widehat{\Omega}^{-1}\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})$,
and $V(\alpha)=\alpha^{2}\Omega^{-1}(I_{M}\otimes\Sigma_{xx})\Omega^{-1}$.
Here, $\Omega$ is defined in (A.1), and it is the probabilistic limit of
$\widehat{\Omega}$, i.e., $\widehat{\Omega}\to_{p}\Omega$ as $N\to\infty$. In
the following steps, we investigate $\widehat{\mu}(\alpha)$ and $V(\alpha)$
subsequently.

Step 1. By the asymptotic form of $\widehat{\Omega}^{-1}$ in (A.3), we can
show that
$\widehat{\mu}(\alpha)=\mbox{\boldmath$1$}_{M}\otimes\theta_{0}-\alpha^{2}\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})+O_{p}(\alpha^{2})$.
This is due to the fact that
$\widehat{H}_{1}\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})=\mbox{\boldmath$1$}_{M}\otimes\theta_{0}$
and
$\widehat{B}\widehat{H}_{1}\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})=\widehat{\Sigma}_{xx}^{*}(\mbox{\boldmath$1$}_{M}\otimes\theta_{0})$.
Next, we can rewrite $\widehat{B}$ as
$\widehat{B}=J^{M-1}\otimes\widehat{\Sigma}_{xx}+\widehat{R}$, where
$J\in\mathbb{R}^{M\times M}$ is the permutation matrix defined in (A.11), and

$\displaystyle\widehat{R}=\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(1)}-\widehat{\Sigma}_{xx}\\\
\widehat{\Sigma}_{xx}^{(2)}-\widehat{\Sigma}_{xx}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(3)}-\widehat{\Sigma}_{xx}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\widehat{\Sigma}_{xx}^{(M)}-\widehat{\Sigma}_{xx}&\boldsymbol{0}\\\
\end{bmatrix}\in\mathbb{R}^{q\times q}.$

By central limit theorem, it can be shown that
$\|\widehat{R}\|=O_{p}(n^{-1/2})$. As derived in Appendix A.6, we have
$H_{2}=Q\otimes I_{p}$, where
$Q=(a,J^{M-1}a,\dots,Ja)^{\top}\in\mathbb{R}^{M\times M}$ for some vector
$a\in\mathbb{R}^{M}$ satisfying
$\mbox{\boldmath$1$}_{M}^{\top}a=\boldsymbol{0}$. Then we can verify that
$\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}_{M}^{\top}J^{M-1}Q=\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}_{M}^{\top}Q=\boldsymbol{0}$.
Consequently, it follows that
$\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}=\widehat{H}_{1}\widehat{R}H_{2}\widehat{R}\widehat{H}_{1}\widehat{R}H_{2}$.
Thus, we should have
$\|\widehat{H}_{1}\widehat{B}H_{2}\widehat{B}\widehat{H}_{1}\widehat{B}H_{2}\|=O_{p}(n^{-3/2})$.
Finally, we obtain that
$\widehat{\mu}(\alpha)=\mbox{\boldmath$1$}_{M}\otimes\theta_{0}+O_{p}(\alpha^{2}n^{-3/2})$.

Step 2. By law of large numbers, we have $\Sigma_{xx}\to_{p}\Sigma_{xx}$ and
$\Sigma_{xx}^{(m)}\to_{p}\Sigma_{xx}$ as $N\to\infty$ for each $1\leq m\leq
M$. Then we can obtain $\widehat{B}\to_{p}B$ and $\widehat{H}_{1}\to_{p}H_{1}$
as $N\to\infty$, where $B=J^{M-1}\otimes\Sigma_{xx}$ and
$H_{1}=P_{1}\Sigma_{xx}^{-1}P_{1}^{\top}$. By replacing $\widehat{B}$ and
$\widehat{H}_{1}$ by $B$ and $H_{1}$ in (A.3) respectively, we can derive that
$\Omega^{-1}=\alpha^{-1}H_{1}+H_{2}-\alpha H_{2}BH_{2}+O(\alpha^{2})$. Here,
we have used the fact that $H_{1}BH_{2}=H_{2}BH_{1}=\boldsymbol{0}$.
Consequently, we have
$V(\alpha)=H_{1}(I_{M}\otimes\Sigma_{xx})H_{1}^{\top}+\alpha^{2}H_{2}(I_{M}\otimes\Sigma_{xx})H_{2}^{\top}+O(\alpha^{4})=M^{-1}(\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}_{M}^{\top})\otimes\Sigma_{xx}^{-1}+\alpha^{2}(QQ^{\top})\otimes\Sigma_{xx}+O(\alpha^{4})$.

Let $e_{m}\in\mathbb{R}^{M}$ be the unit vector with $1$ in the $m$-th entry.
Then we should have
$\widehat{\theta}^{(m)}=(e_{m}\otimes\mbox{\boldmath$1$}_{p})^{\top}\widehat{\theta}^{*}$.
In addition,
$(e_{m}\otimes\mbox{\boldmath$1$}_{p})^{\top}\widehat{\mu}(\alpha)=\theta_{0}+O_{p}(\alpha^{2}n^{-3/2})$,
and
$(e_{m}\otimes\mbox{\boldmath$1$}_{p})^{\top}V(\alpha)(e_{m}\otimes\mbox{\boldmath$1$}_{p})=M^{-1}\Sigma_{xx}^{-1}+\alpha^{2}\|a\|^{2}\Sigma_{xx}+O(\alpha^{4})$,
where $\|a\|$ can be calculated by (A.10) and (A.12) in Appendix A.6. Thus,
$\sqrt{n}\big{\\{}(\widehat{\theta}^{(m)}-\theta_{0}+O_{p}(\alpha^{2}n^{-3/2})\big{\\}}\to_{d}N_{p}\Big{(}\boldsymbol{0},\sigma^{2}\big{\\{}M^{-1}\Sigma_{xx}^{-1}+\alpha^{2}\|a\|^{2}\Sigma_{xx}+O(\alpha^{4})\big{\\}}\Big{)}$.
Consequently,

$\sqrt{N}\Big{\\{}(\widehat{\theta}^{(m)}-\theta_{0}+O_{p}(\alpha^{2}n^{-3/2})\Big{\\}}\to_{d}N_{p}\bigg{(}\boldsymbol{0},\sigma^{2}\Big{\\{}\Sigma_{xx}^{-1}+\alpha^{2}M\|a\|^{2}\Sigma_{xx}+O(\alpha^{4})\Big{\\}}\bigg{)}.$

Further note that $\|a\|$ is equal to the $l_{2}$-norm of the first row of
$H_{2}=Q\otimes I_{p}$. We thus complete the proof.

Appendix A.6 The Explicit form of $Q$

One can verify that $A=\widetilde{A}\otimes I_{p}$, where

$\widetilde{A}=\begin{bmatrix}1&0&0&\cdots&0&\gamma&-(1+\gamma)\\\
-(1+\gamma)&1&0&\cdots&0&0&\gamma\\\ \gamma&-(1+\gamma)&1&\cdots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&\gamma&-(1+\gamma)&1\end{bmatrix}\in\mathbb{R}^{M\times M}$

is a circulant matrix. Then we should have
$H_{2}=P_{2}(A_{22}^{P})^{-1}P_{2}^{\top}=(Z\otimes
I_{p})\big{\\{}(Z^{\top}\otimes I_{p})(\widetilde{A}\otimes I_{p})(Z\otimes
I_{p})\big{\\}}^{-1}(Z^{\top}\otimes I_{p})=Q\otimes I_{p}$, where
$Q=Z(Z^{\top}\widetilde{A}Z)^{-1}Z^{\top}$, and $Z$ is defined in (A.2). It
can be easily verified that $Q\widetilde{A}Q=Q$, which implies $Q$ is the
generalized inverse of $\widetilde{A}$. Then it suffices to find the
generalized inverse of $\widetilde{A}$. Since $\mbox{rank}(\widetilde{A})=M-1$
and $\mbox{\boldmath$1$}_{M}$ is the eigenvector of zero eigenvalue, according
to Theorem 2 in Searle, (1979), $Q$ has the following form

$Q=\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}_{M}^{\top}\\}^{-1}-\frac{1}{(1+\gamma)M^{2}}\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}_{M}^{\top}.$
(A.5)

Since the inverse of a circulant matrix is also a circulant matrix, we must
have
$\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}^{-1}$
be a circulant matrix, which has the form of

$\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}^{-1}=\begin{bmatrix}x_{0}&x_{1}&x_{2}&\cdots&x_{M-1}\\\
x_{M-1}&x_{0}&x_{1}&\cdots&x_{M-2}\\\ x_{M-2}&x_{M-1}&x_{0}&\cdots&x_{M-3}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
x_{1}&x_{2}&x_{3}&\cdots&x_{0}\end{bmatrix}.$

Since
$\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}^{-1}=I_{M}$,
we use the first row of
$\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}$
to multiply all columns of
$\\{\widetilde{A}+(1+\gamma)\mbox{\boldmath$1$}_{M}\mbox{\boldmath$1$}^{\top}_{M}\\}^{-1}$.
Then we have the following $M$ equations,

$\begin{cases}(2+\gamma)x_{0}+(1+\gamma)x_{M-1}+\cdots+(1+\gamma)x_{3}+(1+2\gamma)x_{2}+0x_{1}=1\\\
(2+\gamma)x_{1}+(1+\gamma)x_{0}+\cdots+(1+\gamma)x_{4}+(1+2\gamma)x_{3}+0x_{2}=0\\\
\cdots\\\
(2+\gamma)x_{M-1}+(1+\gamma)x_{M-2}+\cdots+(1+\gamma)x_{2}+(1+2\gamma)x_{1}+0x_{0}=0.\end{cases}$

Sum up all the equations above, we have $(1+\gamma)\sum^{M-1}_{m=0}x_{m}=1/M$.
Use all $M$ equations above minus $(1+\gamma)\sum^{M-1}_{m=0}x_{m}=1/M$, we
then have

$\begin{cases}x_{0}-(1+\gamma)x_{1}+\gamma
x_{2}=(x_{0}-x_{1})-\gamma(x_{1}-x_{2})=1-\frac{1}{M}\\\
x_{1}-(1+\gamma)x_{2}+\gamma
x_{3}=(x_{1}-x_{2})-\gamma(x_{2}-x_{3})=-\frac{1}{M}\\\ \cdots\\\
x_{M-1}-(1+\gamma)x_{0}+\gamma
x_{1}=(x_{M-1}-x_{0})-\gamma(x_{0}-x_{1})=-\frac{1}{M}.\end{cases}$ (A.6)

We then study the cases of $\gamma>0$ and $\gamma=0$ subsequently.

Case 1 $(\gamma>0)$. From equations at (A.6), we can verify that

$x_{k}-x_{k+1}=\frac{1}{\gamma^{k}}(x_{0}-x_{1})+\frac{1}{M}S_{k}-\frac{1}{\gamma^{k}}\text{\
for\ }k=1,2,...,M-1,$ (A.7)

where $S_{m}=\sum^{m}_{k=1}\gamma^{-k}$. Sum up all the equations in (A.7), we
have

$0=\sum^{M-1}_{m=0}(x_{m}-x_{m+1})=\gamma
S_{M}(x_{0}-x_{1})+\frac{1}{M}\sum^{M-1}_{m=1}S_{m}-S_{M-1}.$

Therefore, we can obtain that

$\beta=x_{0}-x_{1}=\frac{MS_{M-1}-\sum^{M-1}_{m=1}S_{m}}{M\gamma S_{M}}.$
(A.8)

Again, by (A.7) we can verify that

$x_{m}=x_{0}-\beta\gamma S_{m}-\frac{1}{M}\sum^{m-1}_{k=1}S_{k}+S_{m-1}\text{\
for\ }m=2,...,M-1.$ (A.9)

Thus we can verify that

$\frac{1}{(1+\gamma)M}=\sum^{M-1}_{m=0}x_{m}=Mx_{0}-a\gamma\sum^{M-1}_{m=1}S_{m}-\frac{1}{M}\sum^{M-1}_{m=2}\sum^{m-1}_{k=1}S_{k}+\sum^{M-2}_{m=1}S_{m},$

which implies

$x_{0}=\frac{1}{(1+\gamma)M^{2}}+\frac{\beta\gamma}{M}\sum^{M-1}_{m=1}S_{m}+\frac{1}{M^{2}}\sum^{M-1}_{m=2}\sum^{m-1}_{k=1}S_{k}-\frac{1}{M}\sum^{M-2}_{m=1}S_{m}.$

Since $Q$ is also a circulant matrix, we assume its first row is
$(a_{0},...,a_{M-1})$. According to (A.5) and (A.9), we have

$\begin{cases}a_{0}=&x_{0}-\frac{1}{(1+\gamma)M^{2}}=\frac{\beta\gamma}{M}\sum^{M-1}_{m=1}S_{m}+\frac{1}{M^{2}}\sum^{M-1}_{m=2}\sum^{m-1}_{k=1}S_{k}-\frac{1}{M}\sum^{M-2}_{m=1}S_{m}\\\
a_{1}=&a_{0}-\beta\\\ a_{m}=&a_{0}-\beta\gamma
S_{m}-\frac{1}{M}\sum^{m-1}_{k=1}S_{k}+S_{m-1}\text{\ for\
}m=2,...,M-1,\end{cases}$ (A.10)

where $\beta$ is defined in (A.8). Then we can write
$Q=(a,J^{M-1}a,\dots,Ja)^{\top}$, where
$a=(a_{0},...,a_{M-1})^{\top}\in\mathbb{R}^{M}$, and $J\in\mathbb{R}^{M\times
M}$ is a permutation matrix defined as

$J=\begin{bmatrix}0&1&0&\cdots&0&0\\\ 0&0&1&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\ 0&0&0&\cdots&0&1\\\
1&0&0&\cdots&0&0\end{bmatrix}.$ (A.11)

Case 2 $(\gamma=0)$. From equations at (A.6), we have

$x_{0}-x_{1}=1-\frac{1}{M},\ x_{1}-x_{2}=-\frac{1}{M},\
\dots,x_{M-1}-x_{0}=-\frac{1}{M}.$

By similar procedures used above, we can obtain that
$x_{0}=(M-1)/(2M)+1/M^{2}$ and $x_{m}=x_{0}-(M-m)/M$ for $m=1,\dots,M-1$. Then
by (A.5), we have

$\begin{cases}a_{0}=&x_{0}-\frac{1}{M^{2}}=\frac{M-1}{2M}\\\
a_{m}=&a_{0}-\frac{M-m}{M}=\frac{2m-M-1}{2M}\text{ for
}m=1,\dots,M-1.\end{cases}$ (A.12)

Additionally, by (A.5) and the fact that
$(1+\gamma)\sum_{m=0}^{M-1}x_{m}=1/M$, it can be easily verified that
$\mbox{\boldmath$1$}_{M}^{\top}a=\sum_{m=0}^{M-1}a_{m}=0$ for any $\gamma\geq
0$.

APPENDIX B

######  Lemma 1.

Suppose $\Sigma_{xx}\in\mathbb{R}^{p\times p}$ is a symmetric positive
definite matrix with eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{p}>0$.
Consider a block matrix defined by

$\displaystyle\Omega=\begin{bmatrix}I_{p}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\gamma
I_{p}&-\Delta\\\
-\Delta&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}&\gamma
I_{p}\\\ \gamma
I_{p}&-\Delta&I_{p}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\boldsymbol{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\cdots&\gamma
I_{p}&-\Delta&I_{p}\end{bmatrix}\in\mathbb{R}^{q\times q},$

where $\Delta=(1+\gamma)I_{p}-\alpha\Sigma_{xx}\in\mathbb{R}^{p\times p}$.
Assume $\gamma\neq 1$, $\alpha\neq 0$ and $\alpha\neq 2(1+\gamma)/\lambda_{j}$
for each $1\leq j\leq p$. Then we have $\Omega$ is invertible.

###### Proof.

To prove the invertibility of $\Omega$, we only need to prove that $0$ is not
an eigenvalue of $\Omega$ under the assumed conditions. To this end, we first
note that $\Omega$ is a block circulant matrix. Hence, we can analyze the
eigenvalues of $\Omega$ with help of the results in Kaveh and Rahami, (2011).
Specifically, we denote $\omega_{m}=(e^{2\pi i/M})^{m-1}\in\mathbb{C}$ for
$1\leq m\leq M$, where $i=\sqrt{-1}$ is the imaginary unit. We can see that
$|\omega_{m}|=1$ and $\omega_{m}^{M}=1$ for each $1\leq m\leq M$. Define
$H(\omega)=I_{p}+\omega^{M-2}\gamma I_{p}-\omega^{M-1}\Delta$ as a map from
$\mathbb{C}$ to $\mathbb{C}^{p\times p}$. Then by the results in Kaveh and
Rahami, (2011), the eigenvalues of $\Omega$ will be the union of the
eigenvalues of matrices $H(\omega_{m})$ for $1\leq m\leq M$. We next show that
$0$ cannot be an eigenvalue of $H(\omega_{m})$ for any $1\leq m\leq M$. Since
$\Sigma_{xx}$ has eigenvalues $\lambda_{1},\dots,\lambda_{p}$, the eigenvalues
of $H(\omega_{m})$ should be
$h(\omega_{m})=1+\omega_{m}^{M-2}\gamma-\omega_{m}^{M-1}(1+\gamma-\alpha\lambda_{j})$
for $1\leq j\leq p$. Suppose $h(\omega_{m})=0$ for some $m$. Then we have
$0=\omega_{m}^{2}h(\omega_{m})=\omega_{m}^{2}+(\alpha\lambda_{j}-1-\gamma)\omega_{m}+\gamma$,
i.e., $\omega_{m}$ is the root of the quadratic equation
$\omega^{2}+(\alpha\lambda_{j}-1-\gamma)\omega+\gamma=0$ with real
coefficients. (1) Suppose $\omega_{m}\notin\mathbb{R}$. Then we must have
$|\omega_{m}|^{2}=\omega_{m}\accentset{{\cc@style\underline{\mskip
10.0mu}}}{\omega}_{m}=\gamma\neq 1$. This is a contradiction. (2) Suppose
$\omega_{m}\in\mathbb{R}$ . Then it must be $-1$ or $1$. But this is
impossible under the conditions $\alpha\neq 0$ and $\alpha\neq
2(1+\gamma)/\lambda_{j}$ for each $1\leq j\leq p$. Thus, we have proved the
result. ∎

######  Lemma 2.

Suppose $\Sigma_{xx}\in\mathbb{R}^{p\times p}$ is a symmetric positive
definite matrix with eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{p}>0$.
Consider a block matrix defined $D=[\Delta,-\gamma
I_{p};I_{p},\boldsymbol{0}]\in\mathbb{R}^{(2p)\times(2p)}$, where
$\Delta=(1+\gamma)I_{p}-\alpha\Sigma_{xx}$, $\alpha>0$ and $\gamma\geq 0$.

  1. (a)

Suppose that $0<\alpha<2(1+\gamma)/\lambda_{1}$. Then, $\gamma\leq\rho(D)<1$
if $0\leq\gamma<1$, $\rho(D)>1$ if $\gamma>1$, and $\rho(D)=1$ if $\gamma=1$.

  2. (b)

If $\gamma=0$, then
$\rho(D)=\max\\{|1-\alpha\lambda_{1}|,|1-\alpha\lambda_{p}|\\}$. When
$\alpha=2/(\lambda_{1}+\lambda_{p})$, $\rho(D)$ attains the minimum value as
$\rho(D)=(\lambda_{1}-\lambda_{p})/(\lambda_{1}+\lambda_{p})$.

  3. (c)

If
$\gamma\in[(1-\sqrt{\alpha\lambda_{j}})^{2},(1+\sqrt{\alpha\lambda_{j}})^{2}]$
for all $1\leq j\leq p$, then $\rho(D)=\sqrt{\gamma}$. When
$\alpha=4/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$ and
$\gamma=(\sqrt{\lambda_{1}}-\sqrt{\lambda_{p}})^{2}/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})^{2}$,
then
$\rho(D)=(\sqrt{\lambda_{1}}-\sqrt{\lambda_{p}})/(\sqrt{\lambda_{1}}+\sqrt{\lambda_{p}})$.

  4. (d)

Suppose $0<\alpha<1/\lambda_{1}$. If $\gamma=0$, then
$\rho(D)=1-\alpha\lambda_{p}$; If $0<\gamma<(1-{\alpha\lambda_{p}})^{2}$, then
$\rho<1-\alpha\lambda_{p}$. When $\gamma=(1-\sqrt{\alpha\lambda_{p}})^{2}$,
$\rho(D)$ attains the minimum value as $\rho(D)=1-\sqrt{\alpha\lambda_{p}}$.

  5. (e)

If $\alpha>2(1+\gamma)/\lambda_{1}$ or $\gamma>1$, then $\rho(D)>1$.

###### Proof.

We prove (a)-(e) subsequently.

Proof of (a). We only consider the case of $0\leq\gamma<1$, since the case of
$\gamma\geq 1$ can be derived similarly. Suppose $\xi\in\mathbb{C}$ is an
arbitrary eigenvalue of $D$, and
$v=(v_{1}^{\top},v_{2}^{\top})^{\top}\in\mathbb{C}^{2p}$ is the corresponding
eigenvector, where $v_{1}\in\mathbb{C}^{p}$ and $v_{2}\in\mathbb{C}^{p}$. By
$Dv=\xi v$, we have $\Delta v_{1}-\gamma v_{2}=\xi v_{1}$ and $v_{1}=\xi
v_{2}$. (1) If $\gamma=0$, one can derive that the nonzero eigenvalues of $D$
are the same as those of $\Delta$. Then by assumed condition
$0<\alpha<2(1+\gamma)/\lambda_{1}$, we should have
$\rho(D)=\rho(\Delta)=\max_{1\leq j\leq p}|1-\alpha\lambda_{j}|<1$ for all
$1\leq j\leq p$. Consequently, $\rho(D)\in[0,1)=[\gamma,1)$. (2) If
$0<\gamma<1$, one can verify that $D$ is nonsingular. Thus $\xi\neq 0$. Then
we have $\Delta v_{1}=(\gamma/\xi+\xi)v_{1}$. In other words, $\gamma/\xi+\xi$
is the eigenvalue of the real symmetric matrix $\Delta$. Hence
$\gamma/\xi+\xi$ must be a real number. (2.1) Suppose that
$\xi\notin\mathbb{R}$, i.e., $\xi=a+bi$ with $b\neq 0$. Since
$\gamma/\xi+\xi=(\gamma/|\xi|^{2})\accentset{{\cc@style\underline{\mskip
10.0mu}}}{\xi}+\xi\in\mathbb{R}$, we should have $|\xi|^{2}=\gamma$. Thus
$|\xi|=\sqrt{\gamma}$. (2.2) Suppose $\xi\in\mathbb{R}$. By Weyl’s inequality
and condition $0<\alpha<2(1+\gamma)/\lambda_{1}$ , we know that the
eigenvalues of $\Delta$ are contained in the interval
$\big{(}-(1+\gamma),1+\gamma\big{)}$. Thus we should have
$-(1+\gamma)<(\gamma/\xi+\xi)<1+\gamma$. (2.2.1) If $\xi>0$, by
$(\gamma/\xi+\xi)<1+\gamma$ we have
$\xi^{2}-(1+\gamma)\xi+\gamma=(\xi-1)(\xi-\gamma)<0$, or equivalently,
$\xi\in(\gamma,1)$. (2.2.2) If $\xi<0$, by $-(1+\gamma)<(\gamma/\xi+\xi)$ we
have $\xi^{2}+(1+\gamma)\xi+\gamma=(\xi+1)(\xi+\gamma)<0$, or equivalently,
$\xi\in(-1,-\gamma)$. Thus, if $\xi\in\mathbb{R}$, then $|\xi|\in(\gamma,1)$.
Since $\sqrt{\gamma}\in(\gamma,1)$ if $0<\gamma<1$, we conclude that
$\rho(D)\in(\gamma,1)$. By (1) and (2) above, we have proved that
$\rho(D)\in[\gamma,1)$ if $0<\alpha<2(1+\gamma)/\lambda_{1}$ and
$0\leq\gamma<1$.

Proof of (b). Similar to the proof of (a), $D$ has the same nonzero
eigenvalues as those of $\Delta$ if $\gamma=0$. Since the eigenvalues of
$\Delta$ are $1-\alpha\lambda_{j}$ for $1\leq j\leq p$, claims in (b) follow
immediately.

Proof of (c). The proof can be found in Polyak, (1964). We give a simple
version for the sake of completeness. If $\gamma=0$, then we must have
$\alpha\lambda_{j}=1$ for all $1\leq j\leq p$ under the assumed conditions.
One can verify that $\Delta$ will be a zero matrix in this case. Consequently,
the eigenvalues of $D$ are all zero, i.e., $\rho(D)=0$. We then consider the
case of $\gamma>0$. First, note that when $\gamma>0$, $D$ is nonsingular,
i.e., $D$ has no zero eigenvalue. Then we construct a matrix function of
$\lambda$ as

$G(\lambda)=\begin{bmatrix}1+\gamma-\alpha\lambda&-\gamma\\\
1&0\end{bmatrix}\in\mathbb{R}^{2\times 2}.$

One can verify that the roots of characteristic polynomial $\det(\xi
I_{2p}-D)$ are the same as the roots of the $p$ characteristic polynomials
$\det(\xi I_{2}-G(\lambda_{j}))$ ($1\leq j\leq p$). Hence the eigenvalues of
$D$ are the union of the eigenvalues of $G(\lambda_{j})$ for $1\leq j\leq p$.
Therefore, we turn to investigate the spectrum of $G(\lambda)$. The
characteristic polynomial of $D$ is $\det(\xi
I_{2}-G(\lambda))=\xi^{2}-(1+\gamma-\alpha\lambda)\xi+\gamma$. We find that
the discriminant of the corresponding quadratic equation
$(1+\gamma-\alpha\lambda)^{2}-4\gamma\leq 0$ if and only if
$\gamma\in[(1-\sqrt{\alpha\lambda})^{2},(1+\sqrt{\alpha\lambda})^{2}]$. In
this case, the roots of this equation have the same magnitude $\sqrt{\gamma}$.
The first claim in (b) follows immediately. The second claim can be obtained
by direct calculations for the given $\alpha$ and $\gamma$. Thus we have
proved (c).

Proof of (d). The case of $\gamma=0$ follows immediately from (b). We then
investigate the case of $0<\gamma<(1-\alpha\lambda_{p})^{2}$. (1) First, if
$\gamma\in\big{[}(1-\sqrt{\alpha\lambda_{p}})^{2},(1-\alpha\lambda_{p})^{2}\big{)}$,
then
$\gamma\in\big{[}(1-\sqrt{\alpha\lambda_{j}})^{2},(1+\sqrt{\alpha\lambda_{j}})^{2}\big{]}$
for all $1\leq j\leq p$ under the condition $0<\alpha<1/\lambda_{1}$. Then by
(c), we know that $\rho(D)=\sqrt{\gamma}$. Hence,
$\rho(D)<1-\alpha\lambda_{p}$ if
$\gamma\in\big{[}(1-\sqrt{\alpha\lambda_{p}})^{2},(1-\alpha\lambda_{p})^{2}\big{)}$.
(2) We then consider the case of $0<\gamma<(1-\sqrt{\alpha\lambda_{p}})^{2}$.
In this case, there must be some $1\leq j_{1}\leq p$ such that
$\gamma\geq(1-\sqrt{\alpha\lambda_{j}})^{2}$ for $j<j_{1}$, and
$\gamma<(1-\sqrt{\alpha\lambda_{j}})^{2}$ for $j\geq j_{1}$. (2.1) For
$j<j_{1}$, it follows the proof of (c) that
$\rho(G(\lambda_{j}))=\sqrt{\gamma}<1-\sqrt{\alpha\lambda}<1-\alpha\lambda_{p}$.
(2.2) For $j\geq j_{1}$, the characteristic equation
$\xi^{2}-(1+\gamma-\alpha\lambda_{j})\xi+\gamma=0$ should have two different
real roots $\xi_{1},\xi_{2}$, since the discriminant of this equation is
larger than $0$. Then by Vieta’s formulas, we have
$\xi_{1}+\xi_{2}=1+\gamma-\alpha\lambda_{j}>0$ and $\xi_{1}\xi_{2}=\gamma>0$.
Thus, both $\xi_{1}$ and $\xi_{2}$ should be positive. By (a), we cam conclude
that $\xi_{1},\xi_{2}$ should be contained in $[\gamma,1)$. Since
$\xi_{1}\xi_{2}=\gamma$, we can further verify that $\xi_{1},\xi_{2}$ must be
strictly larger that $\gamma$, i.e., $\xi_{1},\xi_{2}\in(\gamma,1)$.
Consequently,
$\xi_{1}=(1+\gamma-\alpha\lambda_{j})-\xi_{2}<1-\alpha\lambda_{j}\leq
1-\alpha\lambda_{p}$. Obviously, we also have $\xi_{2}<1-\alpha\lambda_{p}$.
Therefore, we should have $\rho(G(\lambda_{j}))<1-\alpha\lambda_{p}$ for all
$j_{1}<j\leq p$. It follows from (2.1) and (2.2) that
$\rho(D)<1-\alpha\lambda_{p}$ if
$\gamma\in\big{(}0,(1-\sqrt{\alpha\lambda_{p}})^{2}\big{)}$. Hence by (1) and
(2), we have proved that $\rho(D)<1-\alpha\lambda_{p}$ if
$0<\gamma<(1-\alpha\lambda_{p})^{2}$.

From (1) above, we know that $\rho(D)=1-\sqrt{\alpha\lambda_{p}}$ if
$\gamma=(1-\sqrt{\alpha\lambda_{p}})^{2}$, and
$\rho(D)>1-\sqrt{\alpha\lambda_{p}}$ if
$\gamma\in\big{(}(1-\sqrt{\alpha\lambda_{p}})^{2},(1-\alpha\lambda_{p})^{2}\big{)}$.
Hence, to prove that the minimal $\rho(D)$ is achieved by choosing
$\gamma=(1-\sqrt{\alpha\lambda_{p}})^{2}$, we only need to verify that if
$0<\gamma<(1-\sqrt{\alpha\lambda_{p}})^{2}$, then
$\rho(D)>1-\sqrt{\alpha\lambda_{p}}$. In fact, if
$0<\gamma<(1-\sqrt{\alpha\lambda_{p}})^{2}$, one can verify that the equation
$h(\xi)=\xi^{2}-(1+\gamma-\alpha\lambda_{p})\xi+\gamma=0$ has two different
positive real roots $\xi_{1},\xi_{2}$ by similar arguments in (2.2) above. To
prove that one of $\xi_{1}$ and $\xi_{2}$ is larger that
$1-\sqrt{\alpha\lambda_{p}}$, it suffices to show
$h(1-\sqrt{\alpha\lambda_{p}})<0$. This is true, since it can be easily
verified that
$h(1-\sqrt{\alpha\lambda_{p}})=-(1-\sqrt{\alpha\lambda_{p}})^{2}\sqrt{\alpha\lambda_{p}}+\gamma\sqrt{\alpha\lambda_{p}}<0$.
Thus, we complete the proof of (d).

Proof of (e). By arguments used above, one can easily show that, if
$\gamma>1$, then $\rho(D)>1$. Hence, we only prove the result under
$\alpha>2(1+\gamma)/\lambda_{1}$ and $0\leq\gamma\leq 1$. The case of
$\gamma=0$ is also trivial, hence we only deal with the case of $0<\gamma\leq
1$. Consider the characteristic equation $\det(\xi
I_{2}-G(\lambda_{1}))=\xi^{2}+b\xi+\gamma=0$, where
$b=-(1+\gamma-\alpha\lambda_{1})$. If $\alpha>2(1+\gamma)/\lambda_{1}$, we
have $b>1+\gamma$. Hence, the discriminant of this equation
$b^{2}-4\gamma>(1-\gamma)^{2}\geq 0$. Since $0<\gamma\leq 1$, one root of this
equation should be
$(-b-\sqrt{b^{2}-4\gamma})/2<\\{-(1+\gamma)-(1-\gamma)\\}/2=-1$, i.e.,
$\rho(G(\lambda_{1}))>1$. Thus, in either case, we should have $\rho(D)>1$. We
complete the proof. ∎

REFERENCES

## References

  * Abadi et al., (2016) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., et al. (2016).  Tensorflow: Large-scale machine learning on heterogeneous distributed systems.  arXiv preprint arXiv:1603.04467.
  * Battey et al., (2018) Battey, H., Fan, J., Liu, H., Lu, J., and Zhu, Z. (2018).  Distributed testing and estimation under sparse high dimensional models.  The Annals of Statistics, 46(3):1352.
  * Bickel et al., (2008) Bickel, P. J., Levina, E., et al. (2008).  Regularized estimation of large covariance matrices.  The Annals of Statistics, 36(1):199–227.
  * Chang et al., (2017) Chang, X., Lin, S.-B., and Zhou, D.-X. (2017).  Distributed semi-supervised learning with kernel ridge regression.  The Journal of Machine Learning Research, 18(1):1493–1514.
  * (5) Chen, X., Lee, J. D., Li, H., and Yang, Y. (2021a).  Distributed estimation for principal component analysis: An enlarged eigenspace analysis.  Journal of the American Statistical Association, pages 1–12.
  * (6) Chen, X., Lee, J. D., Tong, X. T., Zhang, Y., et al. (2020a).  Statistical inference for model parameters in stochastic gradient descent.  The Annals of Statistics, 48(1):251–273.
  * (7) Chen, X., Liu, W., Mao, X., and Yang, Z. (2020b).  Distributed high-dimensional regression under a quantile loss function.  The Journal of Machine Learning Research, 21(182):1–43.
  * Chen et al., (2019) Chen, X., Liu, W., and Zhang, Y. (2019).  Quantile regression under memory constraint.  The Annals of Statistics, 47(6):3244–3273.
  * (9) Chen, X., Liu, W., and Zhang, Y. (2021b).  First-order newton-type estimator for distributed estimation and inference.  Journal of the American Statistical Association, pages 1–17.
  * Chen and Xie, (2014) Chen, X. and Xie, M.-g. (2014).  A split-and-conquer approach for analysis of extraordinarily large data.  Statistica Sinica, pages 1655–1684.
  * Deng et al., (2009) Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K., and Fei-Fei, L. (2009).  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee.
  * Fan et al., (2019) Fan, J., Wang, D., Wang, K., and Zhu, Z. (2019).  Distributed estimation of principal eigenspaces.  The Annals of Statistics, 47(6):3009.
  * He et al., (2016) He, K., Zhang, X., Ren, S., and Sun, J. (2016).  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778.
  * Horn and Johnson, (2012) Horn, R. A. and Johnson, C. R. (2012).  Matrix analysis.  Cambridge university press.
  * Kaveh and Rahami, (2011) Kaveh, A. and Rahami, H. (2011).  Block circulant matrices and applications in free vibration analysis of cyclically repetitive structures.  Acta Mechanica, 217(1):51–62.
  * Lin et al., (2020) Lin, S.-B., Wang, D., and Zhou, D.-X. (2020).  Distributed kernel ridge regression with communications.  J. Mach. Learn. Res., 21:93–1.
  * Luo and Song, (2020) Luo, L. and Song, P. X.-K. (2020).  Renewable estimation and incremental inference in generalized linear models with streaming data sets.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):69–97.
  * Nesterov et al., (2018) Nesterov, Y. et al. (2018).  Lectures on convex optimization, volume 137.  Springer.
  * Paszke et al., (2017) Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A. (2017).  Automatic differentiation in pytorch.  In NIPS 2017 Autodiff Workshop: The Future of Gradient-based Machine Learning Software and Techniques, Long Beach, CA, US, December 9, 2017,.
  * Polyak, (1964) Polyak, B. T. (1964).  Some methods of speeding up the convergence of iteration methods.  Ussr computational mathematics and mathematical physics, 4(5):1–17.
  * Polyak, (1987) Polyak, B. T. (1987).  Introduction to optimization. optimization software.  Inc., Publications Division, New York, 1.
  * Polyak and Juditsky, (1992) Polyak, B. T. and Juditsky, A. B. (1992).  Acceleration of stochastic approximation by averaging.  SIAM journal on control and optimization, 30(4):838–855.
  * Rao et al., (1973) Rao, C. R., Rao, C. R., Statistiker, M., Rao, C. R., and Rao, C. R. (1973).  Linear statistical inference and its applications, volume 2.  Wiley New York.
  * Robbins and Monro, (1951) Robbins, H. and Monro, S. (1951).  A stochastic approximation method.  The annals of mathematical statistics, pages 400–407.
  * Searle, (1979) Searle, S. (1979).  On inverting circulant matrices.  Linear Algebra and its Applications, 25:77–89.
  * Seber, (2008) Seber, G. A. (2008).  A matrix handbook for statisticians, volume 15.  John Wiley & Sons.
  * Sutskever et al., (2013) Sutskever, I., Martens, J., Dahl, G., and Hinton, G. (2013).  On the importance of initialization and momentum in deep learning.  In International conference on machine learning, pages 1139–1147. PMLR.
  * Tang et al., (2020) Tang, L., Zhou, L., and Song, P. X.-K. (2020).  Distributed simultaneous inference in generalized linear models via confidence distribution.  Journal of Multivariate Analysis, 176:104567.
  * Toulis et al., (2017) Toulis, P., Airoldi, E. M., et al. (2017).  Asymptotic and finite-sample properties of estimators based on stochastic gradients.  The Annals of Statistics, 45(4):1694–1727.
  * Vershynin, (2018) Vershynin, R. (2018).  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press.
  * Wang et al., (2013) Wang, C., Chen, X., Smola, A., and Xing, E. P. (2013).  Variance reduction for stochastic gradient optimization.  In Proceedings of the 26th International Conference on Neural Information Processing Systems-Volume 1, pages 181–189.
  * Wang et al., (2017) Wang, M., Fang, E. X., and Liu, H. (2017).  Stochastic compositional gradient descent: algorithms for minimizing compositions of expected-value functions.  Mathematical Programming, 161(1-2):419–449.
  * Wang et al., (2019) Wang, X., Yang, Z., Chen, X., and Liu, W. (2019).  Distributed inference for linear support vector machine.  The Journal of Machine Learning Research, 20.
  * Yang et al., (2019) Yang, S., Wang, M., and Fang, E. X. (2019).  Multilevel stochastic gradient methods for nested composition optimization.  SIAM Journal on Optimization, 29(1):616–659.
  * Zhu et al., (2021) Zhu, W., Chen, X., and Wu, W. B. (2021).  Online covariance matrix estimation in stochastic gradient descent.  Journal of the American Statistical Association, (just-accepted):1–30.