# Incremental Quasi-Newton Methods with Faster Superlinear Convergence Rates

Zhuanghua Liu1, 2, Luo Luo3, Bryan Kian Hsiang Low1 The corresponding author

###### Abstract

We consider the finite-sum optimization problem, where each component function
is strongly convex and has Lipschitz continuous gradient and Hessian. The
recently proposed incremental quasi-Newton method is based on BFGS update and
achieves a local superlinear convergence rate that is dependent on the
condition number of the problem. This paper proposes a more efficient quasi-
Newton method by incorporating the symmetric rank-1 update into the
incremental framework, which results in the condition-number-free local
superlinear convergence rate. Furthermore, we can boost our method by applying
the block update on the Hessian approximation, which leads to an even faster
local convergence rate. The numerical experiments show the proposed methods
significantly outperform the baseline methods.

##  1 Introduction

We study the following finite-sum minimization problem:

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

where each individual function $f_{i}:{\mathbb{R}}^{d}\to{\mathbb{R}}$ is
strongly convex and has Lipschitz continuous gradient and Hessian. This
formulation is ubiquitous in various machine learning models, including
maximum likelihood estimation (MLE) (Bishop and Nasrabadi 2006; Bottou,
Curtis, and Nocedal 2018) and unsupervised learning problems (Murphy 2012;
Hastie et al. 2009). A notable example of the problem (1) is the empirical
risk minimization in supervised learning, where $n$ is the number of data
examples and $f_{i}(\cdot)$ corresponds to the loss function incurred by each
training instance.

In this paper, we are interested in solving the large-scale finite-sum
problem, that is, the number of components $n$ in formulation (1) is large. In
this scenario, accessing the exact gradient or Hessian over the entire dataset
is too expensive for each iteration. To circumvent this issue, stochastic or
incremental optimization methods were introduced since they only require
computing an estimation of the gradient or Hessian by a single sample (or a
small mini-batch of samples) at each round. The most popular of these methods
is stochastic gradient descent (SGD). It has been widely used in large-scale
optimization problems thanks to its cheap computational cost per iteration
(Bottou, Curtis, and Nocedal 2018). Applying the variance reduction (Schmidt,
Le Roux, and Bach 2017; Johnson and Zhang 2013; Defazio, Bach, and Lacoste-
Julien 2014; Zhang, Mahdavi, and Jin 2013) and acceleration techniques
(Nesterov 2003; Allen-Zhu 2017) can improve the vanilla SGD, and it achieves a
linear convergence rate with optimal incremental first-order oracle complexity
(Woodworth and Srebro 2016).

The second-order methods (Nesterov 2003) incorporate the additional curvature
information in every iteration, and it is possible to establish the local
superliner convergence rate with these methods. For the finite-sum problem
(1), Rodomanov and Kropotov (2016) proposed the Newton incremental method
(NIM), which requires accessing the exact gradient and exact Hessian of one
individual function in each iteration and attains the local superlinear
convergence rate. The classical quasi-Newton methods (Broyden, Dennis Jr, and
Moré 1973; Dennis and Moré 1974; Powell 1971) estimate the second-order
information with first-order oracle calls and still hold the superlinear
convergence rate. However, most of the stochastic variants for quasi-Newton
methods (Lucchi, McWilliams, and Hofmann 2015; Moritz, Nishihara, and Jordan
2016) that employ gradient estimators only achieve linear convergence rates.

Mokhtari, Eisen, and Ribeiro (2018) proposed the Incremental quasi-Newton
(IQN) method by using classical BFGS update (Broyden 1970; Fletcher 1970;
Goldfarb 1970; Shanno 1970), which is the first superlinear convergent quasi-
Newton method without exact second-order oracle call in each iteration.
However, the best-known analysis of IQN (Mokhtari, Eisen, and Ribeiro 2018)
only provided the asymptotic convergence result. Several follow-up works (Gao,
Koppel, and Ribeiro 2020; Lahoti et al. 2023) attempted to characterize the
convergence rate by fusing the greedy quasi-Newton update (Rodomanov and
Nesterov 2021a) into the framework of IQN. Specifically, Lahoti et al. (2023)
proposed sharpened lazy incremental quasi-Newton (SLIQN) by utilizing lazy
propagation strategy and showed it achieves the superlinear convergence rate
of ${\mathcal{O}}\big{(}(1-d^{-1}\varkappa^{-1})^{\left\lceil
t/n\right\rceil^{2}}\big{)}$, where $\varkappa$ is the condition number and
$t$ is the number of iterations. Gao, Koppel, and Ribeiro (2020) proposed the
Incremental Greedy BFGS (IGS) method with the same convergence rate as the
SLIQN, but it requires more expensive per-iteration complexity.

In this work, we propose an efficient quasi-Newton method named the Lazy
Incremental Symmetric Rank-1 (LISR-1) method for the finite-sum minimization
problem. Our approach takes advantage of the well-known symmetric rank-1 (SR1)
update to construct the Hessian estimator with sharper error bound than BFGS
methods, and it also exploits the lazy propagation strategy to maintain a low
per-iteration complexity. We show that LISR-1 achieves a local superlinear
convergence rate of ${\mathcal{O}}\big{(}(1-d^{-1})^{\left\lceil
t/n\right\rceil^{2}}\big{)}$, shaving off the dependency on the condition
number $\varkappa$ compared with the convergence rate achieved by SLIQN and
IGS. Each iteration of LISR-1 requires only $O(1)$ incremental
gradient/Hessian-vector oracle calls and ${\mathcal{O}}(d^{2})$ flops in
matrix operations, matching the existing IQN methods. Furthermore, we extend
LISR-1 by making use of the symmetric rank-$k$ update (Liu, Chen, and Luo
2023) to construct the more accurate Hessian estimator where $k<d$ is the rank
of the update, resulting in the block IQN method called Lazy Incremental
Symmetric Rank-$k$ (LISR-$k$). It enjoys the local convergence rate up to
${\mathcal{O}}\big{(}(1-kd^{-1})^{\left\lceil t/n\right\rceil^{2}}\big{)}$
with additional computational cost of ${\mathcal{O}}(kd^{2})$ flops per-
iteration. The numerical experiments on quadratic programming problems and the
model of regularized logistic regression demonstrate significant improvements
over baseline methods and confirm our theoretical findings.

#### Paper Organization

In Section 2, we provide a literature review for quasi-Newton methods and
their variants for finite-sum optimization problems. In Section 3, we
formalize the notations and assumptions of our problem and introduce the
background of the Broyden family update. In Section 4, we propose our LISR-1
method and provide its convergence analysis. In Section 5, we present the
LISR-$k$ method by incorporating the block-type update. In Section 6, we
demonstrate the numerical experiments to show the improved efficiency of the
proposed methods. Finally, we conclude this work in Section 7. All the proofs
and more experimental results are deferred to the appendix.

Algorithm | Computation Cost | Convergence Rate
---|---|---
IQN (Mokhtari, Eisen, and Ribeiro 2018) | ${\mathcal{O}}(d^{2})$ | asymptotic superlinear
IGS (Gao, Koppel, and Ribeiro 2020) | ${\mathcal{O}}(d^{3})$ | ${\mathcal{O}}\big{(}(1-d^{-1}\varkappa^{-1})^{\left\lceil t/n\right\rceil^{2}}\big{)}$
SLIQN (Lahoti et al. 2023) | ${\mathcal{O}}(d^{2})$ | ${\mathcal{O}}\big{(}(1-d^{-1}\varkappa^{-1})^{\left\lceil t/n\right\rceil^{2}}\big{)}$
LISR-1 (this work) | ${\mathcal{O}}(d^{2})$ | ${\mathcal{O}}\big{(}(1-d^{-1})^{\left\lceil t/n\right\rceil^{2}}\big{)}$
LISR-$k$ (this work) | ${\mathcal{O}}(kd^{2})$ | ${\mathcal{O}}\big{(}(1-kd^{-1})^{\left\lceil t/n\right\rceil^{2}}\big{)}$
Table 1: We compare the per-iteration computation cost and the convergence
rates of incremental fashion quasi-Newton methods. Note that the explicit
convergence rate of the vanilla IQN method still remains a mystery.

##  2 Related Work

In this section, we review related work of quasi-Newton methods and their
variants for large-scale optimization problems.

#### Classical Quasi-Newton Methods

Past decades have witnessed extensive research progress on quasi-Newton
methods. The main advantage of quasi-Newton methods is their capability to
reach a superlinear convergence without computing the exact Hessian or its
inverse. To estimate the second-order information, the classical quasi-Newton
methods are based on the secant equation and the corresponding closeness
criteria between successive Hessian estimations. The choice of closeness
criteria leads to different types of quasi-Newton methods, including Broyden’s
method (Broyden 1965; Broyden, Dennis Jr, and Moré 1973; Gay 1979), the
Davidon-Fletcher-Powell (DFP) method (Davidon 1991; Fletcher and Powell 1963),
the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (Broyden 1970; Fletcher
1970; Goldfarb 1970; Shanno 1970) and the symmetric rank-1 (SR1) method (Conn,
Gould, and Toint 1991). The asymptotic superlinear convergence of quasi-Newton
methods was established in the 1970s (Powell 1971; Dixon 1972a, b; Broyden,
Dennis Jr, and Moré 1973; Dennis and Moré 1974), while the explicit
superlinear rates of quasi-Newton methods were obtained only recently.
Rodomanov and Nesterov (2021a) first proposed greedy quasi-Newton methods and
gave its non-asymptotic superlinear convergence guarantees. Later, Lin, Ye,
and Zhang (2022) provided a sharper analysis for these methods. After that,
Rodomanov and Nesterov (2021b, c); Jin and Mokhtari (2023); Ye et al. (2021)
established the explicit rates for the classical (secant equation-based)
quasi-Newton methods.

#### Block Quasi-Newton Methods

Schnabel (1983) proposed block quasi-Newton methods. These methods construct
the Hessian estimator along multiple directions during each iteration, and
they achieve better empirical performance than classical quasi-Newton methods
like BFGS (O’Leary and Yeremin 1994). After several decades, the superlinear
convergence of these methods was established by Gao and Goldfarb (2018);
Gower, Goldfarb, and Richtárik (2016); Gower and Richtárik (2017). Very
recently, Liu, Chen, and Luo (2023) presented explicit superlinear convergence
rates of block quasi-Newton methods, which explains why the use of multiple
directions benefits the convergence behaviors.

#### Stochastic/Incremental Quasi-Newton Methods

Due to the sheer volume of data in modern machine learning applications,
researchers have been investigating the extension of quasi-Newton methods on
large-scale optimization problems. Several early works established the
stochastic quasi-Newton methods to reduce the computational cost at each
iteration (Byrd et al. 2016; Mokhtari and Ribeiro 2014, 2015; Moritz,
Nishihara, and Jordan 2016; Lucchi, McWilliams, and Hofmann 2015; Chang, Sun,
and Zhang 2019), but these methods cannot obtain the superlinear convergences
like classical quasi-Newton methods. Incremental quasi-Newton methods (IQN)
(Mokhtari, Eisen, and Ribeiro 2018; Gao, Koppel, and Ribeiro 2020; Lahoti et
al. 2023) use the aggregated information to construct a more accurate gradient
and Hessian estimator, which leads to superlinear convergence. We compare the
proposed methods with related work in Table 1.

##  3 Preliminaries

In this section, we formalize the notations and assumptions throughout this
paper, then we introduce the well-known Broyden family updates which are
widely used in quasi-Newton methods.

###  3.1 Notations

We denote $e_{i}\in{\mathbb{R}}^{d}$ as the $i$-th standard basis vector of
$d$-dimensional Euclidean space, where $i\in[d]$. We define the index $i_{t}$
as $t~{}{\rm mod}~{}n$. For vectors $u,v\in{\mathbb{R}}^{d}$, we denote their
inner product by $\langle u,v\rangle\coloneqq u^{\top}v$. We use
$\left\lVert\cdot\right\rVert$ to represent the Euclidean norm of the vector
and the spectral norm of the matrix. Given a positive semi-definite matrix
$A\in{\mathbb{R}}^{d\times d}$ and a vector $u\in{\mathbb{R}}^{d}$, we define
the norm of $u$ with respect to $A$ as $\left\lVert
u\right\rVert_{A}\coloneqq\sqrt{\langle u,Au\rangle}$. We let

$E_{k}(A)=[e_{i_{1}};\ldots;e_{i_{k}}]\in{\mathbb{R}}^{d\times k},$ (2)

where $i_{1},\ldots,i_{k}$ are the indices for the largest $k$ entries in the
diagonal of $A$. We also use ${\rm tr}(\cdot)$ to present the trace of a
square matrix. Additionally, we denote the solution of problem (1) as
$x^{*}\coloneqq\operatorname*{arg\,min}_{x\in{\mathbb{R}}^{d}}f(x)$.

###  3.2 Assumptions

In the remainder of this paper, we always suppose Problem (1) satisfies the
following assumptions.

######  Assumption 3.1.

We suppose each function $f_{i}(\cdot)$ is twice-differentiable, $L$-smooth
and $\mu$-strongly convex, i.e., there exist constants $L>0$ and $\mu>0$ such
that

$\mu I\preceq\nabla^{2}f_{i}(x)\preceq LI$ (3)

for any $x\in{\mathbb{R}}^{d}$.

######  Assumption 3.2.

We suppose each $f_{i}(\cdot)$ has a $\tilde{L}$-Lipschitz continuous Hessian,
i.e., there exists a constant $\tilde{L}$ such that

$\left\lVert\nabla^{2}f_{i}(x)-\nabla^{2}f_{i}(y)\right\rVert\leq\tilde{L}\left\lVert
x-y\right\rVert.$

for any $x,y\in{\mathbb{R}}^{d}$.

The strong convexity and the Lipschitz continuity of Hessian in our
assumptions imply that each $f_{i}(\cdot)$ is strongly self-concordant with
constant $M\coloneqq\tilde{L}\mu^{-3/2}$ (Rodomanov and Nesterov 2021a), i.e,
we have

$\displaystyle\nabla^{2}f_{i}(y)-\nabla^{2}f_{i}(x)\preceq M\left\lVert
y-x\right\rVert_{\nabla^{2}f_{i}(z)}\nabla^{2}f_{i}(w)$

for any $x,y,z,w\in{\mathbb{R}}^{d}$.

