# Nesterov’s Acceleration For Approximate Newton

Haishan Ye
<EMAIL_ADDRESS>
Department of Computer Science and Engineering
Shanghai Jiao Tong University
800 Dong Chuan Road, Shanghai, China 200240 Zhihua Zhang
<EMAIL_ADDRESS>
School of Mathematical Sciences
Peking University
Beijing, China 100871

###### Abstract

Optimization plays a key role in machine learning. Recently, stochastic
second-order methods have attracted much attention due to their low
computational cost in each iteration. However, these algorithms might perform
poorly especially if it is hard to approximate the Hessian well and
efficiently. As far as we know, there is no effective way to handle this
problem. In this paper, we resort to Nesterov’s acceleration technique to
improve the convergence performance of a class of second-order methods called
approximate Newton. We give a theoretical analysis that Nesterov’s
acceleration technique can improve the convergence performance for approximate
Newton just like for first-order methods. We accordingly propose an
accelerated regularized sub-sampled Newton. Our accelerated algorithm performs
much better than the original regularized sub-sampled Newton in experiments,
which validates our theory empirically. Besides, the accelerated regularized
sub-sampled Newton has good performance comparable to or even better than
classical algorithms.

##  1 Introduction

Optimization has become an increasingly popular issue in machine learning.
Many machine learning models can be reformulated as the following optimization
problems:

$\min_{x\in{\mathbb{R}}^{d}}F(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x).$ (1)

where each $f_{i}$ is the loss with respect to (w.r.t.) the $i$-th training
sample. There are many examples such as logistic regressions, smoothed support
vector machines, neural networks, and graphical models.

In the era of big data, large-scale optimization algorithms have become an
important challenge. The stochastic gradient descent algorithm (SGD) has been
widely employed to reduce the computational cost [5, 14, 20]. However, SGD has
poor convergence property. Hence, many variants have been proposed to improve
the convergence rate of SGD [11, 22, 23, 27].

For the first-order methods which only make use of the gradient information,
Nesterov’s acceleration technique is a very useful tool [17]. It greatly
improves the convergence of gradient descent [17], proximal gradient descent
[2, 18], and stochastic gradient with variance reduction [1, 13], etc.

Recently, second-order methods have also received great attention due to their
high convergence rate. However, conventional second-order methods are very
costly because they take heavy computational cost to obtain the Hessian
matrices. To conquer this weakness, one proposed a sub-sampled Newton which
only samples a subset of functions $f_{i}$ randomly to construct a sub-sampled
Hessian [21, 3, 25] . Pilanci and Wainwright [19] applied the sketching
technique to alleviate the computational burden of computing Hessian and
brought up _sketch Newton_. Regularized sub-sampled Newton methods were also
devised to deal with the ill-condition problem [7, 21].

In the latest work, Ye et al. [26] cast these stochastic second-order
procedures into a so-called _approximate Newton_ framework. They showed that
if approximate Hessian $H^{(t)}$ satisfies

$\displaystyle(1-\pi)\nabla^{2}F(x^{(t)})\preceq
H^{(t)}\preceq(1+\pi)\nabla^{2}F(x^{(t)}),$ (2)

where $0<\pi<1$, then approximate Newton converges with rate $\pi$. If
$H^{(t)}$ is a poor approximation like $\pi=1-1/\kappa$, where $\kappa$ is the
condition number of object function $F(x)$, then approximate Newton has the
same convergence rate with gradient descent.

Since approximate Newton converges with a linear rate, it is natural to ask
whether approximate Newton can be accelerated just like gradient descent. If
it can be accelerated, can the convergence rate be promoted to
$1-\sqrt{1-\pi}$ compared to original $\pi$? In this paper, we aim to
introduce Nesterov’s acceleration technique to improve the performance of
second-order methods, specifically approximate Newton.

We summarize our work and contribution as follows:

  * •

First, we introduce Nesterov’s acceleration technique to improve the
convergence rate of the stochastic second-order methods (approximate Newton).
This acceleration is very important especially when $n$ and $d$ are close to
each other and object function in question is ill-conditioned. In these cases,
it is very hard to construct a good approximate Hessian with low cost.

  * •

Our theoretical analysis shows that by Nesterov’s acceleration, the
convergence rate of approximate Newton can be improved to $1-\sqrt{1-\pi}$
from original rate $\pi$ where $0<\pi<1$ when the object function is
quadratic. For general smooth convex functions, we also show that the similar
acceleration also holds when the initial point is close to the optimal point.

  * •

We propose Accelerated Regularized Sub-sampled Newton. Compared with classical
stochastic first-order methods, our algorithm shows competitive or even better
performance. This demonstrates the efficiency of the accelerated second-order
method. Our experimental study shows that Nesterov’s acceleration technique
can improve approximate Newton methods effectively. Our experiments also
reveal a fact that adding curvature information properly can always improve
the algorithm’s convergence performance.

##  2 Notation and Preliminaries

We first introduce notation that will be used in this paper. Then, we give
some properties of object function that will be used.

###  2.1 Notation

Given a matrix $A=[a_{ij}]\in{\mathbb{R}}^{m\times n}$ of rank $\ell$ and a
positive integer $k\leq\ell$, its SVD is given as $A=U\Sigma
V^{T}=U_{k}\Sigma_{k}V_{k}^{T}+U_{\setminus
k}\Sigma_{{\setminus}k}V_{{\setminus}k}^{T}$, where $U_{k}$ and
$U_{{\setminus}k}$ contain the left singular vectors of $A$, $V_{k}$ and
$V_{{\setminus}k}$ contain the right singular vectors of $A$, and
$\Sigma=\mathrm{diag}(\sigma_{1},\ldots,\sigma_{\ell})$ with
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{\ell}>0$ are the nonzero
singular values of $A$. Additionally, $\|A\|\triangleq\sigma_{1}$ is the
spectral norm. If $A$ is positive semidefinite, then $U=V$ and the eigenvalue
decomposition of $A$ is the same to singular value decomposition. It also
holds that $\lambda_{i}(A)=\sigma_{i}(A)$, where $\lambda_{i}(A)$ is the
$i$-th largest eigenvalue of $A$. Let $\lambda_{\max}(A)$ and
$\lambda_{\min}(A)$ denote the largest and smallest eigenvalue of $A$,
respectively.

###  2.2 Assumptions

In this paper, we focus on the problem described in Eqn. (1). Moreover, we
will make the following two assumptions.

#### Assumption 1

The objective function $F$ is $\mu$-strongly convex, that is,

$F(y)\geq F(x)+[\nabla F(x)]^{T}(y-x)+\frac{\mu}{2}\|y-x\|^{2},\mbox{ for
}\;\mu>0.$

#### Assumption 2

$\nabla F(x)$ is $L$-Lipschitz continuous, that is,

$\|\nabla F(x)-\nabla F(y)\|\leq L\|y-x\|,\mbox{ for }\;L>0.$

By Assumptions 1 and 2, we define the condition number of function $F(x)$ as:
$\kappa\triangleq\frac{L}{\mu}$.

Besides, we will also use the nation of Lipschitz continuity of
$\nabla^{2}F(x)$ in this paper. We say $\nabla^{2}F(x)$ is $\hat{L}$-Lipschitz
continuous if

$\|\nabla^{2}F(x)-\nabla^{2}F(y)\|\leq\hat{L}\|y-x\|,\mbox{ for }\;\hat{L}>0.$

###  2.3 Row Norm Squares Sampling

The row norm squares sampling matrix $S=D\Omega\in{\mathbb{R}}^{s\times n}$
w.r.t. $A\in{\mathbb{R}}^{n\times d}$ is determined by sampling probability
$p_{i}$, a sampling matrix $\Omega\in{\mathbb{R}}^{s\times n}$ and a diagonal
rescaling matrix $D\in{\mathbb{R}}^{s\times s}$. The sampling probability is
defined as

$p_{i}=\frac{\|A_{i,:}\|^{2}}{\|A\|_{F}^{2}},$

where $A_{i,:}$ means the $i$-th row of $A$. We construct $S$ as follows. For
every $j=1,\dots,s$, independently and with replacement, pick an index $i$
from the set $\\{1,2\dots,m\\}$ with probability $p_{i}$ and set
$\Omega_{ji}=1$ and $\Omega_{jk}=0$ for $k\neq i$ as well as
$D_{jj}=1/\sqrt{p_{i}s}$.

Row norm squares sampling matrix has the following important property.

###### Theorem 1

Let $S\in{\mathbb{R}}^{s\times n}$ be a row norm squares sampling matrix
w.r.t. $A\in{\mathbb{R}}^{n\times d}$, then it holds that

$\mathbb{E}\|A^{T}S^{T}SA-A^{T}A\|\leq\left(\sqrt{\frac{4\cdot\mathrm{sr}(A)\cdot\log(2d)}{s}}+\frac{2\cdot\mathrm{sr}(A)\cdot\log(2d)}{3s}\right)\cdot\|A\|^{2}$

###  2.4 Randomized sketching matrices

We first give the definition of the $\epsilon$-subspace embedding property.
Then we list some useful types of randomized sketching matrices including
Gaussian sketching [10, 12], random sampling [6], count sketch [4, 16, 15].

###### Definition 2

$S\in{\mathbb{R}}^{s\times m}$ is said to be an $\varepsilon$-subspace
embedding matrix for any fixed matrix $A\in{\mathbb{R}}^{m\times d}$, if
$\|SAx\|^{2}=(1\pm\varepsilon)\|Ax\|^{2}$
$(i.e.,(1-\varepsilon)\|Ax\|^{2}\leq\|SAx\|^{2}\leq(1+\varepsilon)\|Ax\|^{2})$
for all $x\in{\mathbb{R}}^{d}$.

#### Gaussian sketching matrix:

The most classical sketching matrix is the Gaussian sketching matrix
$S\in{\mathbb{R}}^{s\times m}$ with i.i.d normal random entries with mean 0
and variance $1/s$. Because of well-known concentration properties of Gaussian
random matrices [24], they are very attractive. Besides,
$s=O(d/\varepsilon^{2})$ is enough to guarantee the $\varepsilon$-subspace
embedding property for any fixed matrix $A\in{\mathbb{R}}^{m\times d}$. And
$s=O(d/\varepsilon^{2})$ is the tightest bound in known types of sketching
matrices. However, Gaussian random matrices are dense, so it is costly to
compute $SA$.

#### Count sketch matrix:

Count sketch matrix $S\in{\mathbb{R}}^{s\times m}$ is of the form that there
is only one non-zero entry uniformly sampling from $\\{1,-1\\}$ in each column
[4]. Hence it is very efficient to compute $SA$, especially when $A$ is a
sparse matrix. To achieve an $\epsilon$-subspace embedding property for
$A\in{\mathbb{R}}^{m\times d}$, $s=O(d^{2}/\varepsilon^{2})$ is sufficient
[15, 24].

Other types of sketching matrices like Sub-sampled Randomized Hadamard
Transformation and detailed properties of sketching matrices can be found in
the survey [24].

##  3 Accelerated Approximate Newton

In practice, it is common that the problem is ill-condition and the data size
$n$ and data dimension $d$ are close to each other. Therefore, conventional
Sketch Newton, sub-sampled Newton and regularized sub-sampled Newton can not
construct a good approximate Hessian efficiently. However, Ye et al. [26]
showed that a poor approximate Hessian will lead to a slow converge rate. In
the extreme case, these approximate Newton methods have the same convergence
rate with gradient descent. To improve the convergence property of approximate
Newton when the Hessian can only be approximate poorly, we resort to
Nesterov’s acceleration technique and propose _accelerated approxiamte
Newton_. We summarize the algorithmic procedure of accelerated approxiamte
Newton as follows.

First, we construct an approximate Hessian $[H^{(t)}$ satisfies

$\displaystyle(1-\pi)[\nabla^{2}F(y^{(t)})]^{-1}\preceq\mathbb{E}\left([H^{(t)}]^{-1}\right)\preceq[\nabla^{2}F(y^{(t)})]^{-1},$
(3)

where $0<\pi<1$. This condition is a litte stronger than the one of
approximate Newton which satisfies Eqn. (2). We will see that Condition (3)
can be easily satisfies in practice in next section.

And we update sequence $x^{(t)}$ as follows,

$\left\\{\begin{aligned} &y^{(t+1)}=(1+\theta)x^{(t)}-\theta x^{(t-1)},\\\
&x^{(t+1)}=y^{(t+1)}-[H^{(t+1)}]^{-1}\nabla F(y^{(t+1)}),\end{aligned}\right.$
(4)

where $\theta$ is chosen in terms of the value of $\pi$. We can see that the
iteration (4) is much like the update procedure of Nesterov’s accelerated
gradient descent but replacing step size with $H^{-1}$. Besides, when
$\theta=0$, the above update procedure reduces to the update step of
approximate Newton [26]. Thus, we refer a class of methods satisfying Eqn. (3)
and algorithm procedure (4) as _accelerated approximate Newton_.

Simlar to approximate Newton method, we can update $x^{(t+1)}$ with a
direction vector $p^{(t)}$ which is an inexact solution of following problem

$\displaystyle\frac{1}{2}p^{T}H^{(t)}p-p^{T}\nabla F(y^{(t)}).$ (5)

There are different ways to solve problem (5) like conjugate gradient.

###  3.1 Theoretical Analysis

Because a general smooth convex object function can be approximated by
quadratic functions in a region close to optimal point, we show the
convergence properties of Algorithm 1 when applied to a quadratic object
function. The following lemma (Lemma 3) gives the detailed reason why we can
only analyze convergence properties of Algorithm 1 applied to a quadratic
object function.

Furthermore, we will first give the convergence analysis of accelerated
approximate Newton which satisfies Eqn. (3) and of algorithm procedure (4) in
Theorem 4.

###### Lemma 3

Let Assumptions 1 and 2 hold. Suppose that $\nabla^{2}F(x)$ exists and is
continuous in a neighborhood of a minimizer $x^{*}$. Let $H^{(t)}$ be an
approximation of $\nabla^{2}F(y^{(t)})$. Consider the iteration (4). If
$x^{(t)}$ is sufficient close to $x^{*}$ then we have the following result

$\displaystyle[\nabla^{2}F(x^{*})]^{-\frac{1}{2}}\nabla F(x^{(t+1)})$
$\displaystyle=$
$\displaystyle\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)$
$\displaystyle\times\left((1+\theta)\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t)})-\theta\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t-1)})\right)$ $\displaystyle+o(\nabla F(x^{(t)}))+o(\nabla
F(x^{(t-1)})).$

Furthermore, if $\nabla^{2}F(x)$ is $\hat{L}$-Lipschitz continuous, then the
above result holds whenever $x^{(t)}$ satisfies

$\|x^{(t)}-x^{*}\|\leq\frac{o(1)}{\hat{L}}.$

Lemma 3 shows that the convergence property of iteration (4) is mainly
determined by
$I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}$
and $\theta$. Lemma 3 also describes such a fact that if $x^{(t-1)}$ and
$x^{(t)}$ are sufficient close to $x^{*}$, that is, $o(F(x^{(t)}))$ and
$o(F(x^{(t-1)}))$ are very small, then the convex function can be well
approximated by a quadratic function. Therefore, we will demonstrate the
convergence analysis of accelerated approximate Newton on the strongly convex
quadratic functions. Because the Hessian of a quadratic function $F(x)$ is a
constant matrix, that is the Lipschitz constant of $\nabla^{2}F(x)$ is zero.
Hence, the result of Lemma 3 degenerates to

$\displaystyle[\nabla^{2}F(x^{*})]^{-\frac{1}{2}}\nabla F(x^{(t+1)})$
$\displaystyle=$
$\displaystyle\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)$
$\displaystyle\times\left((1+\theta)\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t)})-\theta\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t-1)})\right).$

This equation describes a linear dynamic system which contains the convergence
property of ieration (4).

###### Theorem 4

Let $F(x)$ be a quadratic function with Assuption 1 and 2 holding. Let
$H^{(t)}$ be an approximation of $\nabla^{2}F(y^{(t)})$ satisfying Eqn. (3)
and $\mathbb{E}\left([H^{(t)}]^{-1}\right)$ is a constant matrix for all $t$.
We set $\theta=\frac{1-\sqrt{1-\pi}}{1+\sqrt{1-\pi}}-\epsilon_{0}$ with
$0<\epsilon_{0}<1$. $T_{1}$ is a matrix of the form

$\displaystyle T_{i}=\begin{bmatrix}(1+\theta)\pi&-\theta\pi\\\
1&0\end{bmatrix}.$

Let $S_{1}$ and $q$ be the eigenmatrix and larger eigenvalue of $T_{1}$
respectively. Then $q$ is of the value
$q=1-\sqrt{1-\pi}+O(\sqrt{\epsilon_{0}})$. And $c$ is a constant defined by
$c_{1}=\|S_{1}\|\cdot\|S_{1}^{-1}\|$. Let direction vector $p^{(t)}$ satisfy

$\displaystyle\left\lVert H^{(t)}p^{(t)}-\nabla
F(y^{(t)})\right\rVert\leq\frac{\epsilon_{1}}{12c_{1}\sqrt{\kappa}}\left\lVert\nabla
F(y^{(t)})\right\rVert,$ (6)

where $\kappa$ is the condition number of the Hessian matrix
$\nabla^{2}F(x^{\star})$. Then after $t$ iterations, Algorithm 1 has

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(t+1)})\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]\right\rVert\leq
2c_{1}p^{t}\left\lVert\left[\begin{array}[]{c}M_{*}\nabla F(x^{(1)})\\\
M_{*}\nabla F(x^{(0)})\\\ \end{array}\right]\right\rVert,$

with $p=q+\epsilon_{1}$ and $M_{\star}=[\nabla^{2}F(x^{\star})]^{-1/2}$.

From Theorem 4, we can see that if we choose
$\theta=\frac{1-\sqrt{1-\pi}}{1+\sqrt{1-\pi}}-\epsilon_{0}$ and
$\epsilon_{1}=0$, then the convergence rate of Algorithm 1 is
$1-\sqrt{1-\pi}+O(\sqrt{\epsilon_{0}})$ in contrast to $\pi$ of the
conventional regularized sub-sampled Newton method. If we choose $H=(1/L)I$
and $\theta=\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}}$, where
$L=\lambda_{\max}(\nabla^{2}F(x^{\star}))$ and
$\mu=\lambda_{\min}(\nabla^{2}F(x^{\star}))$, then Theorem 4 shows that the
convergence rate of Nesterov’s acceleration for the least square regression
problem is $1-\sqrt{\mu/L}+O(\sqrt{\epsilon_{0}})$.

Algorithm 1 Accelerated Regularized Sub-sample Newton (ARSSN).

1: Input: $x^{(0)}$ and $x^{(1)}$ are initial points sufficient close to
$x^{*}$. $\theta$ is the acceleration parameter.

2: for $t=1,\dots$ until termination do

3: Select a sample set ${\mathcal{S}}$ of size $|{\mathcal{S}}|$ by random
sampling and construct $H^{(t)}$ of the form (7) satisfying Eqn. (3);

4: $y^{(t+1)}=(1+\theta)x^{(t)}-\theta x^{(t-1)}$;

5: Obtain the direction vector $p^{(t)}$ by solving problem .

6: $x^{(t+1)}=y^{(t+1)}-p^{(t+1)}$.

7: end for

##  4 Accelerated Regularized Sub-sampled Newton

In practice, it is common that data size $n$ and data dimension $d$ are close
to each other. Conventional Sketch Newton method is not suitable because
sketching size $|{\mathcal{S}}|$ will be less than $d$, and the sketched
Hessian is not invertible. Hence, adding a proper regularizer is a good
approach. On the other hand, sub-sampled Newton method needs lots of samples
when the problem is ill-conditioned. Hence, regularized sub-sample Newton
methods are proppsed [21]. However, Ye et al. [26] showed that regularized
sub-sample Newton will converge slowly as the sample size decreases. In the
extreme case, regularized sub-sample Newton has the same convergence rate with
gradient descent.

To conquer the weakness of slow convergence rate of regularized sub-sample
Newton, we resort to Nesterov’s acceleration technique and propose accelerated
regularized sub-sample Newton method in Algorithm 1. Just like conventional
second-order methods, we assume initial point $x^{(0)}$ and $x^{(1)}$ are
sufficient close to optimal point $x^{\star}$ in our algorithms.