Additionally, we let $\varkappa\coloneqq L/\mu$ be the condition number of our
problem which could be very large in practice.

###  3.3 Broyden Family Update

Many popular quasi-Newton methods such as DFP, BFGS, and SR1 belong to the
Broyden family update (Nocedal and Wright 1999, Section 6.3), which is defined
as follows.

######  Definition 3.3.

Let $G\in{\mathbb{R}}^{d\times d}$ and $A\in{\mathbb{R}}^{d\times d}$ be two
positive define matrices satisfying $G\succeq A$. For any non-zero
$u\in{\mathbb{R}}^{d}$ and $\tau\in[0,1]$, if $Gu=Au$, we define ${\rm
Broyd}_{\tau}(G,A,u)\coloneqq G$. Otherwise, we define

$\begin{split}\\!&\\!{\rm Broyd}_{\tau}(G,A,u)\\\
\\!\\!\coloneqq&\tau\left[G-\frac{Auu^{\top}G+Guu^{\top}A}{u^{\top}Au}+\left(\frac{u^{\top}Gu}{u^{\top}Au}+1\right)\frac{Auu^{\top}A}{u^{\top}Au}\right]\\\
&\quad+(1-\tau)\left[G-\frac{(G-A)uu^{\top}(G-A)}{u^{\top}(G-A)u}\right].\end{split}$
(4)

We can recover several well-known quasi-Newton methods by taking the different
values of $\tau$:

  * •

For $\tau=1$, Eq. (4) corresponds to the DFP update

$\begin{split}{\rm
DFP}(G,A,u)\coloneqq&G-\frac{Auu^{\top}G+Guu^{\top}A}{u^{\top}Au}\\\
&\quad+\left(\frac{u^{\top}Gu}{u^{\top}Au}+1\right)\frac{Auu^{\top}A}{u^{\top}Au}.\end{split}$

  * •

For $\tau=\dfrac{u^{\top}Au}{u^{\top}Gu}$, we recover the BFGS update

${\rm BFGS}(G,A,u)\coloneqq
G-\frac{Guu^{\top}G}{u^{\top}Gu}+\frac{Auu^{\top}A}{u^{\top}Au}.$

  * •

For $\tau=0$, we achieve the SR1 update

${\rm SR1}(G,A,u)\coloneqq G-\frac{(G-A)uu^{\top}(G-A)}{u^{\top}(G-A)u}.$ (5)

We can generalize the Broyden family updates with multiple directions (Gao and
Goldfarb 2018; Liu, Chen, and Luo 2023; Gower, Goldfarb, and Richtárik 2016;
Gower and Richtárik 2017). In particular, Liu, Chen, and Luo (2023) establish
the block version of the SR1 update called the symmetric rank-$k$ (SR-$k$)
update, which is defined as follows.

######  Definition 3.4.

Let $A\in{\mathbb{R}}^{d\times d}$ and $G\in{\mathbb{R}}^{d\times d}$ be two
positive-definite matrices satisfying $G\succeq A$. For any full rank matrix
$U\in{\mathbb{R}}^{d\times k}$ with $k<d$, we define
$\textrm{SR-}k(G,A,U)\coloneqq G$ if $GU=AU$. Otherwise, we define

$\begin{split}\textrm{SR-}k(G,A,U)\coloneqq
G-(G\\!-\\!A)U(U^{\top}(G\\!-\\!A)U)^{{\dagger}}U^{\top}(G\\!-\\!A).\end{split}$

######  Remark 3.5.

Note that the SR-$k$ update shown in the above definition is equivalent to the
SR1 update when $k=1$.

##  4 Methodology

In this section, we propose the lazy incremental symmetric rank*1 (LISR-1)
method and provide theoretical analysis to show it enjoys condition number-
free local superlinear convergence.

###  4.1 The Algorithm

Algorithm 1 LISR-1

1: Input: $x^{0}\in{\mathbb{R}}^{d}$ and $\\{B_{i}^{0}\in{\mathbb{R}}^{d\times
d}\\}_{i=1}^{n}$.

2: Initialize $t=0$ and $z_{i}^{0}=x^{0}$ for any $i\in[n]$.

3: while not converged do

4: Update $x^{t+1}$ as per (6).

5: Update $z_{i}^{t+1}$ as per (7).

6: Update $B_{i}^{t+1}$ as per (8)-(11).

7: Increment the iteration counter $t$.

8: end while

9: Output: $x^{t}$.

We first introduce the main intuitions of LISR-1. For each component function
$f_{i}(x)$, we consider its quadratic approximation at point
$z_{i}^{t}\in{\mathbb{R}}^{d}$ as

$\displaystyle\begin{split}f_{i}(x)\approx&{\tilde{f}}_{i}^{t}(x)\\\
=&f_{i}(z_{i}^{t})+\nabla
f_{i}(z_{i}^{t})^{\top}(x-z_{i}^{t})+\frac{1}{2}(x-z_{i}^{t})^{\top}B_{i}^{t}(x-z_{i}^{t}),\end{split}$

where we estimate $\nabla^{2}f_{i}(z_{i}^{t})$ by a positive-definite matrix
$B_{i}^{t}\in{\mathbb{R}}^{d\times d}$. Then we obtain $x^{t+1}$ by minimizing
the average of $\\{\tilde{f}_{i}^{t}(x)\\}_{i=1}^{n}$, which has the closed
form solution

$\displaystyle\begin{split}x^{t+1}=&\operatorname*{arg\,min}_{x\in{\mathbb{R}}^{d}}\frac{1}{n}\sum_{i=1}^{n}\tilde{f}_{i}^{t}(x)\\\
=&\left(\sum_{i=1}^{n}B_{i}^{t}\right)^{-1}\left(\sum_{i=1}^{n}B_{i}^{t}z_{i}^{t}-\sum_{i=1}^{n}\nabla
f_{i}(z_{i}^{t})\right).\end{split}$ (6)

We only update one of $\\{z_{i}^{t}\\}_{i=1}^{n}$ at each iteration in a
cyclic fashion to make the algorithm efficient, that is

$\displaystyle
z_{i}^{t+1}=\begin{cases}x^{t+1},&\text{if}~{}i=i_{t},\\\\[2.84544pt]
z_{i}^{t},&\text{otherwise},\end{cases}$ (7)

where $i_{t}=t~{}{\rm mod}~{}n$ is the index of the component we choose at the
$t$-th iteration.

We also wish to construct the Hessian estimators efficiently and keep a fast
convergence rate. In particular, we introduce the scaling parameter
$\omega^{t+1}$ and apply the SR1 update on one of the individual Hessian
estimators in each iteration:

  * •

For $i=i_{t}$, we let

$\begin{split}\\!\\!\\!B_{i}^{t+1}=\omega^{t+1}\textrm{SR1}(B_{i}^{t},\nabla^{2}f_{i}(z_{i}^{t+1}),\bar{u}(B_{i}^{t},\nabla^{2}f_{i}(z_{i}^{t+1}))),\end{split}$
(8)

where $\bar{u}(\cdot,\cdot)$ is the greedy direction which is defined as

$\displaystyle\bar{u}(G,A)\coloneqq\operatorname*{arg\,max}_{u\in\\{e_{i}\\}_{i=1}^{d}}u^{\top}(G-A)u.$
(9)

  * •

For $i\neq i_{t}$, we let

$\begin{split}B_{i}^{t+1}=\omega^{t+1}B_{i}^{t}.\end{split}$ (10)

Additionally, we set

$\displaystyle\omega^{t}=\begin{cases}\big{(}1+M\sqrt{L}r_{0}\rho^{\left\lceil\frac{t}{n}\right\rceil}\big{)}^{2},&\text{if}~{}n~{}{\rm
mod}~{}t=0,\\\\[2.84544pt] 1,&\text{otherwise},\end{cases}$ (11)

for some $\rho\in(0,1-d^{-1})$ and let $r_{0}$ be an upper bound of
$\left\lVert x_{0}-x^{*}\right\rVert$. This setting implies the step of
scaling is executed once every $n$ iterations.

We present the whole procedure of the proposed LISR-1 in Algorithm 1. We can
verify that the per-iteration cost of our algorithm is ${\mathcal{O}}(d^{2})$
flops. Notice that the main cost of LISR-1 comes from the computation of Eq.
(6), which is dominated by maintaining the inverse of the following sum of
individual Hessian estimators

$\bar{B}^{t+1}\coloneqq\sum_{i=1}^{n}B_{i}^{t+1}.$

We can rewrite the above matrix in the recursive form as

$\bar{B}^{t+1}=\bar{B}^{t}+B_{i_{t}}^{t+1}-B_{i_{t}}^{t}.$ (12)

In the case of $t~{}{\rm mod}~{}n\neq 0$, no scaling is performed since we
have $\omega^{t+1}=1$. Denote $\bar{u}^{t}$ as the abbreviation of
$\bar{u}(B_{i_{t}}^{t},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t}))$, then applying the
Sherman-Morrison formula on Eq. (12) implies

$(\bar{B}^{t+1})^{-1}=(\bar{B}^{t})^{-1}+\frac{(\bar{B}^{t})^{-1}v^{t}(v^{t})^{\top}(\bar{B}^{t})^{-1}}{(\bar{u}^{t})^{\top}v^{t}-(v^{t})^{\top}(\bar{B}^{t})^{-1}v^{t}},$
(13)

where $v^{t}$ is defined as

$v^{t}=(B_{i_{t}}^{t}-\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}))\bar{u}^{t}.$

It is easy to observe that computing the right-hand side of Eq. (13) takes
${\mathcal{O}}(d^{2})$ flops for given $(\bar{B}^{t})^{-1}$ and $v^{t}$. In
the case of $t~{}{\rm mod}~{}n=0$, each Hessian estimator may be scaled by a
factor $\omega^{t}\neq 1$, which results in the additional computational cost
of ${\mathcal{O}}(nd^{2})$ flops. However, the amortized per-iteration
complexity of this step is still ${\mathcal{O}}(d^{2})$ because the scaling
occurs once per $n$ iterations. We provide a more efficient implementation of
LISR-1 in the appendix.

###  4.2 Convergence Analysis

We analyze the convergence of LISR-1 by considering the Euclidean distance to
the optimal solution $x^{*}$. Firstly, the formula (6) indicates the general
result:

######  Lemma 4.1.

The iteration formula (6) satisfies

$\begin{split}\left\lVert
x^{t+1}-x^{*}\right\rVert\leq&\frac{\tilde{L}\Gamma^{t}}{2}\sum_{i=1}^{n}\left\lVert
z_{i}^{t}-x^{*}\right\rVert^{2}\\\ &+\Gamma^{t}\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert,\end{split}$ (14)

for all $t\geq 1$, where
$\Gamma^{t}\coloneqq\big{\|}\big{(}\sum_{i=1}^{n}B_{i}^{t}\big{)}^{-1}\big{\|}$.

######  Remark 4.2.

Notice that the proof of Lemma 4.1 only requires the Lipschitz continuity of
each $\nabla^{2}f_{i}(\cdot)$ and the iteration formula (6). The validity of
this lemma does not rely on the specific choice of Hessian estimators
$\\{B_{i}^{t}\\}_{i=1}^{n}$ and it also can be used to analyze the other
incremental fashion methods (Mokhtari, Eisen, and Ribeiro 2018; Lahoti et al.
2023; Gao, Koppel, and Ribeiro 2020).

In view of Lemma 4.1, the more accurate Hessian estimator
$B_{i}^{t}\approx\nabla^{2}f_{i}(z_{i}^{t})$ can lead to the tighter upper
bound of $\left\lVert x^{t+1}-x^{*}\right\rVert$. Hence, the key to showing
the advantage of the proposed method is bounding the difference between
$B_{i}^{t}$ and $\nabla^{2}f_{i}(z_{i}^{t})$. In particular, we introduce the
quantity

$\nu(G,A)\coloneqq\frac{d\varkappa{\rm tr}(G-A)}{{\rm tr}(A)},$ (15)

to describe the difference between two positive definite matrices
$G\in{\mathbb{R}}^{d\times d}$ and $A\in{\mathbb{R}}^{d\times d}$ such that
$G\succeq A$. Based on the measure $\nu(\cdot,\cdot)$ and Lemma 4.1, we
provide the linear convergence of the distance to solution and the error of
Hessian approximation as follows.

######  Lemma 4.3.

For any $\rho$ satisfying $\rho\in(0,1-d^{-1})$, there exist positive
constants $r_{0}$ and $\sigma_{0}$ such that running LISR*1 (Algorithm 1) with
the initial conditions $\left\lVert x^{0}-x^{*}\right\rVert\leq r_{0}$,
$B_{i}^{0}\succeq\omega^{0}\nabla^{2}f_{i}(z_{i}^{0})$ and
$\nu((\omega^{0})^{-1}{B_{i}^{0}},\nabla^{2}f_{i}(x^{0}))\leq\sigma_{0}$ for
any $i=1,\dots,n$ results in

$\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\rho^{\left\lceil\frac{t+1}{n}\right\rceil}\left\lVert
x^{0}-x^{*}\right\rVert$ (16)

and

$\begin{split}\nu\big{(}(\omega^{t+1})^{-1}B_{i_{t}}^{t+1},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1})\big{)}\leq\left(1-\frac{1}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\delta,\end{split}$
(17)

where $M=\tilde{L}/\mu^{{3}/{2}}$, $\omega_{t}$ follows the definition in (11)
and

$\displaystyle\delta\coloneqq\left(\sigma_{0}+\frac{4MdL^{{3}/{2}}\mu^{-1}r_{0}}{1-(1-d^{-1})^{-1}\rho}\right)\exp\left(\frac{4M\sqrt{L}r_{0}}{1-\rho}\right)$

######  Remark 4.4.

In the proof of Lemma 4.3, we show that the relation
$(\omega^{t+1})^{-1}B_{i_{t}}^{t+1}\succeq\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1})$
holds for each iteration. This guarantees the update rule (6) and the error
measure
$\nu\big{(}(\omega^{t+1})^{-1}B_{i_{t}}^{t+1},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1})\big{)}$
are well-defined.

We establish the mean-superlinear convergence based on Lemma 4.3.
Specifically, we have the following result.

######  Lemma 4.5.

Following the initial conditions of Lemma 4.3, the sequence of iterates
generated by the LISR*1 method (Algorithm 1) satisfies

$\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\left(1-\frac{1}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\cdot\frac{1}{n}\sum_{i=1}^{n}\left\lVert
x^{t+1-i}-x^{*}\right\rVert.$

Using Lemma 4.5, we can achieve the local superlinear convergence rate of the
proposed LISR*1 method by induction. We formally present our main result as
follows.

######  Theorem 4.6.

For the sequence $\\{x_{t}\\}$ generated by LISR*1 (Algorithm 1) with the
initial conditions shown in Lemma 4.3, there exists a sequence
$\\{\zeta^{l}\\}$ such that $\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\zeta^{\left\lfloor{t}/{n}\right\rfloor}$ for any
$t\geq 1$ and it satisfies

$\zeta^{l}\leq r_{0}\left(1-\frac{1}{d}\right)^{\frac{(l+2)(l+1)}{2}}.$ (18)

###  4.3 Discussion

The convergence analysis in the last subsection shows that LISR*1 enjoys the
condition number-free superlinear convergence rate, which is significantly
better than all of the existing incremental fashion quasi-Newton methods (see
Table 1). The improvement is due to that we adopt the greedy SR1 update to
maintain the Hessian estimator in formula (8) and the analysis characterizes
the Hessian approximation error by the measure $\nu(\cdot,\cdot)$ defined in
(15). In contrast, the prior methods IGS (Gao, Koppel, and Ribeiro 2020) and
SLIQN (Lahoti et al. 2023) only consider the general Broyden family update and
characterize the Hessian approximation error by the measure $\sigma(G,A)={\rm
tr}(A^{-1}(G-A))$ for positive definite $G\in{\mathbb{R}}^{d\times d}$ and
$A\in{\mathbb{R}}^{d\times d}$, which leads to additional dependency on
condition number in the superlinear convergence rate.111These work present
their theoretical results by analyzing BFGS update, while their analysis can
be directly applied to the general Broyden family update and achieves the
identical convergence rate. On the other hand, the implementations of these
methods are more complicated than ours. Concretely, IGS requires scaling a
Hessian estimator at each iteration which results in ${\mathcal{O}}(d^{3})$
computational cost, and SLIQN maintains $B_{i_{t}}^{t+1}$ by a combination of
secant equation-based and greedy Broyden family updates while our LISR-1 only
has one step of greedy SR1 update (8).

|  |
---|---|---
(a) $\xi=4$, $\varkappa=3.03\times 10^{2}$ | (b) $\xi=8$, $\varkappa=3.12\times 10^{4}$ | (c) $\xi=12$, $\varkappa=3.12\times 10^{6}$
Figure 1: Normalized error vs. the number of effective passes for the
quadratic programming problem.

##  5 Extension to Block Quasi-Newton Methods

It is possible to incorporate the idea of block quasi-Newton methods into the
framework of the LISR-1. Specifically, we only need to modify Line 6 of
Algorithm 1 by replacing the update rule (8) with

$\begin{split}\\!\\!\\!\\!B_{i_{t}}^{t+1}\\!=\\!\omega^{t+1}\textrm{SR-}k(B_{i_{t}}^{t},\\!\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}),\\!\bar{U}(B_{i_{t}}^{t},\\!\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}))),\end{split}$
(19)

where $\bar{U}(\cdot,\cdot)$ contains greedy directions which is defined as

$\bar{U}(G,A)=E_{k}(G-A).$ (20)

We name the variant of LISR-1 with the above modification as Lazy Incremental
Symmetric Rank-$k$ (LISR-$k$) method.

The LISR-$k$ method requires ${\mathcal{O}}(kd^{2})$ flops in each iteration.
Since we typically set $k$ to be much smaller than $d$, such computational
cost is acceptable. Similar to the previous analysis, the cost of LISR-$k$ is
dominated by maintaining the inverse of the sum of individual Hessian
estimators

$\bar{B}^{t+1}\coloneqq\sum_{i=1}^{n}B_{i}^{t+1},$

which can be written as $\bar{B}^{t}+B_{i_{t}}^{t+1}-B_{i_{t}}^{t}$. The main
difference between the two algorithms is the update on $\bar{B}^{t+1}$ (its
inverse) in the case of $t~{}{\rm mod}~{}n\neq 0$. For the LISR-$k$ method, we
have

$\begin{split}\bar{B}^{t+1}=&\bar{B}^{t}-V^{t}\left((\bar{U}^{t})^{\top}V^{t}\right)^{-1}(V^{t})^{\top},\end{split}$

where we define
$V^{t}=\left(B_{i_{t}}^{t}-\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1})\right)\bar{U}^{t}\in{\mathbb{R}}^{d\times
k}$ and
$\bar{U}^{t}=\bar{U}(B_{i_{t}}^{t},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}))\in{\mathbb{R}}^{d\times
k}$. Applying the Sherman-Morrison formula, we achieve

$\begin{split}(\bar{B}^{t+1})^{-1}=&(\bar{B}^{t})^{-1}+(\bar{B}^{t})^{-1}V^{t}(D^{t})^{-1}(V^{t})^{\top}(\bar{B}^{t})^{-1},\end{split}$
(21)

where
$D^{t}=(\bar{U}^{t})^{\top}V^{t}-(V^{t})^{\top}(\bar{B}^{t})^{-1}V^{t}\in{\mathbb{R}}^{k\times
k}$. It can be observed that constructing $D^{t}$ takes
${\mathcal{O}}(kd^{2})$ flops for given $(\bar{B}_{t})^{-1}$. Additionally,
the complexity of computing $(D^{t})^{-1}$ is not the leading cost since we
take $k\ll d$. Hence, the total cost for computing Eq. (21) is
${\mathcal{O}}(kd^{2})$ flops. Similar to LISR-1, the setting of
$\omega^{t+1}$ guarantees the scaling occurs once every $n$ iterations and its
amortized per-iteration complexity is no more than ${\mathcal{O}}(kd^{2})$
flops.

#### Even Faster Convergence Rate

The rank-$k$ update in the LISR-$k$ leads to sharper upper bounds on the
distance to optimal solution and approximation error of Hessian estimators.
Compared with Lemma 4.5, the LISR-$k$ holds the following tighter upper bounds