In Algorithm 1, $H^{(t)}$ is an approximation of $\nabla^{2}F(y^{(t)})$
constructed by random sampling. $H^{(t)}$ is of the following structure

$\displaystyle H^{(t)}=\hat{H}^{(t)}+\alpha^{(t)}I.$ (7)

$\hat{H}^{(t)}$ is the sub-sampled Hessian. And $\alpha^{(t)}I$ is the
regularizer. $H^{(t)}$ also satisfies Algorithm 1 specifies the way to
construct the approximate Hessian by random sampling and is a kind of
accelerated approximate Newton. Therefore, the convergence property of
Algorithm 1 can be analyzed by Theorem 4.

First, we consider the following case that the Hessian of $\nabla^{2}F(x)$ has
the following structure

$\displaystyle\nabla^{2}F(x)=B(x)^{T}B(x),\>B(x)\in{\mathbb{R}}^{n\times d}.$
(8)

For this case, we construct approximate Hessian’s as

$\displaystyle H^{(t)}=[S^{(t)}B^{(t)}]^{T}S^{(t)}B^{(t)}+\alpha^{(t)}I.$ (9)

$S^{(t)}$ is a random sampling matrix or other sketching matrix. And
$\alpha^{(t)}$ is a regularizer scaler.

###### Theorem 5

Let $F(x)$ be a quadratic function with Assumption 1 and 2 holding.
$\nabla^{2}F(x)$ satisfies Eqn. (8). Given sample size parameter $0<c<1$,
$S^{(t)}\in{\mathbb{R}}^{s\times n}$ is a row norm squares sampling matrix
w.r.t. $B^{(t)}$ with $s=O(c^{-2}\cdot\mathrm{sr}(B^{(t)})\log d)$. And set
regularizer $\alpha=c\|B^{(t)}\|^{2}$. Construct the approximate Hessian
$H^{(t)}=[S^{(t)}B^{(t)}]^{T}S^{(t)}B^{(t)}+\alpha^{(t)}I$. Let
$\pi=\frac{2c\kappa}{1+2c\kappa}$, and set parameters as Theorem 4, Algorithm
1 converges as

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(t+1)})\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]\right\rVert\leq
2c_{1}\left(1-\frac{1}{\sqrt{1+2c\kappa}}+O(\sqrt{\epsilon_{0}})+\epsilon_{1}\right)^{t}\left\lVert\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(1)})\\\ M_{*}\nabla F(x^{(0)})\\\ \end{array}\right]\right\rVert,$

with $M_{\star}=[\nabla^{2}F(x^{\star})]^{-1/2}$ and $\epsilon_{0}$,
$\epsilon_{1}$ defined in Theorem 4.

Proof  For notation convenience, we will omit superscript and just use $S$,
$B$ and $\alpha$ instead of $S^{(t)}$, $B^{(t)}$ and $\alpha^{(t)}$. By
Theorem 1 and omit the high order term, we have

$\mathbb{E}\|B^{T}B-B^{T}S^{T}SB\|\leq c\|B\|^{2}.$

Hence, in expectation, we have

$\displaystyle|y^{T}(B^{T}B-B^{T}S^{T}SB)y|\leq c\|B\|^{2}\cdot\|y\|^{2}$
$\displaystyle\Rightarrow$ $\displaystyle
y^{T}B^{T}S^{T}SBy-c\|B\|^{2}\cdot\|y\|^{2}\leq y^{T}B^{T}By\leq
y^{T}B^{T}S^{T}SBy+c\|B\|^{2}\cdot\|y\|^{2}$ $\displaystyle\Rightarrow$
$\displaystyle y^{T}B^{T}By\leq y^{T}H^{(t)}y\leq
y^{T}B^{T}By+2c\|B\|^{2}\cdot\|y\|^{2}$ $\displaystyle\Rightarrow$
$\displaystyle y^{T}[B^{T}B]^{-1}y\geq y^{T}[H^{(t)}]^{-1}y\geq
y^{T}\left(B^{T}B+2c\|B\|^{2}\cdot I\right)^{-1}y$

Furthermore, we have

$\displaystyle y^{T}B^{T}By+2c\|B\|^{2}\cdot\|y\|^{2}\leq
y^{T}B^{T}By+\frac{2c\|B\|^{2}}{\sigma_{\min}(B)^{2}}y^{T}B^{T}By$
$\displaystyle\Rightarrow$ $\displaystyle B^{T}B+2c\|B\|^{2}\cdot
I\preceq\left(1+2c\kappa\right)B^{T}B$ $\displaystyle\Rightarrow$
$\displaystyle\left(1-\frac{2c\kappa}{1+2c\kappa}\right)B^{T}B\preceq\left(B^{T}B+2c\|B\|^{2}\cdot
I\right)^{-1}$

Thus, we have

$\left(1-\frac{2c\kappa}{1+2c\kappa}\right)\cdot\left[\nabla^{2}F(x^{\star})\right]^{-1}\preceq\mathbb{E}\left([H^{(t)}]^{-1}\right)\preceq\left[\nabla^{2}F(x^{\star})\right]^{-1}.$

Then final convergence property can be obtained by Theorem 4 using
$\pi=\frac{2c\kappa}{1+2c\kappa}$.

Then, we consider the case that each $f_{i}(x)$ and $F(x)$ in (1) have the
following properties:

$\displaystyle\max_{1\leq i\leq n}\|\nabla^{2}f_{i}(x)\|\leq K<\infty,$ (10)
$\displaystyle\lambda_{\min}(\nabla^{2}F(x))\geq\sigma>0.$ (11)

In this case, we do not need the Hessian of the speciall form (9) but need
that each individual Hessian is upper bounded. We construct the approximate
Hessian $H^{(t)}$ just by uniformly sampling as follows

$\displaystyle
H^{(t)}=\frac{1}{|{\mathcal{S}}|}\sum_{j\in\mathcal{S}}\nabla^{2}f_{j}(x^{(t)})+\alpha^{(t)}I$
(12)

###### Theorem 6

Let $F(x)$ be a quadratic function with Assumption 1 and 2 holding. And Eqns
(10) and (11) also hold. We construct the approximate Hessian $H^{(t)}$ as
Eqn. (12) by uniformly sampling with sample size
$|{\mathcal{S}}|=O(c^{-2}K^{2}\log d)$. And set regularizer
$\alpha^{(t)}=c\|\nabla^{2}F(x^{(t)})\|^{2}$. Set parameters as Theorem 4,
Algorithm 1 converges as

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(t+1)})\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]\right\rVert\leq
2c_{1}\left(1-\frac{1}{\sqrt{1+2c\kappa}}+O(\sqrt{\epsilon_{0}})+\epsilon_{1}\right)^{t}\left\lVert\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(1)})\\\ M_{*}\nabla F(x^{(0)})\\\ \end{array}\right]\right\rVert,$

with $M_{\star}=[\nabla^{2}F(x^{\star})]^{-1/2}$ and $\epsilon_{0}$,
$\epsilon_{1}$ defined in Theorem 4.

Proof  Consider $|{\mathcal{S}}|$ i.i.d random matrces
$H_{j}^{(t)},j=1,\dots,|{\mathcal{S}}|$ sampled uniformly. Then, we have
$\mathbb{E}(H_{j}^{(t)})=\nabla^{2}F(x^{(t)})$ for all
$j=1,\dots,|{\mathcal{S}}|$. By (10) and the positive semi-definite property
of $H_{j}^{(t)}$, we have $\lambda_{\max}(H_{j}^{(t)})\leq K$ and
$\lambda_{\min}(H_{j}^{(t)})\geq 0$.

We define random maxtrices $X_{j}=H_{j}^{(t)}-\nabla^{2}F(x^{(t)})$ for all
$j=1,\dots,|{\mathcal{S}}|$. We have $\mathbb{E}[X_{j}]=0$, $\|X_{j}\|\leq 2K$
and $\|X_{j}\|^{2}\leq 4K^{2}$. By the matrix Bernstein inequality, we have

$\mathbb{E}(\|H^{(t)}-\nabla^{2}F(x^{(t)})\|)\leq\sqrt{\frac{2\cdot 4K^{2}\log
2d}{|{\mathcal{S}}|}}+\frac{2K\log 2d}{3|{\mathcal{S}}|}.$

When the sample size $|{\mathcal{S}}|=O(c^{-2}K^{2}\log d)$,
$\|H^{(t)}-\nabla^{2}F(x^{(t)})\|\leq cK$ holds in expectation, that is

$\mathbb{E}\|H^{(t)}-\nabla^{2}F(x^{(t)})\|\leq c\nabla^{2}F(x^{(t)}).$

For notation convenience, we will omit superscript. Hence, in expectation, we
have

$\displaystyle|y^{T}(H-\nabla^{2}F(x))y|\leq
c\|\nabla^{2}F(x)\|\cdot\|y\|^{2}$ $\displaystyle\Rightarrow$ $\displaystyle
y^{T}HYy-c\|\nabla^{2}F(x)\|\cdot\|y\|^{2}\leq y^{T}\nabla^{2}F(x)y\leq
y^{T}Hy+c\|\nabla^{2}F(x)\|\cdot\|y\|^{2}$ $\displaystyle\Rightarrow$
$\displaystyle y^{T}\nabla^{2}F(x)y\leq y^{T}Hy\leq
y^{T}\nabla^{2}F(x)y+2c\|\nabla^{2}F(x)\|\cdot\|y\|^{2}$
$\displaystyle\Rightarrow$ $\displaystyle y^{T}[\nabla^{2}F(x)]^{-1}y\geq
y^{T}[H]^{-1}y\geq y^{T}\left(\nabla^{2}F(x)+2c\|\nabla^{2}F(x)\|\cdot
I\right)^{-1}y$

Furthermore, we have