$\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\left(1-\frac{k}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\frac{1}{n}\sum_{i=1}^{n}\left\lVert
x^{t+1-i}-x^{*}\right\rVert$

and

$\nu((\omega^{t+1})^{-1}B_{i_{t}}^{t+1},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}))\leq\left(1-\frac{k}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\delta.$

Consequently, we can show the mean-superlinear convergence result like Lemma
4.5, and the new result improves the base of convergence rate from $1-d^{-1}$
to $1-kd^{-1}$. Finally, we achieve the main result of the LISR-$k$ method as
follows.

######  Theorem 5.1.

We follow the initial conditions of Lemma 4.3 but initialize $\rho$ with
$\rho\in(0,1-kd^{-1})$. For the sequence of iterates $\\{x^{t}\\}$ generated
by the LISR-$k$ method, there exists sequence $\\{\zeta^{l}\\}$ such that
$\left\lVert
x^{t}-x^{*}\right\rVert\leq\zeta^{\left\lfloor(t-1)/n\right\rfloor}$ for any
$t\geq 1$ and it satisfies

$\zeta^{l}\leq r_{0}\left(1-\frac{k}{d}\right)^{\frac{(l+2)(l+1)}{2}}.$ (22)

######  Remark 5.2.

For $k=1$, the LISR-$k$ method degenerates to the LISR-1 method. For $k\geq
2$, the superlinear convergence rate of LISR-$k$ (Theorem 22) is strictly
tighter than the counterpart of LISR-$1$ (Theorem 18).

##  6 Experiments

We compare the proposed methods LISR-1 and LISR-$k$ with baseline methods
including IQN (Mokhtari, Eisen, and Ribeiro 2018) and SLIQN (Lahoti et al.
2023). We test all methods on the problems of quadratic programming and
regularized logistic regression. For the LISR-$k$ method, we set $k=5$ for all
of the cases. For the fairness of comparison, we run all algorithms from the
same initial point.

###  6.1 Quadratic Function Minimization

We consider the following quadratic function minimization problem

$\min_{x\in{\mathbb{R}}^{d}}f(x)\coloneqq\frac{1}{n}\sum_{i=1}^{n}\Big{(}\frac{1}{2}\langle
x,A_{i}x\rangle+\langle b_{i},x\rangle\Big{)},$ (23)

where $A_{i}\in{\mathbb{R}}^{d\times d}$ is positive definite and
$b_{i}\in{\mathbb{R}}^{d}$. Following the setup of Mokhtari, Eisen, and
Ribeiro (2018), we let each $A_{i}$ be diagonal matrix by setting the first
half of diagonal entries be independent uniformly sampled from
$[1,10^{\xi/2}]$ while the others are independent uniformly sampled from
$[10^{-\xi/2},1]$, where $\xi>0$ is the parameter that affects the condition
number of the problem. For each $b_{i}$, we let its entries be independently
uniformly sampled from $[0,10^{3}]$.

We run the experiments by taking $n=1000$, $d=50$ and $\xi\in\\{4,8,12\\}$,
and we present the results in Figure 1. We observe that the condition number
heavily affects the convergence behaviors of IQN and SLIQN, while the proposed
methods LISR-1 and LISR-$k$ are insensitive to the varying condition numbers.
These results validate our theoretical analysis since we have shown the
superlinear convergence rates of our methods do not depend on the condition
number.

|  |
---|---|---
(a) a9a | (b) w8a | (c) ijcnn1
|  |
(e) mushrooms | (f) phishing | (g) svmguide3
|  |
(h) german.numer | (i) splice | (j) covtype
Figure 2: Normalized error vs. the number of effective passes for the
regularized logistic regression problem on several real-world datasets .

###  6.2 Regularized Logistic Regression

We consider $\ell_{2}$-regularized logistic regression problem

$\begin{split}\min_{x\in{\mathbb{R}}^{d}}f(x)\\!\coloneqq\\!\frac{1}{n}\sum_{i=1}^{n}\ln(1\\!+\\!\exp(-y_{i}\langle
x,z_{i}\rangle))\\!+\\!\frac{\lambda}{2}\\!\left\lVert
x\right\rVert^{2}\\!,\end{split}$ (24)

where $z_{i}\in{\mathbb{R}}^{d}$ is the feature of the $i$-th training sample
and $y_{i}\in\\{1,-1\\}$ is the corresponding labels. We conduct our
experiments on nine real-world datasets (“a9a”, “w8a”, “ijcnn”, “mushrooms”,
“phishing”, “svmguide3”, “german.numer”, “splice” and ‘covtype”) from LIBSVM
repository. We take $\lambda=10^{-3}$ for “a9a”, “mushrooms”, “svmguide3”,
“german.numer”, “covtype” and $\lambda=10^{-4}$ for others.

We present the experimental results in Figure 2. We observe that the proposed
LISR-$k$ significantly outperforms other methods on all datasets. The LISR-1
enjoys a faster convergence rate than IQN and SLIQN when it starts to
converge, while it may be slower at the early stage. We conjecture that IQN
and SLIQN contain the steps of classical quasi-Newton updates. By accessing
the exact gradient information, Rodomanov and Nesterov (2021b, c)
theoretically showed that classical quasi-Newton methods converge faster than
greedy quasi-Newton methods at the early stage. We empirically observe similar
results for incremental quasi-Newton methods, while the rigorous theory for
such a phenomenon is still unclear. On the other hand, the block update in
LISR-$k$ leads to much better Hessian estimators. Hence, the early stage of
LISR-$k$ only contains a few iterations.

##  7 Conclusion

This paper has proposed the efficient incremental quasi-Newton method called
LISR-1 and its extension named LISR-$k$ method for the finite-sum convex
optimization. We have theoretically shown the proposed methods enjoy faster
superlinear convergence rates than the state-of-the-art incremental quasi-
Newton methods. The numerical experiments on quadratic programming and
regularized logistic regression also validate the advantages of the proposed
methods over existing IQN baselines.

In future work, it is interesting to study incremental quasi-Newton methods
for more general settings, such as minimizing nonconvex functions (Wang et al.
2017; Yang et al. 2021). It is also possible to leverage the idea to design
efficient incremental quasi-Newton methods for solving minimax problems (Liu
and Luo 2022; Liu et al. 2022) or nonlinear equations (Liu et al. 2023).

## Acknowledgments

This research/project is supported by the National Research Foundation,
Singapore under its AI Singapore Programme (AISG Award No:
AISG2-PhD-2023-08-043T-J). This research is part of the programme DesCartes
and is supported by the National Research Foundation, Prime Minister’s Office,
Singapore under its Campus for Research Excellence and Technological
Enterprise (CREATE) programme. Luo Luo is supported by National Natural
Science Foundation of China (No. 62206058) and Shanghai Sailing Program
(22YF1402900).

## References

  * Allen-Zhu (2017) Allen-Zhu, Z. 2017.  Katyusha: The first direct acceleration of stochastic gradient methods.  In _Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing_ , 1200–1205.
  * Bishop and Nasrabadi (2006) Bishop, C. M.; and Nasrabadi, N. M. 2006.  _Pattern recognition and machine learning_ , volume 4.  Springer.
  * Bottou, Curtis, and Nocedal (2018) Bottou, L.; Curtis, F. E.; and Nocedal, J. 2018.  Optimization methods for large-scale machine learning.  _SIAM review_ , 60(2): 223–311.
  * Broyden (1965) Broyden, C. G. 1965.  A class of methods for solving nonlinear simultaneous equations.  _Mathematics of computation_ , 19(92): 577–593.
  * Broyden (1970) Broyden, C. G. 1970.  The convergence of single-rank quasi-Newton methods.  _Mathematics of Computation_ , 24(110): 365–382.
  * Broyden, Dennis Jr, and Moré (1973) Broyden, C. G.; Dennis Jr, J. E.; and Moré, J. J. 1973.  On the local and superlinear convergence of quasi-Newton methods.  _IMA Journal of Applied Mathematics_ , 12(3): 223–245.
  * Byrd et al. (2016) Byrd, R. H.; Hansen, S. L.; Nocedal, J.; and Singer, Y. 2016.  A stochastic quasi-Newton method for large-scale optimization.  _SIAM Journal on Optimization_ , 26(2): 1008–1031.
  * Chang, Sun, and Zhang (2019) Chang, D.; Sun, S.; and Zhang, C. 2019.  An accelerated linearly convergent stochastic L-BFGS algorithm.  _IEEE Transactions on neural networks and learning systems_ , 30(11): 3338–3346.
  * Conn, Gould, and Toint (1991) Conn, A. R.; Gould, N. I.; and Toint, P. L. 1991.  Convergence of quasi-Newton matrices generated by the symmetric rank one update.  _Mathematical programming_ , 50(1-3): 177–195.
  * Davidon (1991) Davidon, W. C. 1991.  Variable metric method for minimization.  _SIAM Journal on optimization_ , 1(1): 1–17.
  * Defazio, Bach, and Lacoste-Julien (2014) Defazio, A.; Bach, F.; and Lacoste-Julien, S. 2014.  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _Twenty-seventh Conference on Neural Information Processing Systems_.
  * Dennis and Moré (1974) Dennis, J. E.; and Moré, J. J. 1974.  A characterization of superlinear convergence and its application to quasi-Newton methods.  _Mathematics of computation_ , 28(126): 549–560.
  * Dixon (1972a) Dixon, L. 1972a.  Quasi-Newton algorithms generate identical points.  _Mathematical Programming_ , 2: 383–387.
  * Dixon (1972b) Dixon, L. 1972b.  Quasi Newton techniques generate identical points II: the proofs of four new theorems.  _Mathematical Programming_ , 3: 345–358.
  * Fletcher (1970) Fletcher, R. 1970.  A new approach to variable metric algorithms.  _The computer journal_ , 13(3): 317–322.
  * Fletcher and Powell (1963) Fletcher, R.; and Powell, M. J. 1963.  A rapidly convergent descent method for minimization.  _The computer journal_ , 6(2): 163–168.
  * Gao and Goldfarb (2018) Gao, W.; and Goldfarb, D. 2018.  Block BFGS methods.  _SIAM Journal on Optimization_ , 28(2): 1205–1231.
  * Gao, Koppel, and Ribeiro (2020) Gao, Z.; Koppel, A.; and Ribeiro, A. 2020.  Incremental greedy BFGS: An incremental quasi-Newton method with explicit superlinear rate.  In _Advanced Neural Information Processing System 12th OPT Workshop Optimization on Machine Learning_.
  * Gay (1979) Gay, D. M. 1979.  Some convergence properties of Broyden’s method.  _SIAM Journal on Numerical Analysis_ , 16(4): 623–630.
  * Goldfarb (1970) Goldfarb, D. 1970.  A family of variable-metric methods derived by variational means.  _Mathematics of computation_ , 24(109): 23–26.
  * Gower, Goldfarb, and Richtárik (2016) Gower, R.; Goldfarb, D.; and Richtárik, P. 2016.  Stochastic block BFGS: Squeezing more curvature out of data.  In _International Conference on Machine Learning_ , 1869–1878. PMLR.
  * Gower and Richtárik (2017) Gower, R. M.; and Richtárik, P. 2017.  Randomized quasi-Newton updates are linearly convergent matrix inversion algorithms.  _SIAM Journal on Matrix Analysis and Applications_ , 38(4): 1380–1409.
  * Hastie et al. (2009) Hastie, T.; Tibshirani, R.; Friedman, J. H.; and Friedman, J. H. 2009.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer.
  * Jin and Mokhtari (2023) Jin, Q.; and Mokhtari, A. 2023.  Non-asymptotic superlinear convergence of standard quasi-Newton methods.  _Mathematical Programming_ , 200(1): 425–473.
  * Johnson and Zhang (2013) Johnson, R.; and Zhang, T. 2013.  Accelerating stochastic gradient descent using predictive variance reduction.  In _Twenty-sixth Conference on Neural Information Processing Systems_.
  * Lahoti et al. (2023) Lahoti, A.; Senapati, S.; Rajawat, K.; and Koppel, A. 2023.  Sharpened Lazy Incremental Quasi-Newton Method.  _arXiv preprint arXiv:2305.17283_.
  * Lin, Ye, and Zhang (2022) Lin, D.; Ye, H.; and Zhang, Z. 2022.  Explicit convergence rates of greedy and random quasi-Newton methods.  _The Journal of Machine Learning Research_ , 23(1): 7272–7311.
  * Liu et al. (2022) Liu, C.; Bi, S.; Luo, L.; and Lui, J. C. 2022.  Partial-Quasi-Newton Methods: Efficient Algorithms for Minimax Optimization Problems with Unbalanced Dimensionality.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , 1031–1041.
  * Liu, Chen, and Luo (2023) Liu, C.; Chen, C.; and Luo, L. 2023.  Symmetric Rank-$k$ Methods.  _arXiv preprint arXiv:2303.16188_.
  * Liu et al. (2023) Liu, C.; Chen, C.; Luo, L.; and Lui, J. C. 2023.  Block Broyden’s Methods for Solving Nonlinear Equations.  In _Thirty-seventh Conference on Neural Information Processing Systems_.
  * Liu and Luo (2022) Liu, C.; and Luo, L. 2022.  Quasi-Newton Methods for Saddle Point Problems.  In _Thirty-fifth Conference on Neural Information Processing Systems_.
  * Lucchi, McWilliams, and Hofmann (2015) Lucchi, A.; McWilliams, B.; and Hofmann, T. 2015.  A variance reduced stochastic Newton method.  _arXiv preprint arXiv:1503.08316_.
  * Mokhtari, Eisen, and Ribeiro (2018) Mokhtari, A.; Eisen, M.; and Ribeiro, A. 2018.  IQN: An incremental quasi-Newton method with local superlinear convergence rate.  _SIAM Journal on Optimization_ , 28(2): 1670–1698.
  * Mokhtari and Ribeiro (2014) Mokhtari, A.; and Ribeiro, A. 2014.  RES: Regularized stochastic BFGS algorithm.  _IEEE Transactions on Signal Processing_ , 62(23): 6089–6104.
  * Mokhtari and Ribeiro (2015) Mokhtari, A.; and Ribeiro, A. 2015.  Global convergence of online limited memory BFGS.  _The Journal of Machine Learning Research_ , 16(1): 3151–3181.
  * Moritz, Nishihara, and Jordan (2016) Moritz, P.; Nishihara, R.; and Jordan, M. 2016.  A linearly-convergent stochastic L-BFGS algorithm.  In _Artificial Intelligence and Statistics_ , 249–258. PMLR.
  * Murphy (2012) Murphy, K. P. 2012.  _Machine learning: a probabilistic perspective_.  MIT press.
  * Nesterov (2003) Nesterov, Y. 2003.  _Introductory lectures on convex optimization: A basic course_ , volume 87.  Springer Science & Business Media.
  * Nocedal and Wright (1999) Nocedal, J.; and Wright, S. J. 1999.  _Numerical optimization_.  Springer.
  * O’Leary and Yeremin (1994) O’Leary, D. P.; and Yeremin, A. 1994.  The linear algebra of block quasi-Newton algorithms.  _Linear Algebra and its Applications_ , 212: 153–168.
  * Powell (1971) Powell, M. J. 1971.  On the convergence of the variable metric algorithm.  _IMA Journal of Applied Mathematics_ , 7(1): 21–36.
  * Rodomanov and Kropotov (2016) Rodomanov, A.; and Kropotov, D. 2016.  A superlinearly-convergent proximal Newton-type method for the optimization of finite sums.  In _International Conference on Machine Learning_ , 2597–2605. PMLR.
  * Rodomanov and Nesterov (2021a) Rodomanov, A.; and Nesterov, Y. 2021a.  Greedy quasi-Newton methods with explicit superlinear convergence.  _SIAM Journal on Optimization_ , 31(1): 785–811.
  * Rodomanov and Nesterov (2021b) Rodomanov, A.; and Nesterov, Y. 2021b.  New results on superlinear convergence of classical quasi-Newton methods.  _Journal of optimization theory and applications_ , 188: 744–769.
  * Rodomanov and Nesterov (2021c) Rodomanov, A.; and Nesterov, Y. 2021c.  Rates of superlinear convergence for classical quasi-Newton methods.  _Mathematical Programming_ , 1–32.
  * Schmidt, Le Roux, and Bach (2017) Schmidt, M.; Le Roux, N.; and Bach, F. 2017.  Minimizing finite sums with the stochastic average gradient.  _Mathematical Programming_ , 162: 83–112.
  * Schnabel (1983) Schnabel, R. B. 1983.  Quasi-Newton methods using multiple secant equations.  _Computer Science Technical Reports_ , 244(41): 06.
  * Shanno (1970) Shanno, D. F. 1970.  Conditioning of quasi-Newton methods for function minimization.  _Mathematics of computation_ , 24(111): 647–656.
  * Wang et al. (2017) Wang, X.; Ma, S.; Goldfarb, D.; and Liu, W. 2017.  Stochastic quasi-Newton methods for nonconvex stochastic optimization.  _SIAM Journal on Optimization_ , 27(2): 927–956.
  * Woodworth and Srebro (2016) Woodworth, B. E.; and Srebro, N. 2016.  Tight complexity bounds for optimizing composite objectives.  In _Twenty-ninth Conference on Neural Information Processing Systems_.
  * Yang et al. (2021) Yang, M.; Milzarek, A.; Wen, Z.; and Zhang, T. 2021.  A stochastic extra-step quasi-Newton method for nonsmooth nonconvex optimization.  _Mathematical Programming_ , 1–47.
  * Ye et al. (2021) Ye, H.; Lin, D.; Zhang, Z.; and Chang, X. 2021.  Explicit superlinear convergence rates of the SR1 algorithm.  _arXiv preprint arXiv:2105.07162_.
  * Zhang, Mahdavi, and Jin (2013) Zhang, L.; Mahdavi, M.; and Jin, R. 2013.  Linear convergence with condition number independent access of full gradients.  In _Twenty-sixth Conference on Neural Information Processing Systems_.

The appendix is organized as below, Section A introduces several key lemmas
which are essential for the convergence analysis of the LISR-1 and LISR-$k$.
Section B presents the proof of lemmas and theorems introduced in Section 4
and Section 5. Section C shows the efficient implementation of the LISR-1 and
LISR-$k$ method which achieve a computation complexity of
${\mathcal{O}}(d^{2})$ at every iteration. Section D presents more details of
the dataset. It also presents additional comparison results on the quadratic
function minimization task and general function minimization task.

##  Appendix A Established Results

In this section, we revisit some key lemmas which are essential for the
analysis of the proposed algorithms.

######  Lemma A.1 (Rodomanov and Nesterov (2021a)).

Let $f$ be a strongly self-concordant function with some constant $M$, and let
$r\coloneqq\left\lVert y-x\right\rVert_{\nabla^{2}f(x)}$ for any
$x,y\in{\mathcal{R}}^{d}$. Then,

$\frac{\nabla^{2}f(x)}{1+Mr}\preceq\nabla^{2}f(y)\preceq(1+Mr)\nabla^{2}f(x).$
(25)

Also, for $J\coloneqq\int_{0}^{1}\nabla^{2}f(x+t(y-x))dt$, we have

$\frac{\nabla^{2}f(x)}{1+\frac{Mr}{2}}\preceq
J\preceq\left(1+\frac{Mr}{2}\right)\nabla^{2}f(x),$ (26)
$\frac{\nabla^{2}f(y)}{1+\frac{Mr}{2}}\preceq
J\preceq\left(1+\frac{Mr}{2}\right)\nabla^{2}f(y).$ (27)

######  Lemma A.2 (Banach’s Lemma).

Let $A\in{\mathbb{R}}^{d\times d}$ be a matrix such that its norm satisfies
$\left\lVert A\right\rVert<1$. Then the matrix $(I+A)$ is invertible and

$\frac{1}{1+\left\lVert
A\right\rVert}<\left\lVert(I+A)^{-1}\right\rVert<\frac{1}{1-\left\lVert
A\right\rVert}.$

######  Lemma A.3 (Ye et al. (2021)).

Consider positive definite matrices $A,G\in{\mathbb{R}}^{d\times d}$ and
suppose $G_{+}\coloneqq SR1(G,A,u)$, where $u\neq 0$. Suppose for some
$\eta>1$, we have

$A\preceq G\preceq\eta A,$

then, for any $u\neq 0$, we also have:

$A\preceq G_{+}\preceq G\preceq\eta A.$

######  Lemma A.4 (Lin, Ye, and Zhang (2022)).

Consider positive definite matrices $A,G\in{\mathbb{R}}^{d\times d}$ such that
$A\preceq G$. Suppose that $G_{+}\coloneqq\textrm{SR1}(G,A,\tilde{u}(G,A))$,
where $\tilde{u}(G,A)$ (equation (9)) is the greedy vector of $G$ with respect
to $A$. Then, the following holds:

$\tau(G_{+},A)\leq\left(1-\frac{1}{d}\right)\tau(G,A).$

######  Lemma A.5 (Liu, Chen, and Luo (2023)).

For any positive-definite matrices $A\in{\mathbb{R}}^{d\times d}$ and
$G\in{\mathbb{R}}^{d\times d}$ with $A\preceq G\preceq\eta A$ for some
$\eta\geq 1$, we let $G_{+}=\mathrm{SR}\textrm{-}k(G,A,U)$ for some full rank
matrix $U\in{\mathbb{R}}^{d\times k}$. Then it holds that

$A\preceq G_{+}\preceq\eta A.$ (28)

######  Lemma A.6 (Liu, Chen, and Luo (2023)).

Suppose that

$G_{+}=\mathrm{SR}\textrm{-}k(G,A,U),$

where $G\succeq A\in{\mathbb{R}}^{d\times d}$ and select
$U\in{\mathbb{R}}^{d\times k}$ by equation (20). Then we have

$\tau(G_{+},A)\leq\left(1-\frac{k}{d}\right)\tau(G,A).$ (29)

##  Appendix B Convergence Analysis of the LISR-1/LISR-$k$ method

In this section, we present the theoretical analysis of the Theorem 18 and
Theorem 22. To facilitate the analysis of main theorems, we introduce the
following quantity given two PSD matrices $G$ and $A$ such that $G\succeq A$:

$\tau(G,A)=\rm tr(G-A).$

###  B.1 Supporting Lemmas

The norm of the difference between two PSD matrices can be upper bounded as
below:

######  Lemma B.1.

For all positive definite matrices $A,G\in{\mathbb{R}}^{d\times d}$, if
$A\preceq G$ and $A\preceq LI$, then

$\left\lVert
G-A\right\rVert\leq\frac{dL\tau(G,A)}{\mathrm{tr}(A)}=\mu\nu(G,A).$

###### Proof.

For any PSD matrix $X$, let $\lambda_{max}(X)$ denote the $i$-th largest
eigenvalue of $X$, we have

$\displaystyle\left\lVert
G-A\right\rVert=\lambda_{1}(G-A)\leq\sum_{i=1}^{d}\lambda_{i}(G-A)=tr(G-A)=\tau(G,A).$

Furthermore, we can bound $\mathrm{tr}(A)$ by:

$\mathrm{tr}(A)\leq dL.$

Therefore, we have:

$\frac{dL\tau(G,A)}{\mathrm{tr}(A)}\geq\frac{dL\left\lVert
G-A\right\rVert}{dL}=\left\lVert G-A\right\rVert.$

∎

Now we show how the Hessian approximation metric $\nu(B,\nabla^{2}f(x))$
changes after applying one step of greedy SR1 update:

######  Lemma B.2.

Let $f:{\mathbb{R}}^{d}\to{\mathbb{R}}$ be a real-valued function that is
$\mu$-strongly convex, $L$-smooth, and $M$-strongly self-concordant. Let
$x,x_{+}\in{\mathbb{R}}^{d}\setminus\\{0\\}$ and $B$ be a matrix such that
$B\geq\nabla^{2}f(x)$. Define the constant $r\coloneqq\left\lVert
x_{+}-x\right\rVert_{\nabla^{2}f(x)}$ and the matrix
$P\coloneqq\left(1+\frac{Mr}{2}\right)^{2}B$. Consider the following SR-$k$
updates:

$B_{+}\coloneqq\mathrm{SR}\textrm{-}k(P,\nabla^{2}f(x_{+}),\bar{U}(P,\nabla^{2}f(x_{+}))).$

Here the vector $\bar{U}(P,\nabla^{2}f(x_{+}))$ is the greedy matrx defined in
equation (20). Then $B_{+}\succeq\nabla^{2}f(x_{+})$ and

$\frac{\tau(B_{+},\nabla^{2}f(x_{+}))}{\mathrm{tr}(\nabla^{2}f(x_{+}))}\leq\left(1-\frac{k}{d}\right)\left(1+\frac{Mr}{2}\right)^{4}\left(\frac{\tau(B,\nabla^{2}f(x))}{\mathrm{tr}(\nabla^{2}f(x))}+2Mr\right).$

###### Proof.

Since $B\succeq\nabla^{2}f(x)$, we have:

$\displaystyle P=$
$\displaystyle(1+\frac{Mr}{2})^{2}B\succeq(1+\frac{Mr}{2})^{2}\nabla^{2}f(x)\succeq(1+\frac{Mr}{2})K\succeq\nabla^{2}f(x_{+}),$

where the matrix $K\coloneqq\int_{\tau=0}^{1}f(x+\tau(x_{+}-x))d\tau$. The
first inequality follows from the problem assumption and the last two
inequalities are due to Lemma 27.

Apply Lemma 29, and we have:

$\displaystyle\tau(B_{+},\nabla^{2}f(x_{+}))\preceq$
$\displaystyle\left(1-\frac{k}{d}\right)\tau(P,\nabla^{2}f(x_{+})).$

The last term can be bounded as

$\displaystyle\mathrm{tr}(P-\nabla^{2}f(x_{+}))$
$\displaystyle\overset{(\ref{matrix_approx_1})}{\leq}$
$\displaystyle\mathrm{tr}\left(\left(1+\frac{Mr}{2}\right)^{2}B-\frac{\nabla^{2}f(x)}{1+Mr}\right)$
$\displaystyle=$
$\displaystyle\left(1+\frac{Mr}{2}\right)^{2}\tau(B,\nabla^{2}f(x))+\left(\left(1+\frac{Mr}{2}\right)^{2}-\frac{1}{1+Mr}\right)\mathrm{tr}(\nabla^{2}f(x))$
$\displaystyle\leq$
$\displaystyle\left(1+\frac{Mr}{2}\right)^{2}\tau(B,\nabla^{2}f(x))+\left(\left(1+\frac{Mr}{2}\right)^{2}+Mr-1\right)\mathrm{tr}(\nabla^{2}f(x))$
$\displaystyle=$
$\displaystyle\left(1+\frac{Mr}{2}\right)^{2}\tau(B,\nabla^{2}f(x))+Mr\left(2+\frac{Mr}{4}\right)\mathrm{tr}(\nabla^{2}f(x))$
$\displaystyle\overset{(\ref{matrix_approx_1})}{\leq}$
$\displaystyle\left(1+\frac{Mr}{2}\right)^{2}(1+Mr)\left(\frac{\tau(B,\nabla^{2}f(x))}{\mathrm{tr}(\nabla^{2}f(x))}+2Mr\right)\mathrm{tr}(\nabla^{2}f(x_{+})),$

By $1+Mr\leq(1+\frac{Mr}{2})^{2}$, and rearrange the terms we have:

$\frac{\tau(B_{+},\nabla^{2}f(x_{+}))}{\mathrm{tr}(\nabla^{2}f(x_{+}))}\leq
c\left(1+\frac{Mr}{2}\right)^{4}\left(\frac{\tau(B,\nabla^{2}f(x))}{\mathrm{tr}(\nabla^{2}f(x))}+2Mr\right),$

where $c\coloneqq 1-kd^{-1}$. ∎

Applying the result of Lemma B.2, we can upper bound the Hessian approximation
error with the following lemma:

######  Lemma B.3.

Let $f:{\mathbb{R}}^{d}\to{\mathbb{R}}$ be a real-valued function that is
$\mu$-strongly convex, $L$-smooth, and $M$-strongly self-concordant. Let
$\tilde{x}\in{\mathbb{R}}^{d}$ be some fixed vector and $0\leq\gamma<1$ be
some fixed constant such that the sequence $\\{x^{t}\\}$, for all $t\in[T]$
satisfies

$\left\lVert x^{t}-\tilde{x}\right\rVert\leq\gamma^{t}\left\lVert
x^{0}-\tilde{x}\right\rVert.$ (30)

Define the constant $r_{t}\coloneqq\left\lVert
x^{t}-x^{t-1}\right\rVert_{\nabla^{2}f(x^{t-1})}$ for every $t$. Let $B^{0}$
be a matrix such that it satisfies $B_{0}\succeq\nabla^{2}f(x^{0})$. Consider
the following SR-$k$ update:

$\displaystyle B^{t}\coloneqq$
$\displaystyle\mathrm{SR}\textrm{-}k(P^{t-1},\nabla^{2}f(x^{t}),\tilde{U}(P^{t-1},\nabla^{2}f(x^{t}))),$

where $P^{t-1}\coloneqq(1+\frac{Mr_{t}}{2})^{2}B^{t-1}$, and
$\tilde{U}(P^{t-1},\nabla^{2}f(x^{t}))$ is the greedy matrix (20). Then the
following holds for all $t\in[T]$:

$\nu(B^{t},\nabla^{2}f(x^{t}))\leq\left(1-\frac{k}{d}\right)^{t}e^{\frac{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert}{1-\gamma}}\left(\nu(B^{0},\nabla^{2}f(x^{0}))+\left\lVert
x^{0}-\tilde{x}\right\rVert\frac{4dML^{\frac{3}{2}}\mu^{-1}}{1-(1-k/d)^{-1}\gamma}\right).$
(31)

###### Proof.

From Lemma B.2, it can be shown that $B^{t}\succeq\nabla^{2}f(x^{t})$ for
$t\in[T]$. Therefore, $\tau(B^{t},\nabla^{2}f(x^{t}))$ and
$\nu(B^{t},\nabla^{2}f(x^{t}))$ are both well defined.

First, we can apply the triangle inequality on $r_{t}$ to obtain:

$\displaystyle r_{t}=\left\lVert
x^{t}-x^{t-1}\right\rVert_{\nabla^{2}f(x^{t-1})}\leq$
$\displaystyle\sqrt{L}\left\lVert x^{t}-x^{t-1}\right\rVert$
$\displaystyle\leq$ $\displaystyle\sqrt{L}(\left\lVert
x^{t}-\tilde{x}\right\rVert+\left\lVert x^{t-1}-\tilde{x}\right\rVert)\leq
2\sqrt{L}\gamma^{t-1}\left\lVert x^{0}-\tilde{x}\right\rVert,$

Where the first inequality is due to the $L$-smoothness of function
$f(\cdot)$. The last inequality follows from equation (30). Applying Lemma B.2
for the SR-$k$ update, we have:

$\displaystyle\frac{\tau(B^{t},\nabla^{2}f(x^{t}))}{\mathrm{tr}(\nabla^{2}f(x^{t}))}\leq$
$\displaystyle
c\left(1+\frac{Mr_{t-1}}{2}\right)^{4}\left(\frac{\tau(B^{t-1},\nabla^{2}f(x^{t-1}))}{\mathrm{tr}(\nabla^{2}f(x^{t-1}))}+2Mr_{t-1}\right)$
$\displaystyle\leq$ $\displaystyle
ce^{2Mr_{t-1}}\left(\frac{\tau(B^{t-1},\nabla^{2}f(x^{t-1}))}{\mathrm{tr}(\nabla^{2}f(x^{t-1}))}+2Mr_{t-1}\right)$
$\displaystyle\leq$ $\displaystyle ce^{4M\sqrt{L}\gamma^{t-1}\left\lVert
x^{0}-\tilde{x}\right\rVert}\left(\frac{\tau(B^{t-1},\nabla^{2}f(x^{t-1}))}{\mathrm{tr}(\nabla^{2}f(x^{t-1}))}+4M\sqrt{L}\gamma^{t-1}\left\lVert
x^{0}-\tilde{x}\right\rVert\right),$

where $c=1-kd^{-1}$ for the SR-$k$ update. The second inequality follows from
$1+x\leq e^{x}$. Expand the recursion, then:

$\displaystyle\frac{\tau(B^{t},\nabla^{2}f(x^{t}))}{\mathrm{tr}(\nabla^{2}f(x^{t}))}\leq$
$\displaystyle ce^{4M\sqrt{L}\gamma^{t-1}\left\lVert
x^{0}-\tilde{x}\right\rVert}\left(\frac{\tau(B^{t-1},\nabla^{2}f(x^{t-1}))}{\mathrm{tr}(\nabla^{2}f(x^{t-1}))}+4M\sqrt{L}\gamma^{t-1}\left\lVert
x^{0}-\tilde{x}\right\rVert\right)$ $\displaystyle\leq$ $\displaystyle
c^{2}e^{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert(\gamma^{t-1}+\gamma^{t-2})}\frac{\tau(B^{t-2},\nabla^{2}f(x^{t-2}))}{\mathrm{tr}(\nabla^{2}f(x^{t-2}))}+$
$\displaystyle 4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\left(c\gamma^{t-1}e^{4M\sqrt{L}\gamma^{t-1}\left\lVert
x^{0}-\tilde{x}\right\rVert}+c^{2}\gamma^{t-2}e^{4M\sqrt{L}(\gamma^{t-2}+\gamma^{t-1})\left\lVert
x^{0}-\tilde{x}\right\rVert}\right)$ $\displaystyle\leq$ $\displaystyle
c^{t}e^{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\sum_{j=0}^{t-1}\gamma^{j}}\frac{\tau(B^{0},\nabla^{2}f(x^{0}))}{\mathrm{tr}(\nabla^{2}f(x^{0}))}+$
$\displaystyle 4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\left(\sum_{j=0}^{t-1}c^{t-j}\gamma^{j}e^{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\sum_{i=1}^{t-j}\gamma^{t-i}}\right)$
$\displaystyle\leq$ $\displaystyle c^{t}e^{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\sum_{j=0}^{+\infty}\gamma^{j}}\frac{\tau(B^{0},\nabla^{2}f(x^{0}))}{\mathrm{tr}(\nabla^{2}f(x^{0}))}+$
$\displaystyle 4M\sqrt{L}\left\lVert x^{0}-\tilde{x}\right\rVert
e^{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert\sum_{i=1}^{+\infty}\gamma^{i}}\left(\sum_{j=0}^{t-1}c^{t-j}\gamma^{j}\right)$
$\displaystyle\leq$ $\displaystyle c^{t}e^{\frac{4M\sqrt{L}\left\lVert
x^{0}-\tilde{x}\right\rVert}{1-\gamma}}\left(\frac{\tau(B^{0},\nabla^{2}f(x^{0}))}{\mathrm{tr}(\nabla^{2}f(x^{0}))}+\left\lVert
x^{0}-\tilde{x}\right\rVert\frac{4M\sqrt{L}}{1-c^{-1}\gamma}\right).$

Multiply both sides of the above inequality by $d\kappa$ and we get the
desired result. ∎

###  B.2 Proof of Lemma 4.1

For the completeness of the paper presentation, we show the proof of Lemma 4.1
below:

###### Proof.

For all $t\geq 0$, we define $H^{t}\coloneqq(\sum_{i=1}^{n}B_{i}^{t})^{-1}$.
From the update for $x^{t+1}$, we have:

$\displaystyle\left\lVert x^{t+1}-x^{*}\right\rVert=$ $\displaystyle
H^{t}\left(\sum_{i=1}^{n}B_{i}^{t}z_{i}^{t}-\sum_{i=1}^{n}\nabla
f_{i}(z_{i}^{t})\right)-x^{*}$ $\displaystyle=$ $\displaystyle
H^{t}\left(\sum_{i=1}^{n}B_{i}^{t}(z_{i}^{t}-x^{*})-\sum_{i=1}^{n}\nabla
f_{i}(z_{i}^{t})\right)$ $\displaystyle=$ $\displaystyle
H^{t}\left(\sum_{i=1}^{n}B_{i}^{t}(z_{i}^{t}-x^{*})-\sum_{i=1}^{n}(\nabla
f_{i}(z_{i}^{t})-\nabla f_{i}(x^{*}))\right)$ $\displaystyle=$ $\displaystyle
H^{t}\left(\sum_{i=1}^{n}B_{i}^{t}(z_{i}^{t}-x^{*})-\sum_{i=1}^{n}\int_{0}^{1}\nabla^{2}f(x^{*}+(z_{i}^{t}-x^{*})v)(z_{i}^{t}-x^{*})dv\right)$
$\displaystyle=$ $\displaystyle
H^{t}\left(\sum_{i=1}^{n}(B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t}))(z_{i}^{t}-x^{*})+\sum_{i=1}^{n}(\nabla^{2}f_{i}(z_{i}^{t})-\int_{0}^{1}\nabla^{2}f(x^{*}+(z_{i}^{t}-x^{*})v))(z_{i}^{t}-x^{*})dv\right),$