$\displaystyle y^{T}\nabla^{2}F(x)y+2c\|\nabla^{2}F(x)\|\cdot\|y\|^{2}\leq
y^{T}\nabla^{2}F(x)y+\frac{2c\|\nabla^{2}F(x)\|}{\sigma_{\min}(\nabla^{2}F(x))}y^{T}\nabla^{2}F(x)y$
$\displaystyle\Rightarrow$
$\displaystyle\nabla^{2}F(x)+2c\|\nabla^{2}F(x)\|\cdot
I\preceq\left(1+2c\kappa\right)\nabla^{2}F(x)$ $\displaystyle\Rightarrow$
$\displaystyle\left(1-\frac{2c\kappa}{1+2c\kappa}\right)\nabla^{2}F(x)\preceq\left(\nabla^{2}F(x)+2c\|\nabla^{2}F(x)\|\cdot
I\right)^{-1}$

Thus, we have

$\left(1-\frac{2c\kappa}{1+2c\kappa}\right)\cdot\left[\nabla^{2}F(x^{\star})\right]^{-1}\preceq\mathbb{E}\left([H^{(t)}]^{-1}\right)\preceq\left[\nabla^{2}F(x^{\star})\right]^{-1}.$

Then final convergence property can be obtained by Theorem 4 using
$\pi=\frac{2c\kappa}{1+2c\kappa}$.

###  4.1 Fast Sub-problem Solver

In practice, we will sub-sample a small subset of samples, that is sample size
$s$ is much smaller than data dimension $d$. And Eqn. (9) becomes

$\displaystyle H^{(t)}=\tilde{B}^{T}\tilde{B}+\alpha^{(t)}I$ (13)

where $\tilde{B}=S^{(t)}B^{(t)}\in{\mathbb{R}}^{s\times d}$ with $s\ll d$.
Therefore, by Sherman-Woodbury identity formula, we can get the exact solution
of Eqn. (5) by

$\displaystyle p^{(t)}=[H^{(t)}]^{-1}\nabla
F(y^{(t)})=(\alpha^{(t)})^{-1}\nabla
F(y^{(t)})-\tilde{B}^{T}(\tilde{B}\tilde{B}^{T}+\alpha^{(t)}I)^{-1}\tilde{B}\nabla
F(y^{(t)})/\alpha^{(t)}$ (14)

And $p^{(t)}$ can be obtain in $O(ds^{2}+s^{3})$ time. Since Theorem 5 shows
that the direction vector can be an approximate solution of Eqn. (5), we can
get an inexact solution in a more efficient manner.

The main computational burden of Eqn. (14) is the matrix product of
$\tilde{B}$ and $\tilde{B}^{T}$. Rather than Sherman-Woodbury identity
formula, we can approximate the matrix inversion and get rid of the matrix
multiplication by conjugate gradient. The main computational cost of conjugate
gradient to compute $p^{(t)}$ is the matrix vector multiplication. This is
very suitable for sparse dataset. However, its computational complexity
depends on $\sqrt{\kappa(H^{(t)})}$ linearly. This may be expensive when the
condition number of $H^{(t)}$ is large.

Therefore, instead of using conjugate gradient method directly, we resort to
preconditioned conjugate gradient method (Algorithm 5) to obtain an
approximation of $(\tilde{B}\tilde{B}^{T}+\alpha I)^{-1}g$ in Eqn. (14) with
$g=\tilde{B}\cdot\alpha^{-1}\nabla F(y^{(t)})$. Because of
$\tilde{B}\in{\mathbb{R}}^{s\times d}$ with $s\ll d$, we can use sketching
tools to construct a good preconditioner very efficiently. We describe the
detailed algorithm in Algorithm 2. And we have the following result.

###### Theorem 7

Let the approximate Hessian $H^{(t)}$ be of form Eqn. (13) with
$\tilde{B}\in{\mathbb{R}}^{s\times d}$ and $s\ll d$. Set iteration number $T$
as

$T=\log\left(\frac{\alpha^{-\frac{3}{2}}\left\lVert
H^{(t)}\right\rVert^{\frac{3}{2}}c_{1}\sqrt{\kappa}}{\epsilon_{1}}\right).$

where $c_{1}$ and $\kappa$ are defined in Theorem 4. $p^{(t)}$ is returned
from Algorithm 2, then we have

$\left\lVert H^{(t)}p^{(t)}-\nabla
F(y^{(t)})\right\rVert\leq\frac{\epsilon_{1}}{12c_{1}\sqrt{\kappa}}\left\lVert\nabla
F(y^{(t)})\right\rVert.$

And the computation complexity of computing $p^{(t)}$ is $O(Tds+s^{3})$.

Proof  First, by the property of subspace embedding property (Definition 2)
and $G$ is a $1/3$-subspace embedding matrix w.r.t $\tilde{B}$, we have

$\left(1-\frac{1}{3}\right)\tilde{B}\tilde{B}^{T}\preceq\tilde{B}GG^{T}\tilde{B}^{T}\preceq\left(1+\frac{1}{3}\right)\tilde{B}\tilde{B}^{T}.$

With a simple transformation, we obtain

$\left(1-\frac{1}{2}\right)\tilde{B}GG^{T}\tilde{B}^{T}\preceq\tilde{B}\tilde{B}^{T}\preceq\left(1+\frac{1}{2}\right)\tilde{B}GG^{T}\tilde{B}^{T}.$

Thus sketced matrix $P$ is a good preconditioner. And by the convergence
property of preconditioned conjugate gradient method (Lemma 8), after
Algorithm 5 runing $T$ iterations, we have

$\displaystyle\left\lVert
A^{-1}g-q\right\rVert_{A}\leq\left(\frac{1}{2}\right)^{T}\left\lVert
A^{-1}g\right\rVert_{A}$ $\displaystyle\Rightarrow$
$\displaystyle\sqrt{\lambda_{\min}(A)}\cdot\left\lVert
A^{-1}g-q\right\rVert\leq\left(\frac{1}{2}\right)^{T}\sqrt{g^{T}A^{-1}g}$

where $A=\tilde{B}\tilde{B}^{T}+\alpha I$ and
$g=\tilde{B}\cdot\alpha^{-1}\nabla F(y^{(t)})$. We also have

$\displaystyle g^{T}A^{-1}g=\alpha^{-2}\cdot\nabla
F(y^{(t)})^{T}\tilde{B}^{T}(\tilde{B}\tilde{B}^{T}+\alpha
I)^{-1}\tilde{B}\nabla F(y^{(t)})\leq\alpha^{-2}\lVert\nabla
F(y^{(t)})\rVert^{2}.$

The second inequality is because $\tilde{B}^{T}(\tilde{B}\tilde{B}^{T}+\alpha
I)^{-1}\tilde{B}$ is positive definite and its largest eigenvalue is smaller
than $1$. Thus, we have

$\displaystyle\left\lVert
A^{-1}g-q\right\rVert\leq\left(\frac{1}{2}\right)^{T}\alpha^{-\frac{3}{2}}\left\lVert\nabla
F(y^{(t)})\right\rVert.$

Furthermore, we have

$\displaystyle\left\lVert H^{(t)}p^{(t)}-\nabla F(y^{(t)})\right\rVert\leq$
$\displaystyle\left\lVert H^{(t)}\right\rVert\cdot\left\lVert
p^{(t)}-[H^{(t)}]^{-1}\nabla F(y^{(t)})\right\rVert$ $\displaystyle=$
$\displaystyle\left\lVert H^{(t)}\right\rVert\cdot\left\lVert\alpha^{-1}\nabla
F(y^{(t)})-\tilde{B}^{T}q-\left(\alpha^{-1}\nabla
F(y^{(t)})-\tilde{B}^{T}A^{-1}g\right)\right\rVert$ $\displaystyle\leq$
$\displaystyle\left(\frac{1}{2}\right)^{T}\alpha^{-\frac{3}{2}}\left\lVert
H^{(t)}\right\rVert\cdot\left\lVert\tilde{B}\right\rVert\cdot\left\lVert\nabla
F(y^{(t)})\right\rVert$ $\displaystyle\leq$
$\displaystyle\left(\frac{1}{2}\right)^{T}\alpha^{-\frac{3}{2}}\left\lVert
H^{(t)}\right\rVert^{\frac{3}{2}}\cdot\left\lVert\nabla
F(y^{(t)})\right\rVert$

Since we set iteration number $T$ to be

$T=\log\left(\frac{\alpha^{-\frac{3}{2}}\left\lVert
H^{(t)}\right\rVert^{\frac{3}{2}}c_{1}\sqrt{\kappa}}{\epsilon_{1}}\right),$

we obtain the result.

Because the main operation in preconditioned conjugate gradient is matrix-
vector production, the total cost is $O(Tds+s^{3})$.

From Theorem 7, we can see that iteration number $T$ is only of order
$O(\log\kappa)$. Therefore the time of computing a $p^{(t)}$ by Algorithm 2 is
cheap than by Eqn. (14) when $O(\log\kappa)<s$. And this is common in
practice.

Algorithm 2 Fast sub-problem solver.

1: Input: Matrix $\tilde{B}$, regularizer $\alpha$, gradient $\nabla
F(y^{(t)})$, and iteration number $T$;

2: Compute $g=\tilde{B}\cdot\alpha^{-1}\nabla F(y^{(t)})$. Construct a
$1/3$-subspace embedding matrix $G$ w.r.t $\tilde{B}$ and compute
preconditioner matrix $P=\tilde{B}GG^{T}\tilde{B}$.

3: Using $A=\tilde{B}\tilde{B}^{T}+\alpha I$, $T$, and $P$ as input of
Algorithm 5, we get a approximate solution vector $q$ by Algorithm 5.

4: Output: $p^{(t)}=\alpha^{-1}\nabla F(y^{(t)})-\tilde{B}^{T}q$.

##  5 Experiments

In Section 3.1, we have shown that Nesterov’s acceleration technique can
improve the performance of approximate Newton method theoretically. In this
section, we will validate our theory empirically. In particular, we first
compare accelerated regularized sub-sampled Newton (Algorithm 1 refered as
ARSSN) with regularized sub-sampled Newton (Algorithm 3 refered as RSSN) on
the ridge regression whose objective function is a quadratic function to
validate the theoretical analysis in Section 3.1. Then we conduct more
experiments on a popular machine learning problem called Ridge Logistic
Regression, and compare accelerated regularized sub-sampled Newton with other
classical algorithms.

(a) $|{\mathcal{S}}|=1\%$ (b) $|{\mathcal{S}}|=5\%$

(c) $|{\mathcal{S}}|=10\%$

Figure 1: Experiment on the least square regression with $|{\mathcal{S}}|$
being different sample size

###  5.1 Experiments on the Ridge Regression

The ridge regression is defined as follows

$\displaystyle F(x)=\|Ax-b\|^{2}+\lambda\|x\|^{2},$ (15)

where $\lambda$ is a regularizer controls the condition number of the Hessian.
In our experiment In our experiments, we choose dataset ‘gisette’ just as
depicted in Table 1. And we set the regularizer $\lambda=1$.

In experiments, we set the sample size $|{\mathcal{S}}|$ to be $1\%n$, $5\%n$
and $10\%n$. The regularizer $\alpha$ of Algorithm 1 is properly chosen
according to $|{\mathcal{S}}|$. ARSSN and RSSN share the same
$|{\mathcal{S}}|$. And the acceleration parameter $\theta$ of is fixed and
appropriately selected. We report the experiments result in Figure 1.

From Figure 1, we can see that ARSSN and RSSN have significant difference in
convergence rate and ARSSN is much faster. This validates the analysis in
Section 3.1. Besides, we can also observe that ARSSN runs faster as sample
size $|{\mathcal{S}}|$ increases. When $|{\mathcal{S}}|=10\%n$, ARSSN takes
only about $3000$ iterations to achieve an $10^{-14}$ error while it needs
about $6000$ iterations to achieve the same precision when
$|{\mathcal{S}}|=5\%n$.