where the second last equality follows from the Fundamental Theorem of
Calculus. Taking the norm on both sides and applying the triangle inequality,
we have

$\displaystyle\left\lVert x^{t+1}-x^{*}\right\rVert\leq$
$\displaystyle\left\lVert H^{t}\right\rVert\left(\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert+\sum_{i=1}^{n}\left\lVert\int_{0}^{1}(\nabla^{2}f_{i}(z_{i}^{t})-\nabla^{2}f(x^{*}+(z_{i}^{t}-x^{*})v))(z_{i}^{t}-x^{*})dv\right\rVert\right)$
$\displaystyle\leq$ $\displaystyle\left\lVert
H^{t}\right\rVert\left(\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert+\sum_{i=1}^{n}\int_{0}^{1}\left\lVert(\nabla^{2}f_{i}(z_{i}^{t})-\nabla^{2}f(x^{*}+(z_{i}^{t}-x^{*})v))(z_{i}^{t}-x^{*})\right\rVert
dv\right)$ $\displaystyle\leq$ $\displaystyle\left\lVert
H^{t}\right\rVert\left(\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert+\sum_{i=1}^{n}\int_{0}^{1}\left\lVert\nabla^{2}f_{i}(z_{i}^{t})-\nabla^{2}f(x^{*}+(z_{i}^{t}-x^{*})v)\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert dv\right)$ $\displaystyle\leq$
$\displaystyle\left\lVert H^{t}\right\rVert\left(\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert+\sum_{i=1}^{n}\tilde{L}\int_{0}^{1}(1-v)dv\left\lVert
z_{i}^{t}-x^{*}\right\rVert^{2}\right)$ $\displaystyle\leq$
$\displaystyle\left\lVert H^{t}\right\rVert\left(\sum_{i=1}^{n}\left\lVert
B_{i}^{t}-\nabla^{2}f_{i}(z_{i}^{t})\right\rVert\left\lVert
z_{i}^{t}-x^{*}\right\rVert+\frac{\tilde{L}}{2}\sum_{i=1}^{n}\left\lVert
z_{i}^{t}-x^{*}\right\rVert^{2}\right),$