Table 1: Datasets summary(sparsity$=\frac{\\#\text{Non-Zero Entries}}{n\times d}$) Dataset |  $n$ |  $d$ |  sparsity |  source
---|---|---|---|---
gisette |  $5,000$ |  $6,000$ |  dense |  libsvm dataset
sido0 |  $12,678$ |  $4,932$ |  dense |  [9]
svhn |  $19,082$ |  $3,072$ |  dense |  libsvm dataset
rcv1 |  $20,242$ |  $47,236$ |  $0.16\%$ |  libsvm dataset
real-sim |  $72,309$ |  $20,958$ |  $0.24\%$ |  libsvm dataset
avazu |  $2,085,163$ |  $999,975$ |  $0.0015\%$ |  libsvm dataset

###  5.2 Experiments on the Ridge Logistic Regression

We conduct experiments on the Ridge Logistic Regression problem whose
objective is

$F(x)=\frac{1}{n}\sum_{i=1}^{n}\log[1+\exp(-b_{i}\langle
a_{i},x\rangle)]+\frac{\lambda}{2}\|x\|^{2},$ (16)

where $a_{i}\in{\mathbb{R}}^{d}$ is the $i$-th input vector, and
$b_{i}\in\\{-1,1\\}$ is the corresponding label.

We conduct our experiments on six datasets: ‘gisette’, ‘protein’, ‘svhn’,
‘rcv1’, ‘sido0’, and ‘real-sim’. The first three datasets are dense and the
last three ones are sparse. We give the detailed description of the datasets
in Table 1. Notice that the size and dimension of dataset are close to each
other, so the sketch Newton method [19, 25] can not be used. In our
experiments, we try different settings of the regularizer $\lambda$, as $1/n$,
$10^{-1}/n$, and $10^{-2}/n$ to represent different levels of regularization.

We compare Algorithm 1 (ARSSN) with RSSN (Algorithm 3), AGD ([17]) and SVRG
([11]) which are classical and popular optimization methods in machine
learning.

In our experiments, the sample size $|{\mathcal{S}}|$ and regularizer $\alpha$
of ARSSN are chosen according to Theorem 5 while those of RSSN are chosen by
Theorem 9. For a fixed $|{\mathcal{S}}|$, a proper $\alpha$ can be found after
several tries.

In our experiment, the sub-sampled Hessian $H^{(t)}$ constructed in Algorithm
1 can be written as

$H^{(t)}=\tilde{A}^{T}\tilde{A}+(\alpha+\lambda)I,$

where $\tilde{A}\in{\mathbb{R}}^{\ell\times d}$, where $\ell<n$. We can resort
to Woodbury’ identity to compute the inverse of $H^{(t)}$. If $\tilde{A}$ is
sparse, we use conjugate gradient (Algorithm 4 in Appendix) to obtain an
approximation of $[H^{(t)}]^{-1}\nabla F(x^{(t)})$ which exploits the sparsity
of $\tilde{A}$. In our experiments on sparse datasets, we set $tol=0.1\|\nabla
F(x^{(t)})\|$ for conjugate gradient (Algorithm 4 ).

For the acceleration parameter $\theta$, it is hard to get the best value for
ARSSN just like AGD. However, our theoretical analysis implies that for large
sample size $|{\mathcal{S}}|$, a small $\theta$ should chosen. In our
experiments, we set $\theta^{(t)}=\frac{t}{t+16}$ in Algorithm 1 for the dense
datasets and $\theta^{(t)}=\frac{t}{t+30}$ for the sparse datasets. We set
$x^{(0)}=0$ for all the datasets and all the algorithms.

We report our result in Figure 2, 3, 4, 5, 6 and 7. We can see that ARSSN
converges much faster than RSSN when these two algorithms have the same sample
size. This shows Nesterov’s acceleration technique can promote the performance
of regularized sub-sampled Newton effectively. We can also observe that ARSSN
outperforms AGD significantly even when the sample size ${\mathcal{S}}$ is
$1\%n$ or even less. This validates the fact that adding curvature information
is an effective way to improve the ability of accelerated gradient descent.

Compared with SVRG, we can see that ARSSN also has better performance.
Specifically, ARSSN performs much better than SVRG when $\lambda$ is small
like $\lambda=10^{-2}/n$. When $\lambda=10^{-1}$, ARSSN has comparable
performance with SVRG on ‘sido0’, ‘rcv1’, ‘real-sim’, and ‘avazu’ while ARSSN
performs much better on ‘gisette’ and ‘svhn’. This means that ARSSN is an
efficient algorithm compared with classical stochastic first order method.

The experiments also reveal the fact that ARSSN has great advantages over
other algorithms when the problem is ill-conditioned. On ‘gisette’, ‘sido0’,
‘svhn’, other algorithms have very poor performance on the case that
$\lambda=10^{-2}/n$. But ARSSN shows good convergence property.

###  5.3 Conclusion of Empirical Study

The above experiments show that Nesterov’s acceleration is an effective way to
promote the convergence rate of approximate Newton methods. The experiments
also show that adding some curvature information always help AGD to obtain a
faster convergence rate.

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 2: Experiment on ‘gisette’

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 3: Experiment on ‘sido0’

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 4: Experiment on ‘svhn’

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 5: Experiment on ‘rcv1’

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 6: Experiment on ‘real-sim’

(a) $\lambda=1/n$ (b) $\lambda=10^{-1}/n$

(c) $\lambda=10^{-2}/n$

Figure 7: Experiment on ‘avazu’

##  6 Conclusion

In this paper, we have exploited Nesterov’s acceleration technique to promote
the performance of second-order methods and propose Accelerated Regularized
Sub-sample Newton. We have presented the theoretical analysis on the
convergence properties of accelerated second-order methods, showing that
accelerated approximate Newton has higher convergence rate, especially when
the approximate Hessian is not a good approximation. Based on our theory, we
have developed ARSSN. Our experiments have shown that our ARSSN performs much
better than the conventional RSSN, which meets our theory well. ARSSN also has
several advantages over other classical algorithms, demonstrating the
efficiency of accelerated second-order methods.

## References

  * Allen-Zhu [2016] Zeyuan Allen-Zhu.  Katyusha: The first direct acceleration of stochastic gradient methods.  _arXiv preprint arXiv:1603.05953_ , 2016.
  * Beck and Teboulle [2009] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM journal on imaging sciences_ , 2(1):183–202, 2009.
  * Byrd et al. [2011] Richard H Byrd, Gillian M Chin, Will Neveitt, and Jorge Nocedal.  On the use of stochastic hessian information in optimization methods for machine learning.  _SIAM Journal on Optimization_ , 21(3):977–995, 2011.
  * Clarkson and Woodruff [2013] Kenneth L Clarkson and David P Woodruff.  Low rank approximation and regression in input sparsity time.  In _Proceedings of the forty-fifth annual ACM symposium on Theory of computing_ , pages 81–90. ACM, 2013.
  * Cotter et al. [2011] Andrew Cotter, Ohad Shamir, Nati Srebro, and Karthik Sridharan.  Better mini-batch algorithms via accelerated gradient methods.  In _Advances in neural information processing systems_ , pages 1647–1655, 2011.
  * Drineas et al. [2006] Petros Drineas, Michael W Mahoney, and S Muthukrishnan.  Sampling algorithms for l 2 regression and applications.  In _Proceedings of the seventeenth annual ACM-SIAM symposium on Discrete algorithm_ , pages 1127–1136. Society for Industrial and Applied Mathematics, 2006.
  * Erdogdu and Montanari [2015] Murat A Erdogdu and Andrea Montanari.  Convergence rates of sub-sampled newton methods.  In _Advances in Neural Information Processing Systems_ , pages 3034–3042, 2015.
  * Golub and Van Loan [2012] Gene H Golub and Charles F Van Loan.  _Matrix computations_ , volume 3.  JHU Press, 2012.
  * [9] I. Guyon.  Sido: A phamacology dataset.  URL http://www.causality.inf.ethz.ch/data/SIDO.html.
  * Halko et al. [2011] N Halko, P G Martinsson, and J A Tropp.  Finding Structure with Randomness : Probabilistic Algorithms for Matrix Decompositions.  _SIAM Review_ , 53(2):217–288, 2011.
  * Johnson and Zhang [2013] Rie Johnson and Tong Zhang.  Accelerating stochastic gradient descent using predictive variance reduction.  In _Advances in Neural Information Processing Systems_ , pages 315–323, 2013.
  * Johnson and Lindenstrauss [1984] William B. Johnson and Joram Lindenstrauss.  Extensions of Lipschitz mappings into a Hilbert space.  _Contemporary mathematics_ , 26(189-206), 1984.
  * Lan and Zhou [2015] Guanghui Lan and Yi Zhou.  An optimal randomized incremental gradient method.  _arXiv preprint arXiv:1507.02000_ , 2015.
  * Li et al. [2014] Mu Li, Tong Zhang, Yuqiang Chen, and Alexander J Smola.  Efficient mini-batch training for stochastic optimization.  In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 661–670. ACM, 2014.
  * Meng and Mahoney [2013] Xiangrui Meng and Michael W Mahoney.  Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression.  In _Proceedings of the forty-fifth annual ACM symposium on Theory of computing_ , pages 91–100. ACM, 2013.
  * Nelson and Nguyên [2013] Jelani Nelson and Huy L Nguyên.  Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In _Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on_ , pages 117–126. IEEE, 2013.
  * Nesterov [1983] Yurii Nesterov.  A method of solving a convex programming problem with convergence rate o (1/k2).  In _Soviet Mathematics Doklady_ , volume 27, pages 372–376, 1983\.
  * Nesterov et al. [2007] Yurii Nesterov et al.  Gradient methods for minimizing composite objective function, 2007.
  * Pilanci and Wainwright [2017] Mert Pilanci and Martin J. Wainwright.  Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence.  _SIAM Journal on Optimization_ , 27(1):205–245, 2017.
  * Robbins and Monro [1951] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  _The annals of mathematical statistics_ , pages 400–407, 1951.
  * Roosta-Khorasani and Mahoney [2016] Farbod Roosta-Khorasani and Michael W Mahoney.  Sub-sampled newton methods ii: Local convergence rates.  _arXiv preprint arXiv:1601.04738_ , 2016.
  * Roux et al. [2012] Nicolas L Roux, Mark Schmidt, and Francis R Bach.  A stochastic gradient method with an exponential convergence _rate for finite training sets.  In _Advances in Neural Information Processing Systems_ , pages 2663–2671, 2012.
  * Schmidt et al. [2013] Mark Schmidt, Nicolas Le Roux, and Francis Bach.  Minimizing finite sums with the stochastic average gradient.  _arXiv preprint arXiv:1309.2388_ , 2013.
  * Woodruff [2014] David P Woodruff.  Sketching as a tool for numerical linear algebra.  _Foundations and Trends® in Theoretical Computer Science_ , 10(1–2):1–157, 2014.
  * Xu et al. [2016] Peng Xu, Jiyan Yang, Farbod Roosta-Khorasani, Christopher Ré, and Michael W Mahoney.  Sub-sampled newton methods with non-uniform sampling.  In _Advances in Neural Information Processing Systems_ , pages 3000–3008, 2016.
  * Ye et al. [2017] Haishan Ye, Luo Luo, and Zhihua Zhang.  Approximate newton methods and their local convergence.  In _International Conference on Machine Learning_ , pages 3931–3939, 2017.
  * Zhang et al. [2013] Lijun Zhang, Mehrdad Mahdavi, and Rong Jin.  Linear convergence with condition number independent access of full gradients.  In _Advance in Neural Information Processing Systems 26 (NIPS)_ , pages 980–988, 2013.

##  A Conjugate Gradient Descent

Conjugate Gradient is a classical method to solve the linear system of
equation

$Ax=b,$

where $A$ is a $d\times d$ symmetric positive definite matrix. Algorithm 4
gives the detailed implementation to solve above linear system. Conjugate
gradient has the following convergence properties [8]:

$\|x_{k}-x_{*}\|_{A}\leq
2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^{k}\|x_{0}-x_{*}\|_{A}.$

As we can see, the convergence behavior of conjugate gradient method depends
on the condition of $A$. Hence, it suffers from poor convergence rate When $A$
is ill-conditioned. Precondition method is an important way to improve the
convergence properties. In the following Lemma, we give the sufficient
condition of a preconditioner.

###### Lemma 8

If $A$ and $B$ are $d\times d$ symmetric positive matrices, and
$(1-\varepsilon)B\preceq A\preceq(1+\varepsilon)B$, where $0<\epsilon_{0}<1$,
then for the optimization problem $\min_{x}\|Ax-b\|$, we have

$\|x_{1}-x_{*}\|_{A}\leq\varepsilon\|x_{*}\|_{A},$

where $x_{*}=A^{-1}b$ and $x_{1}=B^{-1}b$. Besides, if we set
$r_{k}=Ax_{k}-b$, $x_{kr}=B^{-1}r_{k}$ and $x_{k+1}=x_{k}+x_{kr}$, then
$\|x_{k+1}-x_{*}\|_{A}\leq\epsilon_{0}\|x_{k}-x_{*}\|_{A}$.

Proof  Because $A\preceq(1+\varepsilon)B$, we have

$\displaystyle\lambda_{\max}(B^{-1}A)=$
$\displaystyle\lambda_{\max}(B^{-1/2}AB^{-1/2})$ $\displaystyle=$
$\displaystyle\max_{x\neq 0}\frac{x^{T}B^{-1/2}AB^{-1/2}x}{x^{T}x}$
$\displaystyle=$ $\displaystyle\max_{y\neq 0}\frac{y^{T}Ay}{y^{T}By}$
$\displaystyle\leq$ $\displaystyle 1+\varepsilon,$

where the last equality is by setting $y=B^{-1/2}x$. Similarly, we have
$\lambda_{\min}(B^{-1}A)\geq 1-\varepsilon$.

Since $B^{-1}A$ and $A^{1/2}B^{-1}A^{1/2}$ are similar, the eigenvalues of
$A^{1/2}B^{-1}A^{1/2}$ are all between $1-\varepsilon$ and $1+\varepsilon$.
Therefore, we have

$\|A^{1/2}B^{-1}A^{1/2}-I\|\leq\varepsilon,$

where $I$ is the $d\times d$ identity matrix.

Hence, we have

$\displaystyle\|B^{-1}b-x_{*}\|_{A}=$
$\displaystyle\|A^{1/2}(B^{-1}Ax_{*}-x_{*})\|$ $\displaystyle=$
$\displaystyle\|A^{1/2}B^{-1}A^{1/2}(A^{1/2}x_{*})-A^{1/2}x_{*}\|$
$\displaystyle\leq$
$\displaystyle\|A^{1/2}B^{-1}A^{1/2}-I\|\cdot\|A^{1/2}x_{*}\|$
$\displaystyle=$ $\displaystyle\varepsilon\|x_{*}\|_{A}.$

For the convergence property, we have

$\displaystyle\|x_{k+1}-x_{*}\|_{A}=$
$\displaystyle\|B^{-1}(Ax_{k}-b)+x_{k}-A^{-1}b\|_{A}$ $\displaystyle=$
$\displaystyle\|B^{-1}(Ax_{k}-b)-A^{-1}(Ax_{k}-b)\|_{A}$ $\displaystyle\leq$
$\displaystyle\epsilon_{0}\|x_{k}-A^{-1}b\|_{A}$ $\displaystyle=$
$\displaystyle\varepsilon\|x_{k}-x_{*}\|_{A}.$

The inequality is due to the property $\|B^{-1}b-x_{*}\|_{A}\leq\|x_{*}\|_{A}$
and $b$ replaced with $Ax_{k}-b$ and $x_{*}$ replaced with $A^{-1}(Ax_{k}-b)$.

Let $A$ and $b$ be replaced by $B^{-1/2}AB^{-1/2}$ and $B^{-1/2}b$
respectively. And with some transformation, the iterations described in Lemma
8 can be transformed into preconditioned conjugate gradient method described
in Algorithm 5. From Lemma 8, we can see that preconditioned conjugate
gradient converges linearly with a constant rate using a good preconditioner.

##  B Regularized Sub-sampled Newton

The regularized Sub-sampled Newton method is depicted in Algorithm 3, and we
now give its local convergence properties in the following theorem [26].

###### Theorem 9

Let $F(x)$ satisfy Assumption 1 and 2. Assume Eqns. (10) and (11) hold, and
let $0<\delta<1$, $0\leq\epsilon_{1}<1$ and $0<\alpha$ be given. Assume
$\beta$ is a constant such that $0<\beta<\alpha+\frac{\sigma}{2}$, the
subsampled size $|{\mathcal{S}}|$ satisfies
$|{\mathcal{S}}|\geq\frac{16K^{2}\log(2d/\delta)}{\beta^{2}}$, and $H^{(t)}$
is constructed as in Algorithm 3. Define

$\displaystyle\epsilon_{0}=\max\left(\frac{\beta-\alpha}{\sigma+\alpha-\beta},\frac{\alpha+\beta}{\sigma+\alpha+\beta}\right),$

which implies that $0<\epsilon_{0}<1$. And we define
$\|x\|_{M^{*}}=\|[M^{*}]^{-\frac{1}{2}}x\|$. Then Algorithm 3 has the
following convergence properties:

  1. 1.

There exists a sufficient small value $\gamma$, $0<\nu(t)<1$, and
$0<\eta(t)<1$ such that when $\|x^{(t)}-x^{*}\|\leq\gamma$, each iteration
satisfies

$\|\nabla
F(x^{(t+1)})\|_{M^{*}}\leq\left(\epsilon_{0}+\frac{2\eta(t)}{1-\epsilon_{0}}\right)\frac{1+\nu(t)}{1-\nu(t)}\|\nabla
F(x^{(t)})\|_{M^{*}}.$

Besides, $\nu(t)$ and $\eta(t)$ will go to $0$ as $x^{(t)}$ goes to $x^{*}$.

  2. 2.

If $\nabla^{2}F(x^{(t)})$ is also Lipschitz continuous with parameter
$\hat{L}$ and $x^{(t)}$ satisfies

$\|x^{(t)}-x^{*}\|\leq\frac{\mu}{\hat{L}\kappa}\nu(t),$

where $0<\nu(t)<1$, then it holds that

$\displaystyle\|\nabla F(x^{(t+1)})\|_{M^{*}}\leq$
$\displaystyle\epsilon_{0}\frac{1+\nu(t)}{1-\nu(t)}\|\nabla
F(x^{(t)})\|_{M^{*}}+\frac{2}{(1-\epsilon_{0})^{2}}\frac{\hat{L}\kappa}{\mu\sqrt{\mu}}\frac{(1+\nu(t))^{2}}{1-\nu(t)}\|\nabla
F(x^{(t)})\|_{M^{*}}^{2}.$

Algorithm 3 Regularized Sub-sample Newton (RSSN).

1: Input: $x^{(0)}$, $0<\delta<1$, regularizer parameter $\alpha$, sample size
$|{\mathcal{S}}|$ ;

2: for $t=0,1,\dots$ until termination do

3: Select a sample set ${\mathcal{S}}$, of size $|{\mathcal{S}}|$ and
$H^{(t)}=\frac{1}{|{\mathcal{S}}|}\sum_{j\in\mathcal{S}}\nabla^{2}f_{j}(x^{(t)})+\alpha
I$;

4: Update $x^{(t+1)}=x^{(t)}-\left[H^{(t)}\right]^{-1}\nabla F(x^{(t)})$;

5: end for

##  C Proof of Lemma 3

Proof  By Taylor’s theorem, we have

$\displaystyle\nabla F(x^{(t+1)})$ $\displaystyle=$ $\displaystyle\nabla
F(y^{(t+1)})+\nabla^{2}F(y^{(t+1)})(-p^{(t+1)})+o(p^{(t+1)})$ $\displaystyle=$
$\displaystyle\nabla F(y^{(t+1)})-\nabla^{2}F(y^{(t+1)})[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})+o(p^{(t+1)})$ $\displaystyle=$ $\displaystyle\nabla
F(y^{(t+1)})-\nabla^{2}F(x^{*})[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})-(\nabla^{2}F(x^{*})-\nabla^{2}F(y^{(t+1)}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})+o(p^{(t+1)})$ $\displaystyle=$
$\displaystyle\left[\nabla^{2}F(x^{*})\right]^{\frac{1}{2}}\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(y^{(t+1)})$
$\displaystyle+(\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})+o(p^{(t+1)}).$