where the second inequality follows from the result that if
$g:{\mathbb{R}}\to{\mathbb{R}}^{d}$ is a continuous function, then
$\left\lVert\int_{0}^{1}g(v)dv\right\rVert\leq\int_{0}^{1}\left\lVert
g(v)\right\rVert dv$, and the fourth inequality follows the assumption that
the Hessian of $f_{i}$ is $\tilde{L}$-Lipschitz. ∎

###  B.3 Proof of Lemma 4.3

We present and prove a more generalized version of Lemma 4.3:

######  Lemma B.4.

We initialize each Hessian approximation
$B_{i}^{0}=(1+M\sqrt{L}r_{0})^{2}E_{i}^{0}$ where $E_{i}^{0}$ is some PSD
matrix that satisfies $E_{i}^{0}\succeq\nabla^{2}f_{i}(x^{0})$. For any $\rho$
such that $0<\rho<1-\frac{k}{d}$, there exists positive constants $r_{0}$ and
$\sigma_{0}$ such that if $\left\lVert x^{0}-x^{*}\right\rVert\leq r_{0}$ and
$\nu(E_{i}^{0},\nabla^{2}f_{i}(x^{0}))\leq\sigma_{0}$, for all $i\in[n]$, the
sequence of iterates generated by the LISR-$k$ satisfies

$\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\rho^{\left\lceil\frac{t+1}{n}\right\rceil}\left\lVert
x^{0}-x^{*}\right\rVert.$ (32)

Furthermore, it holds that

$\nu((\omega^{t+1})^{-1}B_{i_{t}}^{t+1},\nabla^{2}f_{i_{t}}(z_{i_{t}}^{t+1}))\leq\left(1-\frac{k}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\delta,$
(33)

where $\delta\coloneqq
e^{\frac{4M\sqrt{L}r_{0}}{1-\rho}}\big{(}\sigma_{0}+r_{0}\frac{4MdL^{\frac{3}{2}}\mu^{-1}}{1-(1-kd^{-1})^{-1}\rho}\big{)}$,
$M=\tilde{L}/\mu^{\frac{3}{2}}$, $\omega^{t}=(1+\alpha_{\left\lceil
t/n\right\rceil})^{2}$ if $t$ is a multiple of $n$ and 1 otherwise, and the
sequence $\\{\alpha_{k}\\}$ is defined as $\alpha_{k}\coloneqq
M\sqrt{L}r_{0}\rho^{k}$, $\forall k\geq 0$.

###### Proof.

Define $c\coloneqq 1-kd^{-1}$.

For a given $\rho$ that satisfies $0<\rho<c$, we choose $r_{0},\delta>0$ such
that they satisfy

$\frac{\tilde{L}\mu^{-1}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}\delta+ML^{\frac{3}{2}}\mu^{-1}r_{0}(2+M\sqrt{L}r_{0})\leq\frac{\rho}{1+\rho}.$
(34)

Base case: At $t=1$, from Lemma 4.1, we have

$\left\lVert
x^{1}-x^{*}\right\rVert\leq\frac{\tilde{L}\Gamma^{0}}{2}\sum_{i=1}^{n}\left\lVert
z_{i}^{0}-x^{*}\right\rVert^{2}+\Gamma^{0}\sum_{i=1}^{n}\left\lVert
B_{i}^{0}-\nabla^{2}f_{i}(z_{i}^{0})\right\rVert\left\lVert
z_{i}^{0}-x^{*}\right\rVert.$

For the initialization, we have that $B_{i}^{0}=(1+\alpha_{0})^{2}E_{i}^{0}$
and $z_{i}^{0}=x^{0}$, for all $i\in[n]$, and $\left\lVert
x^{0}-x^{*}\right\rVert\leq r_{0}$. Substituting these in the above
expression, we obtain:

$\displaystyle\left\lVert x^{1}-x^{*}\right\rVert\leq$
$\displaystyle\Gamma^{0}\left(n\frac{\tilde{L}r_{0}}{2}+\sum_{i=1}^{n}\left\lVert(1+\alpha_{0})^{2}E_{i}^{0}-\nabla^{2}f_{i}(z_{i}^{0})\right\rVert\right)\left\lVert
x^{0}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{0}\left(n\frac{\tilde{L}r_{0}}{2}+(1+\alpha_{0})^{2}\sum_{i=1}^{n}\left\lVert
E_{i}^{0}-\nabla^{2}f_{i}(z_{i}^{0})\right\rVert+\alpha_{0}(\alpha_{0}+2)\sum_{i=1}^{n}\left\lVert\nabla^{2}f_{i}(z_{i}^{0})\right\rVert\right)\left\lVert
x^{0}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{0}\left(n\frac{\tilde{L}r_{0}}{2}+(1+\alpha_{0})^{2}\sum_{i=1}^{n}\frac{dL\tau(E_{i}^{0},\nabla^{2}f_{i}(x^{0}))}{\mathrm{tr}(\nabla^{2}f_{i}(x^{0}))}+nL\alpha_{0}(\alpha_{0}+2)\right)\left\lVert
x^{0}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{0}\left(n\frac{\tilde{L}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}n\mu\sigma_{0}+nML^{\frac{3}{2}}r_{0}(M\sqrt{L}r_{0}+2)\right)\left\lVert
x^{0}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{0}\left(n\frac{\tilde{L}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}n\mu\delta+nML^{\frac{3}{2}}r_{0}(M\sqrt{L}r_{0}+2)\right)\left\lVert
x^{0}-x^{*}\right\rVert.$

The second inequality follows from the triangle inequality and the third
inequality is due to Lemma B.1 and the initialization condition
$\tau(E_{i}^{0},\nabla^{2}f_{i}(z_{i}^{0}))\leq\sigma_{0}$.

We now upper bound $\Gamma^{0}$. Define
$X^{0}\coloneqq\frac{1}{n}\sum_{i=1}^{n}B_{i}^{0}$ and
$Y^{0}\coloneqq\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}(z_{i}^{0})$. Then we
have

$\frac{1}{n}\sum_{i=1}^{n}\left\lVert
B_{i}^{0}-\nabla^{2}f_{i}(z_{i}^{0})\right\rVert\geq\left\lVert
X^{0}-Y^{0}\right\rVert=\left\lVert(Y^{0})((Y^{0})^{-1}X^{0}-I)\right\rVert\geq\mu\left\lVert(Y^{0})^{-1}X^{0}-I\right\rVert,$

where the second inequality follows from each $f_{i}$ is $\mu$-strongly
convex. By tracking the steps for deriving the bound of $\left\lVert
x^{1}-x^{*}\right\rVert$, we have

$\sum_{i=1}^{n}\left\lVert
B_{i}^{0}-\nabla^{2}f_{i}(z_{i}^{0})\right\rVert\leq(1+M\sqrt{L}r_{0})^{2}n\mu\delta+nML^{\frac{3}{2}}r_{0}(M\sqrt{L}r_{0}+2).$

Combining the above two inequalities, we obtain

$\left\lVert(Y^{0})^{-1}X^{0}-I\right\rVert\leq(1+M\sqrt{L}r_{0})^{2}\delta+\frac{ML^{\frac{3}{2}}r_{0}}{\mu}(M\sqrt{L}r_{0}+2)<\frac{\rho}{1+\rho}.$

We can now upper bound $\Gamma^{0}$ using Lemma A.2. Since the matrix
$(Y^{0})^{-1}X^{0}-I\succ 0$, we can deduce that

$\left\lVert(X^{0})^{-1}Y^{0}\right\rVert=\left\lVert(I+((Y^{0})^{-1}X^{0})-I)^{-1}\right\rVert\leq\frac{1}{1-\left\lVert(Y^{0})^{-1}X^{0}-I\right\rVert}\leq
1+\rho.$

Recall that $\mu I\preceq Y^{0}$. Consequently, we can bound $\Gamma^{0}$ with

$\Gamma^{0}=\frac{1}{n}\left\lVert(X^{0})^{-1}\right\rVert\leq\frac{1+\rho}{n\mu}.$

Therefore, using equation (34), we have

$\left\lVert
x^{1}-x^{*}\right\rVert\leq(1+\rho)\left(\frac{\tilde{L}\mu^{-1}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}\delta+ML^{\frac{3}{2}}\mu^{-1}r_{0}(M\sqrt{L}r_{0}+2)\right)\left\lVert
x^{0}-x^{*}\right\rVert\leq\rho\left\lVert x^{0}-x^{*}\right\rVert.$

To complete the base step, we now upper bound
$\nu(\omega_{1}^{-1}B_{1}^{1},\nabla^{2}f_{1}(z_{1}^{1}))$, where
$\omega_{1}=1$. Applying Lemma 31 with parameters as $T=1$, $\tilde{x}=x^{*}$,
$P^{0}=(1+\alpha_{0})^{2}I_{1}^{0}=B_{1}^{0}$, we get

$\displaystyle\nu(B^{1},\nabla^{2}f(x^{1}))\leq$ $\displaystyle
ce^{\frac{4M\sqrt{L}r_{0}}{1-\rho}}\left(\nu(B^{0},\nabla^{2}f(x^{0}))+\frac{4Md\kappa\sqrt{L}r_{0}}{1-c^{-1}\rho}\right)$
$\displaystyle\leq$ $\displaystyle
ce^{\frac{4M\sqrt{L}r_{0}}{1-\rho}}\left(\sigma^{0}+\frac{4MdL^{\frac{3}{2}}\mu^{-1}r_{0}}{1-c^{-1}\rho}\right)$
$\displaystyle\leq$ $\displaystyle c\delta.$

This completes the proof for the base case.

Induction Hypothesis: Let equation (16) and equation (17) hold for all
$t\in[jn+m-1]$ for some $j\geq 0$ and $0\leq m<n$.

Induction Step: We then prove that equation (16) and equation (17) also hold
for $t=jn+m$. Recall that the tuples are updated in a deterministic cyclic
order, and at the current time $t$, we are in the $j$th cycle and have updated
the $m$th tuple. Therefore, it is easy to note that
$z_{i}^{jn+m}=z_{i}^{jn+i}$, for all $i\in[m]$, and
$z_{i}^{jn+m}=z_{i}^{jn-n+i}$ for all $i\in[n]\setminus[m]$. From the
induction hypothesis, we have

$\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert\leq\begin{cases}\rho^{\left\lceil\frac{jn+i}{n}\right\rceil}\left\lVert
x^{0}-x^{*}\right\rVert,&i\in[m]\\\
\rho^{\left\lceil\frac{(j-1)n+i}{n}\right\rceil}\left\lVert
x^{0}-x^{*}\right\rVert,&i\in[n]\setminus[m]\end{cases}$

#### Upper Bound of $\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$.

We will establish an upper bound on $\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$. Since $B_{i_{t}}^{t}$
is updated in a different manner for $t~{}{\rm mod}~{}n\neq 0$ and $t~{}{\rm
mod}~{}n=0$, we analyze these two cases separately.

(a). $t~{}{\rm mod}~{}n\neq 0$

Since $t=jn+m+1$, this case is equivalent to considering $0\leq m<n-1$. From
the structure of the cyclic updates and equation (8), we can observe that
$B_{i}^{jn+m}=B_{i}^{jn+i}$, for all $i\in[m]$,
$B_{i}^{jn+m}=(1+M\sqrt{L}r_{0}\rho^{j})B_{i}^{jn-n+i}$, for all
$i\in[n-1]\setminus[m]$, and $B_{i}^{jn+m}=B_{i}^{jn}$, for $i=n$.

For all $i\in[m]$, from the induction hypothesis, we have