For $\nabla F(y^{(t+1)})$, we have

$\displaystyle\nabla F(y^{(t+1)})=$ $\displaystyle\nabla F(x^{(t)}+\theta
s^{(t)})$ $\displaystyle=$ $\displaystyle\nabla
F(x^{(t)})+\theta\nabla^{2}F(x^{(t)})(s^{(t)})+o(s^{(t)})$ $\displaystyle=$
$\displaystyle\nabla F(x^{(t)})+\theta(\nabla F(x^{(t)})-\nabla
F(x^{(t-1)}))+o(s^{(t)})$ $\displaystyle=$ $\displaystyle(1+\theta)\nabla
F(x^{(t)})-\theta\nabla F(x^{(t-1)})+o(s^{(t)}).$

Besides, we have

$\displaystyle o(s^{(t)})=$ $\displaystyle
o(x^{(t)}-x^{(t-1)})=o(x^{(t)}-x^{*}-(x^{(t-1)}-x^{*}))$ $\displaystyle=$
$\displaystyle o(\nabla F(x^{(t)}))+o(\nabla F(x^{(t-1)})),$

where the last equality is because $\nabla F(x)$ is $L$-Lipschitz continuous.
Hence, it holds that

$\displaystyle\nabla F(y^{(t+1)})=(1+\theta)\nabla F(x^{(t)})-\theta\nabla
F(x^{(t-1)})+o(\nabla F(x^{(t)})+\nabla F(x^{(t-1)}))$ (17)

For $(\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})$, we show that it is of order $o(\nabla F(y^{(t+1)}))$ as
follows. First, if $\nabla^{2}F(x)$ is not Lipschitz continuous, then there
exists a $\gamma$ such that when $\|y-x^{*}\|\leq\gamma$, it holds that

$\|\nabla^{2}F(y)-\nabla^{2}F(x^{*})\|=o(1).$

Such $\gamma$ exists because $\nabla^{2}F(x)$ is continuous near optimal point
$x^{*}$. Hence,
$(\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})$ is of order $o(\nabla F(y^{(t+1)}))$ when $y^{(t+1)}$ is
sufficient close to $x^{*}$.

If $\nabla^{2}F(x)$ is $\hat{L}$-Lipschitz continuous and $F(x)$ is
$\mu$-strongly convex, then we have

$\|\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*})\|\leq\hat{L}\|y^{(t+1)}-x^{*}\|\leq\frac{\hat{L}}{\mu}\|\nabla
F(y^{(t+1)})\|.$

Then, it holds that

$\|(\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})\|=O(\|\nabla F(y^{(t+1)})\|^{2}).$

Besides, because of $p^{(t+1)}=[H^{(t+1)}]^{-1}\nabla F(y^{(t+1)})$, we have

$o(p^{(t+1)})=o(\nabla F(y^{(t+1)})).$

Combining Eqn. (17), we have

$\displaystyle(\nabla^{2}F(y^{(t+1)})-\nabla^{2}F(x^{*}))[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})=o(\nabla F(x^{(t)}))+o(\nabla F(x^{(t-1)}))$

and

$\displaystyle o(p^{(t+1)})=o(\nabla F(x^{(t)}))+o(\nabla F(x^{(t-1)})).$

Hence, we have the following result

$\displaystyle[\nabla^{2}F(x^{*})]^{-\frac{1}{2}}\nabla F(x^{(t+1)})$
$\displaystyle=$
$\displaystyle\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)\left((1+\theta)\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t)})-\theta\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(x^{(t-1)})\right)$ $\displaystyle+o(\nabla F(x^{(t)}))+o(\nabla
F(x^{(t-1)})).$

Algorithm 4 Conjugate Gradient Descent Method.

1: Input: $A$, $b$,$x_{0}$, and $tol$;

2: Set $r_{0}=Ax_{0}-b$, $p_{0}=-r_{0}$, $k=0$;

3: while $\|r_{k}\|>tol$ do

4: Calculate $\alpha_{k}=\frac{r_{k}^{T}r_{k}}{p_{k}^{T}Ap_{k}}$;

5: Calculate $x_{k+1}=x_{k}+\alpha_{k}p_{k}$ and
$r_{k+1}=r_{k}+\alpha_{k}Ap_{k}$;

6: Calculate $\beta_{k+1}=\frac{r_{k+1}^{T}r_{k+1}}{r_{k}^{T}r_{k}}$ and
$p_{k+1}=-r_{k+1}+\beta_{k+1}p_{k}$;

7: $k=k+1$;

8: end while

9: Output: $x_{k}$.

Algorithm 5 Preconditioned Conjugate Gradient Descent Method.

1: Input: $x_{0}$, iteration number $T$, and preconditioner $P$;

2: Set $r_{0}=Ax_{0}-b$, solve $Py_{0}=r_{0}$, set $p_{0}=-y_{0}$, $k=0$;

3: while $k<T$ do

4: Calculate $\alpha_{k}=\frac{r_{k}^{T}r_{k}}{p_{k}^{T}Ap_{k}}$;

5: Calculate $x_{k+1}=x_{k}+\alpha_{k}p_{k}$ and
$r_{k+1}=r_{k}+\alpha_{k}Ap_{k}$;

6: Solve $Py_{k+1}=r_{k+1}$;

7: Calculate $\beta_{k+1}=\frac{r_{k+1}^{T}y_{k+1}}{r_{k}^{T}y_{k}}$ and
$p_{k+1}=-y_{k+1}+\beta_{k+1}p_{k}$;

8: $k=k+1$;

9: end while

10: Output: $x_{k}$.

##  D Proof of Theorem 4

Proof  of Theorem 4 By Taylor’s theorem and the property of qudratic function,
we have

$\displaystyle\nabla F(x^{(t+1)})$ $\displaystyle=$ $\displaystyle\nabla
F(y^{(t+1)})+\nabla^{2}F(y^{(t+1)})(-p^{(t+1)})$ $\displaystyle=$
$\displaystyle\nabla F(y^{(t+1)})-\nabla^{2}F(y^{(t+1)})[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})$ $\displaystyle=$ $\displaystyle\nabla
F(y^{(t+1)})-\nabla^{2}F(x^{*})[H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})+\nabla^{2}F(x^{*})\left([H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})-p^{(t)}\right)$ $\displaystyle=$
$\displaystyle\left[\nabla^{2}F(x^{*})\right]^{\frac{1}{2}}\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}[H^{(t+1)}]^{-1}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)\left[\nabla^{2}F(x^{*})\right]^{-\frac{1}{2}}\nabla
F(y^{(t+1)})$ $\displaystyle+\nabla^{2}F(x^{*})\left([H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})-p^{(t)}\right).$

For $\nabla F(y^{(t+1)})$, we have

$\displaystyle\nabla F(y^{(t+1)})=$ $\displaystyle\nabla F(x^{(t)}+\theta
s^{(t)})$ $\displaystyle=$ $\displaystyle\nabla
F(x^{(t)})+\theta\nabla^{2}F(x^{(t)})(s^{(t)})+o(s^{(t)})$ $\displaystyle=$
$\displaystyle\nabla F(x^{(t)})+\theta(\nabla F(x^{(t)})-\nabla F(x^{(t-1)}))$
$\displaystyle=$ $\displaystyle(1+\theta)\nabla F(x^{(t)})-\theta\nabla
F(x^{(t-1)}).$

For notational convinience, we use $M_{*}$ to denote
$[\nabla^{2}F(x^{*})]^{-\frac{1}{2}}$, we have

$\displaystyle\mathbb{E}\left(M_{*}\nabla F(x^{(t+1)})\right)$
$\displaystyle=$
$\displaystyle\left(I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(t+1)}]^{-1}\right)[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\right)\left((1+\theta)M_{*}\nabla
F(x^{(t)})-\theta M_{*}\nabla F(x^{(t-1)})\right)$
$\displaystyle+[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})-p^{(t)}\right).$ (18)

Because sketching matrices share the same distribution,
$\mathbb{E}\left([H^{(t+1)}]^{-1}\right)$ is a constant matrix. we use the
following notation for convenience

$\displaystyle
K=I-[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(t+1)}]^{-1}\right)[\nabla^{2}F(x^{*})]^{\frac{1}{2}}.$

And $K$ has the following spectral decomposition

$\displaystyle K=U\Lambda U^{T},$

with $\lambda_{1}\geq\lambda_{2}\dots\lambda_{d}$ We can reformulate the Eqn.
(18) as

$\displaystyle\left[\begin{array}[]{c}\mathbb{E}\left(M_{*}\nabla
F(x^{(t+1)})\right)\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]=$
$\displaystyle\left[\begin{array}[]{cc}(1+\theta)K,&-\theta K\\\ I,&0\\\
\end{array}\right]\cdot\left[\begin{array}[]{c}M_{*}\nabla F(x^{(t)})\\\
M_{*}\nabla F(x^{(t-1)})\\\ \end{array}\right]$
$\displaystyle+[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(t+1)}]^{-1}\nabla
F(y^{(t+1)})-p^{(t)}\right).$

Define the matrix

$\displaystyle T=\left[\begin{array}[]{cc}(1+\theta)K,&-\theta K\\\ I,&0\\\
\end{array}\right].$

Thus, we obtain

$\displaystyle\left[\begin{array}[]{c}\mathbb{E}\left(M_{*}\nabla
F(x^{(t+1)})\right)\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]$
$\displaystyle=$ $\displaystyle T^{t}\cdot\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(1)})\\\ M_{*}\nabla F(x^{(0)})\\\
\end{array}\right]+\sum_{i=0}^{t}T^{t-i}[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\nabla
F(y^{(i+1)})-p^{(i)}\right).$

Let $\Pi$ be the $2d\times 2d$ matrix with entries

$\Pi_{i,j}=\left\\{\begin{aligned} &1\qquad i\>\mathrm{odd},\;j=i\\\ &1\qquad
i\;\mathrm{even},\;j=2n+i\\\ &0\qquad\mathrm{otherwise}\end{aligned}\right.$

Then, by conjugation, we have

$\displaystyle\Pi\left[\begin{array}[]{cc}U,&0\\\
0,&U\end{array}\right]^{T}\left[\begin{array}[]{cc}(1+\theta)K,&-\theta K\\\
I,&0\\\ \end{array}\right]\left[\begin{array}[]{cc}U,&0\\\
0,&U\end{array}\right]\Pi^{T}$ $\displaystyle=$
$\displaystyle\Pi\left[\begin{array}[]{cc}(1+\theta)\Lambda,&-\theta\Lambda\\\
I,&0\end{array}\right]\Pi^{T}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}T_{1}&0&\cdots&0\\\ 0&T_{2}&\cdots&0\\\
\vdots&&\ddots&\vdots\\\ 0&0&\cdots&T_{d}\end{bmatrix}$

with

$\displaystyle
T_{i}=\begin{bmatrix}(1+\theta)\lambda_{i}&-\theta\lambda_{i}\\\
1&0\end{bmatrix}.$

That is $T$ is similar to the block diagonal matrix with $2\times 2$ diagonal
blocks $T_{i}$. And the eigenvalues of $T_{i}$ are the roots of

$\displaystyle a^{2}-a(1+\theta)\lambda_{i}+\theta\lambda_{i}=0.$

Because of

$\displaystyle(1-\pi)\nabla
F(x^{\star})\preceq\mathbb{E}\left[H^{-1}\right]\preceq\nabla F(x^{\star}),$

we have

$0\leq\lambda_{i}\leq\pi.$

Setting $\theta=\frac{1-\sqrt{1-\pi}}{1+\sqrt{1-\pi}}-\epsilon_{0}$, the
larger eigenvalue of $T_{1}$ is $q=1-\sqrt{1-\pi}+O(\sqrt{\epsilon_{0}})$. And
Let $S_{1}$ be the matrix consisting of eigenvectors of $T_{1}$. We denote

$c_{1}=\|S_{1}\|\cdot\|S_{1}^{-1}\|.$

In fact, $T_{1}$ dominates the convergence rate of $T^{t}$. Then we have

$\displaystyle\|T^{t}\|\leq c_{1}\cdot(1-\sqrt{1-\pi})^{t}=c_{1}q^{t}.$

Thus, we have

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(t+1)})\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]\right\rVert$
$\displaystyle\leq$ $\displaystyle
c_{1}q^{t}\cdot\left\lVert\left[\begin{array}[]{c}M_{*}\nabla F(x^{(1)})\\\
M_{*}\nabla F(x^{(0)})\\\
\end{array}\right]\right\rVert+c_{1}\sum_{i=0}^{t}q^{t-i}\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\nabla
F(y^{(i+1)})-p^{(i)}\right)\right\rVert$ $\displaystyle=$ $\displaystyle
c_{1}q^{t}\cdot
A_{0}+c_{1}\sum_{i=0}^{t}q^{t-i}\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\nabla
F(y^{(i+1)})-p^{(i)}\right)\right\rVert,$

where $A_{0}=\left\lVert\left[\begin{array}[]{c}M_{*}\nabla F(x^{(1)})\\\
M_{*}\nabla F(x^{(0)})\\\ \end{array}\right]\right\rVert$.

Due to Eqn. (3), (6) and using notation

$\tilde{\epsilon}_{1}=\frac{\epsilon_{1}}{12c_{1}\sqrt{\kappa}},$

we also have

$\displaystyle\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\nabla
F(y^{(i+1)})-p^{(i)}\right)\right\rVert$ $\displaystyle=$
$\displaystyle\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\left(\nabla
F(y^{(i+1)})-H^{(i+1)}p^{(i)}\right)\right)\right\rVert$ $\displaystyle\leq$
$\displaystyle\tilde{\epsilon}_{1}\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\right)\right\rVert\cdot\left\lVert\mathbb{E}\left(\nabla
F(y^{(i+1)})\right)\right\rVert$ $\displaystyle\leq$
$\displaystyle\tilde{\epsilon}_{1}\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\right)[\nabla^{2}F(x^{*})]^{\frac{1}{2}}M_{\star}\right\rVert\cdot\left\lVert
M_{\star}^{-1}\mathbb{E}\left(M_{\star}\nabla F(y^{(ii1)})\right)\right\rVert$
$\displaystyle\leq$
$\displaystyle\tilde{\epsilon}_{1}(1+\pi)\sqrt{\kappa}\left\lVert\mathbb{E}\left(M_{\star}\nabla
F(y^{(i+1)})\right)\right\rVert$ $\displaystyle=$
$\displaystyle\tilde{\epsilon}_{1}(1+\pi)\sqrt{\kappa}\left\lVert[1+\theta,-\theta]\mathbb{E}\left[\begin{array}[]{c}\nabla
F(x^{(i)})\\\ \nabla F(x^{(ti-1)})\end{array}\right]\right\rVert$
$\displaystyle\leq$ $\displaystyle
6\tilde{\epsilon}_{1}\sqrt{\kappa}\left\lVert\mathbb{E}\left[\begin{array}[]{c}\nabla
F(x^{(i)})\\\ \nabla F(x^{(i-1)})\end{array}\right]\right\rVert$

We use induction to prove our convergence rate, it is trivial for $t=1$
because

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(2)})\\\ M_{*}\nabla F(x^{(1)})\\\ \end{array}\right]\right\rVert\leq$
$\displaystyle c_{1}qA_{0}+\frac{\epsilon_{1}}{2c_{1}}A_{0}$ $\displaystyle=$
$\displaystyle c_{1}A_{0}\left(q+\frac{\epsilon_{1}}{2c_{1}^{2}}\right)$
$\displaystyle\leq$ $\displaystyle c_{1}A_{0}(q+\epsilon_{1})$
$\displaystyle<$ $\displaystyle 2c_{1}pA_{0}.$

The second inequality is because $c_{1}$ is no smaller than $1$. We assume
that

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}\nabla
F(x^{(t+1)})\\\ \nabla F(x^{(t)})\end{array}\right]\right\rVert\leq
Cp^{t},\;\mathrm{with}\;p=q+\epsilon_{1}\;\mathrm{and}\;C=2c_{1}A_{0}.$

Then, we have

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}M_{*}\nabla
F(x^{(t+1)})\\\ M_{*}\nabla F(x^{(t)})\\\ \end{array}\right]\right\rVert$
$\displaystyle\leq$ $\displaystyle c_{1}q^{t}\cdot
A_{0}+c_{1}\sum_{i=0}^{t}q^{t-i}\left\lVert[\nabla^{2}F(x^{*})]^{\frac{1}{2}}\mathbb{E}\left([H^{(i+1)}]^{-1}\nabla
F(y^{(i+1)})-p^{(i)}\right)\right\rVert$ $\displaystyle\leq$ $\displaystyle
c_{1}q^{t}\cdot
A_{0}+6\tilde{\epsilon}_{1}c_{1}\sqrt{\kappa}\sum_{i=1}^{t}q^{t-i}\left\lVert\mathbb{E}\left[\begin{array}[]{c}\nabla
F(x^{(i)})\\\ \nabla F(x^{(i-1)})\end{array}\right]\right\rVert$
$\displaystyle\leq$ $\displaystyle c_{1}q^{t}\cdot
A_{0}+6C\tilde{\epsilon}_{1}c_{1}\sqrt{\kappa}\sum_{i=1}^{t}q^{t-i}p^{i-1}$
$\displaystyle=$ $\displaystyle c_{1}q^{t}\cdot
A_{0}+6C\tilde{\epsilon}_{1}c_{1}\sqrt{\kappa}\sum_{i=1}^{t}q^{i}p^{i-i-1}$
$\displaystyle\leq$ $\displaystyle c_{1}q^{t}\cdot
A_{0}+6C\tilde{\epsilon}_{1}c_{1}\sqrt{\kappa}\frac{p^{t}}{p-q}$
$\displaystyle\leq$
$\displaystyle\left(c_{1}A_{0}\left(\frac{q}{p}\right)^{t}+\frac{C}{2}\right)\cdot
p^{t}$ $\displaystyle\leq$
$\displaystyle\left(c_{1}A_{0}+c_{1}A_{0}\right)\cdot p^{t}=Cp^{t}$

Therefore, for all $t=1,\dots,$, we have

$\displaystyle\left\lVert\mathbb{E}\left[\begin{array}[]{c}\nabla
F(x^{(t+1)})\\\ \nabla F(x^{(t)})\end{array}\right]\right\rVert\leq
Cp^{t}=2c_{1}p^{t}\left\lVert\left[\begin{array}[]{c}M_{*}\nabla F(x^{(1)})\\\
M_{*}\nabla F(x^{(0)})\\\ \end{array}\right]\right\rVert.$