$\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert=\left\lVert
B_{i}^{jn+i}-\nabla^{2}f_{i}(z_{i}^{jn+i})\right\rVert\leq\frac{Ld\tau(B_{i}^{jn+i},\nabla^{2}f_{i}(z_{i}^{jn+i}))}{\mathrm{tr}(\nabla^{2}f_{i}(z_{i}^{jn+i}))}\leq
c^{\left\lceil\frac{jn+i}{n}\right\rceil}\mu\delta.$

The first inequality follows from Lemma B.1.

For all $i\in[n-1]\setminus[m]$, we have:

$\displaystyle\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert=$
$\displaystyle\left\lVert(1+M\sqrt{L}r_{0}\rho^{j})^{2}B_{i}^{jn-n+i}-\nabla^{2}f_{i}(x^{jn-n+i})\right\rVert$
$\displaystyle\leq$ $\displaystyle(1+M\sqrt{L}r_{0}\rho^{j})^{2}\left\lVert
B_{i}^{jn-n+i}-\nabla^{2}f_{i}(x^{jn-n+i})\right\rVert$
$\displaystyle+M\sqrt{L}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j})\left\lVert\nabla^{2}f_{i}(x^{jn-n+i})\right\rVert$
$\displaystyle\leq$ $\displaystyle(1+M\sqrt{L}r_{0}\rho^{j})^{2}\left\lVert
B_{i}^{jn-n+i}-\nabla^{2}f_{i}(x^{jn-n+i})\right\rVert+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j})$
$\displaystyle\leq$
$\displaystyle(1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}).$

The first inequality follows from the triangle inequality. The second
inequality is due to the assumption that each $f_{i}$ is $L$-smooth. The last
inequality is deduced from the induction hypothesis.

For $i=n$, we have

$\displaystyle\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert=$
$\displaystyle\left\lVert B_{i}^{jn}-\nabla^{2}f_{i}(z_{i}^{jn})\right\rVert$
$\displaystyle=$
$\displaystyle\left\lVert\omega_{jn}(\omega_{jn}^{-1}B_{i}^{jn}-\nabla^{2}f_{i}(z_{i}^{jn}))+(\omega_{jn}-1)\nabla^{2}f_{i}(z_{i}^{jn})\right\rVert$
$\displaystyle\leq$
$\displaystyle\omega_{jn}\left\lVert\omega_{jn}^{-1}B_{i}^{jn}-\nabla^{2}f_{i}(z_{i}^{jn})\right\rVert+(\omega_{jn}-1)\left\lVert\nabla^{2}f_{i}(z_{i}^{jn})\right\rVert$
$\displaystyle\leq$
$\displaystyle(1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}).$

The last inequality follows from the induction hypothesis and the assumption
that each $f_{i}$ is $L$-smooth.

We can now bound the quantity $\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$ as follows:

$\begin{split}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq&mc^{\left\lceil\frac{jn+i}{n}\right\rceil}\mu\delta+(n-m)((1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta\\\
&+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}))\\\
\leq&m\mu\delta+(n-m)((1+M\sqrt{L}r_{0})^{2}\mu\delta+ML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0}))\\\
\leq&n(1+M\sqrt{L}r_{0})^{2}\mu\delta+nML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0}).\end{split}$
(35)

(b). $t~{}{\rm mod}~{}n=0$

Since $t=jn+m+1$, we can infer that $m=n-1$. Here we have
$B_{i}^{jn}=(1+M\sqrt{L}r_{0}\rho^{j})^{2}B_{i}^{jn-n+i}$, for all
$i\in[n-1]$, and $B_{n}^{jn}$ would be used as it is.

Similar to the case $t\mod n\neq 0$, for $i\in[n-1]$, we get

$\displaystyle\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq(1+M\sqrt{L}r_{0}\rho^{j})^{2}\mu\delta
c^{j}+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}),$

For $i=n$, we have

$\displaystyle\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq(1+M\sqrt{L}r_{0}\rho^{j})^{2}\mu\delta
c^{j}+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}),$

We can now bound the term $\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$ in the following
manner:

$\begin{split}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq&n(1+M\sqrt{L}r_{0}\rho^{j})^{2}\mu\delta
c^{j}+nML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j})\\\
\leq&n(1+M\sqrt{L}r_{0})^{2}\mu\delta+nML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0}).\end{split}$
(36)

#### Upper Bound of $\left\lVert z_{m+1}^{jn+m+1}-x^{*}\right\rVert$

We now give the bound of $\left\lVert z_{m+1}^{jn+m+1}-x^{*}\right\rVert$
using the bound of $\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$:

$\displaystyle\left\lVert z_{m+1}^{jn+m+1}-x^{*}\right\rVert\leq$
$\displaystyle\frac{\tilde{L}\Gamma^{jn+m}}{2}\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert^{2}$
$\displaystyle+\Gamma^{jn+m}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{jn+m}\left(n\frac{\tilde{L}r_{0}}{2}+n(1+M\sqrt{L}r_{0})^{2}\mu\delta+nML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0})\right)\rho^{j}\left\lVert
x^{0}-x^{*}\right\rVert.$

The first inequality is due to Lemma 4.1. The last inequality follows from
equation (35), equation (36), and the induction hypothesis.

To bound $\Gamma^{jn+m}$, let
$X^{jn+m}\coloneqq\frac{1}{n}\sum_{i=1}^{n}B_{i}^{jn+m}$ and
$Y^{jn+m}\coloneqq\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}(z_{i}^{jn+m})$.
Observe that:

$\frac{1}{n}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\geq\left\lVert
X^{jn+m}-Y^{jn+m}\right\rVert\geq\mu\left\lVert(Y^{jn+m})^{-1}X^{jn+m}-I\right\rVert.$

The last inequality follows from the assumption that each $f_{i}$ is
$\mu$-strongly convex. Furthermore, using the equation (35) and equation (36),
we can see that:

$\frac{1}{n}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq(1+M\sqrt{L}r_{0})^{2}\mu\delta+ML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0}).$

Combining the above two inequalities, we can infer that:

$\left\lVert(Y^{jn+m})^{-1}X^{jn+m}-I\right\rVert\leq(1+M\sqrt{L}r_{0})^{2}\delta+ML^{\frac{3}{2}}\mu^{-1}r_{0}(2+M\sqrt{L}r_{0})\leq\frac{\rho}{1+\rho}.$

We can now bound $\Gamma^{jn+m}$ using Banach’s Lemma following the same way
we derive for the base case. We get that

$n\Gamma^{jn+m}=\left\lVert(X^{jn+m})^{-1}\right\rVert\leq\frac{1+\rho}{\mu}.$

Consequently, we have

$\begin{split}\left\lVert
z_{m+1}^{jn+m+1}-x^{*}\right\rVert\leq&\Gamma^{jn+m}\left(n\frac{\tilde{L}r_{0}}{2}+n(1+M\sqrt{L}r_{0})^{2}\mu\delta+nML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0})\right)\rho^{j}\left\lVert
x^{0}-x^{*}\right\rVert\\\
\leq&\frac{1+\rho}{n\mu}\left(n\frac{\tilde{L}r_{0}}{2}+n(1+M\sqrt{L}r_{0})^{2}\mu\delta+nML^{\frac{3}{2}}r_{0}(2+M\sqrt{L}r_{0})\right)\rho^{j}\left\lVert
x^{0}-x^{*}\right\rVert\\\ \leq&\rho^{j+1}\left\lVert
x^{0}-x^{*}\right\rVert.\end{split}$ (37)

The last inequality follows from equation (34).

#### Upper Bound of
$\nu(\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}))$

To provide the upper bound of the given metric, our first step is to establish
that
$\tau(\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}))$
is well-defined, by showing that:

$\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1}\succeq\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}).$
(38)

To establish the above relation, we need two observations. The first relation
is that

$\omega_{jn+m}^{-1}B_{m+1}^{jn+m}\succeq\nabla^{2}f_{m+1}(z_{m+1}^{jn+m})\succeq\left(1+\frac{Mr_{jn+m+1}}{2}\right)^{-1}K^{jn+m},$

where the first inequality follows from the induction hypothesis and the
second inequality follows from Lemma 27.

In addition, using the deductions from Lemma 31, we have:

$r_{jn+m+1}\leq\frac{2\alpha_{j}}{M}.$ (39)

Now we consider two cases depending on $m$.

(a): $0\leq m<n-1$ Since all the $B_{i}$ were scaled by a factor
$(1+\alpha_{j})^{2}$ at the end of the cycle $j-1$, i.e., at $t=jn$, we have
$B_{m+1}^{jn+m}=(1+\alpha_{j})^{2}B_{m+1}^{(j-1)n+m+1}$. Furthermore, from the
induction hypothesis, we have:

$\omega_{((j-1)n+m+1)_{+}}^{-1}B_{m+1}^{((j-1)n+m+1)_{+}}\succeq\nabla^{2}f_{m+1}(z_{m+1}^{((j-1)n+m+1)_{+}}).$

where we define $(x)_{+}=\max(x,0)$.

Since $0\leq m<n-1$, we have $\omega_{(j-1)n+m+1}=1$. Also,
$z_{m+1}^{jn+m}=z_{m+1}^{((j-1)n+m+1)_{+}}$. Therefore,

$\displaystyle
B_{m+1}^{jn+m}=(1+\alpha_{j})^{2}B_{m+1}^{((j-1)n+m+1)_{+}}\succeq$
$\displaystyle(1+\alpha_{j})^{2}\nabla^{2}f_{m+1}(z_{m+1}^{jn+m})$
$\displaystyle\succeq$
$\displaystyle(1+\alpha_{j})^{2}\left(1+\frac{Mr_{jn+m+1}}{2}\right)^{-1}K^{jn+m+1}$
$\displaystyle\succeq$ $\displaystyle(1+\alpha_{j})K^{jn+m+1}.$

The second inequality follows from Lemma 27. The last inequality can be
inferred from equation (39).

(b). $m=n-1$

Since the current index $m+1$ was last updated at time $t=jn$, we have
$B_{m+1}^{jn+m}=B_{m+1}^{jn}$ and $z_{m+1}^{jn+m}=z_{m+1}^{jn}$. Further, the
induction hypothesis yields that
$\omega_{jn}^{-1}B_{m+1}^{jn}\succeq\nabla^{2}f_{m+1}(z_{m+1}^{jn})$. Also,
$\omega_{jn}=(1+\alpha_{j})^{2}$ by definition. Therefore,

$\displaystyle
B_{m+1}^{jn+m}=\omega_{jn}(\omega_{jn}^{-1}B_{m+1}^{jn})\succeq$
$\displaystyle(1+\alpha_{j})^{2}\nabla^{2}f_{m+1}(z_{m+1}^{jn+m})$
$\displaystyle\succeq$
$\displaystyle(1+\alpha_{j})^{2}(1+\frac{Mr_{jn+m+1}}{2})^{-1}K^{jn+m+1}$
$\displaystyle\succeq$ $\displaystyle(1+\alpha_{j})K^{jn+m+1}.$

To sum up, for both cases $0\leq m<n$ and $m=n-1$, we have established that
$B_{m+1}^{jn+m}\succeq(1+\alpha_{j})K^{jn+m+1}$.

Applying Lemma 27 to $B_{m+1}^{jn+m}\succeq(1+\alpha_{j})K^{jn+m+1}$, we have:

$\displaystyle
B_{m+1}^{jn+m}\succeq(1+\alpha_{j})K^{jn+m+1}\succeq\left(1+\frac{Mr_{jn+m+1}}{2}\right)K^{jn+m+1}\succeq\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}).$

For Lemma 4.3, applying Lemma 27 and Lemma A.3, we have:

$\displaystyle\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1}=$
$\displaystyle\mathrm{SR1}(B_{m+1}^{jn+m},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}),\bar{u}(B_{m+1}^{jn+m},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}))$
$\displaystyle\succeq$ $\displaystyle\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}).$

For Lemma B.4, applying Lemma 27 and Lemma 28, it follows that:

$\displaystyle\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1}=$
$\displaystyle\mathrm{SR}\textrm{-}k(B_{m+1}^{jn+m},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}),\bar{U}(B_{m+1}^{jn+m},\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}))$
$\displaystyle\succeq$ $\displaystyle\nabla^{2}f_{m+1}(z_{m+1}^{jn+m+1}).$

Now define the sequence $\\{y_{k}\\}$, for $k=0,\dots,j+1$, such that
$\\{y_{k}\\}=\\{x^{0},z_{m+1}^{m+1},\dots,z_{m+1}^{jn+m+1}\\}$. From the
induction hypothesis and the upper bound of (37), we can infer that
$\left\lVert y_{k}-x^{*}\right\rVert\leq\rho^{k}\left\lVert
x^{0}-x^{*}\right\rVert$, for all $k$. Since $\\{y_{k}\\}$ comes from the
application of SR1 updates (SR-$k$ updates), the sequence satisfies the
conditions (30) from Lemma 31, this implies that

$\displaystyle\nu(\omega_{jn+m+1}^{-1}B_{m+1}^{jn+m+1},\nabla^{2}f_{m+1}(y_{j+1})\leq
c^{j+1}\delta.$

Since $y_{j+1}=z_{m+1}^{jn+m+1}$, this completes the proof. ∎

###  B.4 Proof of Lemma 4.5

We present and prove a more generalized version of Lemma 4.5:

######  Lemma B.5.

The sequence of iterates generated by the LISR-$k$ method satisfies

$\left\lVert
x^{t+1}-x^{*}\right\rVert\leq\left(1-\frac{k}{d}\right)^{\left\lceil\frac{t+1}{n}\right\rceil}\frac{1}{n}\sum_{i=1}^{n}\left\lVert
x^{t+1-i}-x^{*}\right\rVert.$

###### Proof.

Define $c\coloneqq 1-k/d$.

We prove the Lemma for a generic iteration $t=jn+m+1$, for some $j\geq 0$ and
$0\leq m<n$. We restate a few observations derived in the proof of Lemma 4.3.
First, the upper bound of
$\Gamma^{jn+m}=\left\lVert(\sum_{i=1}^{n}B_{i}^{jn+m})^{-1}\right\rVert$ is:

$\Gamma^{jn+m}\leq\frac{1+\rho}{n\mu}.$ (40)

Next, we recall the upper bound of $\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert$ is:

$\displaystyle\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq\begin{cases}c^{j+1}\mu\delta,&i\in[m]\\\
(1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}),&i\in[n]\setminus[m]\end{cases}$

Both cases can be summarized with a common upper bound:

$\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\leq(1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j}),\quad
i\in[n].$ (41)

Finally, we can also establish that

$\begin{split}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert\leq&\rho^{j+1}\left\lVert
x^{0}-x^{*}\right\rVert\leq\rho^{j}\left\lVert x^{0}-x^{*}\right\rVert,\quad
i\in[m]\\\ \left\lVert z_{i}^{jn+m}-x^{*}\right\rVert\leq&\rho^{j}\left\lVert
x^{0}-x^{*}\right\rVert,\quad i\in[n]\setminus[m]\end{split}$ (42)

From Lemma 4.1, we have

$\displaystyle\left\lVert z_{m+1}^{jn+m+1}-x^{*}\right\rVert\leq$
$\displaystyle\Gamma^{jn+m}\frac{\tilde{L}}{2}\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert^{2}+\Gamma^{jn+m}\sum_{i=1}^{n}\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\Gamma^{jn+m}\sum_{i=1}^{n}\left(\frac{\tilde{L}}{2}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert+\left\lVert
B_{i}^{jn+m}-\nabla^{2}f_{i}(z_{i}^{jn+m})\right\rVert\right)\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$ $\displaystyle\leq$
$\displaystyle\frac{1+\rho}{n\mu}\left(\frac{\tilde{L}\rho^{j}r_{0}}{2}+(1+M\sqrt{L}r_{0}\rho^{j})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0}\rho^{j})\right)\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$ $\displaystyle\overset{\rho<1}{\leq}$
$\displaystyle\frac{1+\rho}{n\mu}\left(\frac{\tilde{L}\rho^{j}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}c^{j}\mu\delta+ML^{\frac{3}{2}}r_{0}\rho^{j}(2+M\sqrt{L}r_{0})\right)\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$ $\displaystyle\leq$ $\displaystyle
c^{j}(1+\rho)\left(\frac{\tilde{L}\mu^{-1}r_{0}}{2}+(1+M\sqrt{L}r_{0})^{2}\delta+ML^{\frac{3}{2}}\mu^{-1}r_{0}(2+M\sqrt{L}r_{0})\right)\frac{1}{n}\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert$
$\displaystyle\overset{(\ref{rho_condition})}{\leq}$ $\displaystyle
c^{j}\rho(\frac{1}{n}\sum_{i=1}^{n}\left\lVert
z_{i}^{jn+m}-x^{*}\right\rVert)$ $\displaystyle\leq$ $\displaystyle
c^{j+1}(\frac{1}{n}\sum_{i=1}^{n}\left\lVert x^{t-i}-x^{*}\right\rVert).$

The third inequality is due to equation (40), (41) and (42). The fifth
inequality follows from $\rho<1-d^{-1}$ in the condition of Lemma 4.3 for
Lemma 4.5 and $\rho<1-kd^{-1}$ in the condition of Lemma B.4 for Lemma B.5. ∎

###  B.5 Proof of Theorem 18 and Theorem 22

###### Proof.

The proof is the same as Theorem 1 by (Lahoti et al. 2023) with $c\coloneqq\mu
d^{-1}L^{-1}$ replaced by $c\coloneqq d^{-1}$ for Theorem 18 and $c\coloneqq
kd^{-1}$ for Theorem 22. ∎

Algorithm 2 The LISR-1 method

Input: Initial point $x^{0}$, a predefined sequence
$\\{\alpha_{k}\\}_{k=0}^{\infty}$, and a sequence of initial matrices
$\\{E_{i}^{0}\\}_{i=1}^{n}$ s.t. $E_{i}^{0}\succeq\nabla^{2}f_{i}(z_{i}^{0})$,
$\forall i\in[n]$.
Output: The last iterate $x^{t}$.

1: Initialize $z_{i}^{0}=x^{0}$, $\forall i\in[n]$; Initialize
$B_{i}^{0}=(1+\alpha_{0})^{2}E_{i}^{0}$, $\forall i\in[n]$;

2: Maintain running auxiliary variables $x$, $\bar{u}$.

3: Let $t=0$.

4: while not converged do

5: Current index to be updated is $i_{t}\leftarrow t\mod n$;

6: Update $x^{t}$ as $x\leftarrow(\bar{B})^{-1}(\phi-g)$.

7: Compute $\omega_{t}$;// $\omega_{t}=(1+\alpha_{\left\lceil
t/n\right\rceil})^{2}$ if $t\mod n=0$ and 1 otherwise

8: Update $\bar{u}$ as
$\bar{u}\leftarrow\operatorname*{arg\,max}_{u\in\\{e_{i}\\}_{i=1}^{d}}(\langle
u,(1+\alpha_{\left\lceil t/n\right\rceil})^{2}B_{i_{t}}u\rangle-\langle
u,\nabla^{2}f_{i_{t}}(x^{t})u\rangle)$;

9: Update $Q\leftarrow\mathrm{SR}1((1+\alpha_{\left\lceil
t/n\right\rceil})^{2}B_{i_{t}},\nabla^{2}f_{i_{t}}(x^{t}),\bar{u})$;

10: Update $B^{\textrm{old}}\leftarrow B_{i_{t}}$;

11: Update $B_{i_{t}}\leftarrow Q$;

12: Update $z_{i_{t}}^{t}$ as $z_{i_{t}}^{t}=x^{t}$.

13: Update $\phi$ as
$\phi\leftarrow\omega_{t}(\phi-B^{\textrm{old}}z_{i_{t}})+B_{i_{t}}x^{t}$;

14: Update $g$ as $g\leftarrow g+(\nabla f_{i_{t}}(x^{t})-\nabla
f_{i_{t}}(z_{i_{t}}))$;

15: Update $(\bar{B})^{-1}$ as
$(\bar{B})^{-1}\leftarrow\omega_{t}^{-1}\mathrm{Sherman}\textrm{-}\mathrm{Morrison}(\bar{B}^{-1},v,v,\langle\bar{u},v\rangle)$;
// $v=((1+\alpha_{\left\lceil
t/n\right\rceil})^{2}B_{i_{t}}-\nabla^{2}f_{i_{t}}(z_{i_{t}}))\bar{u}$.

16: Increment the iteration counter $t$.

17: end while

18: return $x^{t}$.

Algorithm 3 The LISR-$k$ method

Input: Initial point $x^{0}$, a predefined sequence
$\\{\alpha_{k}\\}_{k=0}^{\infty}$, and a sequence of initial matrices
$\\{E_{i}^{0}\\}_{i=1}^{n}$ s.t. $E_{i}^{0}\succeq\nabla^{2}f_{i}(z_{i}^{0})$,
$\forall i\in[n]$.
Output: The last iterate $x^{t}$.

1: Initialize $z_{i}^{0}=x^{0}$, $\forall i\in[n]$; Initialize
$B_{i}^{0}=(1+\alpha_{0})^{2}E_{i}^{0}$, $\forall i\in[n]$;

2: Maintain running auxiliary variables $x$, $\bar{u}$.

3: Let $t=0$.

4: while not converged do

5: Current index to be updated is $i_{t}\leftarrow t\mod n$;

6: Update $x^{t}$ as $x\leftarrow(\bar{B})^{-1}(\phi-g)$.

7: Compute $\omega_{t}$;// $\omega_{t}=(1+\alpha_{\left\lceil
t/n\right\rceil})^{2}$ if $t\mod n=0$ and 1 otherwise

8: Update $\bar{U}$ as $\bar{U}\leftarrow E_{k}((1+\alpha_{\left\lceil
t/n\right\rceil})^{2}B_{i_{t}}-\nabla^{2}f_{i_{t}}(x^{t}))$;

9: Update $Q\leftarrow\mathrm{SR}\textrm{-}k((1+\alpha_{\left\lceil
t/n\right\rceil})^{2}B_{i_{t}},\nabla^{2}f_{i_{t}}(x^{t}),\bar{U})$;

10: Update $B^{\textrm{old}}\leftarrow B_{i_{t}}$;

11: Update $B_{i_{t}}\leftarrow Q$;

12: Update $z_{i_{t}}^{t}$ as $z_{i_{t}}^{t}=x^{t}$.

13: Update $\phi$ as
$\phi\leftarrow\omega_{t}(\phi-B^{\textrm{old}}z_{i_{t}})+B_{i_{t}}x^{t}$;

14: Update $g$ as $g\leftarrow g+(\nabla f_{i_{t}}(x^{t})-\nabla
f_{i_{t}}(z_{i_{t}}))$;

15: Update $(\bar{B})^{-1}$ as
$(\bar{B})^{-1}\leftarrow\omega_{t}^{-1}\mathrm{Sherman}\textrm{-}\mathrm{Morrison}(\bar{B}^{-1},V,V,\bar{U}^{\top}V)$;
// $V=((1+\alpha_{\left\lceil
t/n\right\rceil})^{2}B_{i_{t}}-\nabla^{2}f_{i_{t}}(z_{i_{t}}))\bar{U}$.

16: Increment the iteration counter $t$.

17: end while

18: return $x^{t}$.

Algorithm 4 Sherman-Morrison($A^{-1}$, U, V, W)

Input: $A^{-1}\in{\mathbb{R}}^{d\times d}$, $U\in{\mathbb{R}}^{d\times k}$,
$V\in{\mathbb{R}}^{d\times k}$, $W\in{\mathbb{R}}^{k\times k}$.
Output: The updated $A^{-1}$.

1: return $A^{-1}+A^{-1}U\left(W-U^{\top}A^{-1}V\right)^{-1}V^{\top}A^{-1}$.

##  Appendix C Efficient Implementation of the LISR-1/LISR-$k$ method

Recall that although Algorithm 1 enjoys per-iteration complexity of
${\mathcal{O}}(d^{2})$, it applies scaling to all Hessian approximations at
the iteration $t$ satisfying $t\mod n=0$. It results in the total FLOPs of
${\mathcal{O}}(nd^{2})$ at that iteration. We can improve the algorithm by
carrying out the scaling of $B_{i}$’s lazily. In particular, by only scaling
the Hessian approximations when they are used, we can improve the per-
iteration complexity of LISR-1 (LISR-$k$) to ${\mathcal{O}}(d^{2})$ for any
iteration. The resulting algorithm for the LISR-1 (LISR-$k$) is specified in
the Algorithm 2 (Algorithm 3). For Algorithm 2, it leverages the Sherman-
Morrison Formula (Algorithm 4 with $k=1$). On the other hand, Algorithm 3
applies the general Sherman-Morrison Formula (Algorithm 4 with $k>1$).

##  Appendix D Experiment Details

In this section, we present more details of the experiment. First, we give an
overview of the datasets used in the logistic regression task. Then we study
the impact of different hyperparameters on the quadratic function minimization
task. Finally, we present the comparison result of the LISR-$k$ methods with
different choices of $k$.

###  D.1 Datasets

We present more details of real-world datasets used in Table 2. Specifically,
each line shows the name of the dataset, the number of data instances, the
dimension of each data, the regularization parameter we chose for the logistic
regression objective function (24), and the resulting conditional number. All
datasets can be found on the LIBSVM website. We preprocess the datasets by
changing $\\{0,1\\}$ labels to $\\{\pm 1\\}$.

Table 2: Details of the real-world datasets. Dataset | $n$ | $d$ | $\lambda$ | $\kappa$
---|---|---|---|---
a9a | 32,561 | 123 | $1\times 10^{-3}$ | $4.71\times 10^{7}$
w8a | 49,749 | 300 | $1\times 10^{-4}$ | $2.93\times 10^{8}$
ijcnn1 | 49,749 |  22 | $1\times 10^{-4}$ | $1.37\times 10^{7}$
mushrooms |  8,124 | 112 | $1\times 10^{-3}$ | $2.08\times 10^{7}$
phishing | 11,055 |  68 | $1\times 10^{-4}$ | $1.79\times 10^{7}$
svmguide3 |  1,243 |  21 | $1\times 10^{-3}$ | $6.67\times 10^{5}$
german.numer |  1,000 |  24 | $1\times 10^{-3}$ | $2.11\times 10^{6}$
splice |  1,000 |  60 | $1\times 10^{-4}$ | $4.34\times 10^{6}$
covtype | 50,000 |  54 | $1\times 10^{-3}$ | $3.93\times 10^{7}$

###  D.2 Effect of different hyperparameters on the Quadratic Function
Minimization Task

|  |
---|---|---
(a) $n=10$, $d=100$ | (b) $n=100$, $d=100$ | (c) $n=1000$, $d=100$
|  |
(e) $n=100$, $d=80$ | (f) $n=100$, $d=200$ | (g) $n=100$, $d=500$
Figure 3: Comparison of the proposed methods with baselines for the quadratic
function minimization problem.

We present the comparison result of our proposed method and baseline methods
on the quadratic function minimization task in Figure 3. It can be observed
that tuning different choices of hyperparameters $n$ and $d$ does not have
much impact on the superlinear convergence of the LISR-1 and LISR-$k$ methods.
However, the time that superlinear convergence kicks in is correlated with the
dimension of the underlying problem. Interestingly, the normalized error of
the IQN method starts to diverge after reaching certain accuracy.

###  D.3 Effect of $k$ on the Logistic Regression Task

|  |
---|---|---
(a) a9a | (b) w8a | (c) ijcnn1
|  |
(e) mushrooms | (f) phishing | (g) svmguide3
|  |
(h) german.numer | (i) splice | (j) covtype
Figure 4: Comparison of the LISR-$k$ method with different choices of $k$ for the general function minimization.  |  |
---|---|---
(a) a9a | (b) w8a | (c) ijcnn1
|  |
(e) mushrooms | (f) phishing | (g) svmguide3
|  |
(h) german.numer | (i) splice | (j) covtype
Figure 5: Comparison of the LISR-$k$ method with different choices of $k$ for
the general function minimization.

In this subsection, we present the comparison result of the LISR-$k$ method
with different choices of $k$ on real-world datasets shown in Table 2. The
comparison result of normalized loss against number of the effective passes
can be viewed in Figure 4. It is not surprising to see that the LISR-$k$
generally converges at a faster rate if $k$ is set larger.

Additionally, we plot the comparison result of the normalized loss against
seconds in Figure 5. It can be seen that the LISR-$k$ method converges fastest
if $k$ is set to 10 or 20. Although the per-iteration cost of LISR-$k$ is
larger, the overall time to converge is still shorter than LISR-1 thanks to
the faster convergence rate.