# Sharpened Lazy Incremental Quasi-Newton Method

Aakash Lahoti∗
Carnegie Mellon University
<EMAIL_ADDRESS>&Spandan Senapati∗
University of Southern California
<EMAIL_ADDRESS>Ketan Rajawat
Indian Institute of Technology Kanpur
<EMAIL_ADDRESS>&Alec Koppel
J.P. Morgan AI Research
<EMAIL_ADDRESS>
###### Abstract

The problem of minimizing the sum of $n$ functions in $d$ dimensions is
ubiquitous in machine learning and statistics. In many applications where the
number of observations $n$ is large, it is necessary to use incremental or
stochastic methods, as their per-iteration cost is independent of $n$. Of
these, Quasi-Newton (QN) methods strike a balance between the per-iteration
cost and the convergence rate. Specifically, they exhibit a superlinear rate
with $\mathcal{O}(d^{2})$ cost in contrast to the linear rate of first-order
methods with $\mathcal{O}(d)$ cost and the quadratic rate of second-order
methods with $\mathcal{O}(d^{3})$ cost. However, existing incremental methods
have notable shortcomings: Incremental Quasi-Newton (IQN) only exhibits
asymptotic superlinear convergence. In contrast, Incremental Greedy BFGS (IGS)
offers explicit superlinear convergence but suffers from poor empirical
performance and has a per-iteration cost of $\mathcal{O}(d^{3})$. To address
these issues, we introduce the Sharpened Lazy Incremental Quasi-Newton Method
(SLIQN) that achieves the best of both worlds: an explicit superlinear
convergence rate, and superior empirical performance at a per-iteration
$\mathcal{O}(d^{2})$ cost. SLIQN features two key changes: first, it
incorporates a hybrid strategy of using both classic and greedy BFGS updates,
allowing it to empirically outperform both IQN and IGS. Second, it employs a
clever constant multiplicative factor along with a lazy propagation strategy,
which enables it to have a cost of $\mathcal{O}(d^{2})$. Additionally, our
experiments demonstrate the superiority of SLIQN over other incremental and
stochastic Quasi-Newton variants and establish its competitiveness with
second-order incremental methods.

###  1 INTRODUCTION

**footnotetext: Equal Contribution. Work done when AL and SS were at the
Indian Institute of Technology Kanpur.11footnotetext: Accepted at the 27th
International Conference on Artificial Intelligence and Statistics (AISTATS)
2024.

We consider the finite sum minimization problem,

$\displaystyle{\bm{x}}^{\star}=\arg\min_{{\bm{x}}\in{\mathbb{R}}^{d}}\frac{1}{n}\sum_{i=1}^{n}f_{i}({\bm{x}}),$
($\mathcal{P}$)

where each $f_{i}$ is $\mu$-strongly convex, $L$-smooth and has a Lipschitz
continuous Hessian. The canonical example of ($\mathcal{P}$) is the empirical
risk minimization problem in supervised learning, where ${\bm{x}}$ denotes
model parameters, $n$ is the number of training samples, and $f_{i}$ denotes
the loss incurred by the $i^{\text{th}}$ sample. Other instances of this
problem arise in maximum likelihood estimation (Li et al., (2021, 2020)),
control theory (Wu et al., (2018)), and unsupervised learning (Song and Ermon,
(2020)). In many applications, ($\mathcal{P}$) is both high-dimensional (large
$d$) and data-intensive (large $n$).

When $n$ is large, it becomes infeasible to process the entire dataset at
every iteration, thus making classical algorithms such as gradient descent or
Newton’s method impractical. Consequently, stochastic and incremental variants
of these algorithms have been widely adopted for such problems, because their
per-iteration complexity is independent of $n$. While first-order methods like
stochastic gradient descent (SGD) enjoy a low per-iteration complexity of
$\mathcal{O}(d)$, their convergence rates, even with enhancements like
variance reduction or acceleration (Defazio et al., (2014); Johnson and Zhang,
(2013)), remain linear at best. In contrast, second-order methods like Newton
Incremental Method (NIM) by Rodomanov and Kropotov, (2016), achieve a
superlinear rate but at a ${\mathcal{O}}(d^{3})$ cost, which is prohibitively
large for high-dimensional problem settings.

Stochastic and incremental Quasi-Newton (QN) methods strike a balance between
the computational efficiency of SGD and the fast convergence rate of NIM.
Specifically, the Incremental Quasi-Newton (IQN) method from Mokhtari et al.,
(2018), was the first QN method to achieve a superlinear convergence rate with
a per-iteration complexity of ${\mathcal{O}}(d^{2})$. However, the analysis
presented in Mokhtari et al., (2018) was asymptotic and did not include an
explicit rate of convergence. Typically, explicit rates are preferred over
asymptotic ones as they enable a more fine-grained comparison among
algorithms. For instance, with $\rho\in(0,1)$, both
$\mathcal{O}(\rho^{t^{2}})$ and $\mathcal{O}(\rho^{t\ln(t)})$ qualify as
superlinear rates. However, the former, $\mathcal{O}(\rho^{t^{2}})$, is faster
than the latter, $\mathcal{O}(\rho^{t\ln(t)})$. Furthermore, the mathematical
expression of $\rho$ informs the rate’s dependence on parameters like
condition number and dimension.

The Incremental Greedy BFGS (IGS) method by Gao et al., (2020) aimed to
address this issue by incorporating the greedy updates, first introduced in
Rodomanov and Nesterov, 2021a , into the IQN framework. While IGS achieves an
explicit superlinear convergence rate of ${\mathcal{O}}(e^{-t^{2}/n^{2}})$, it
suffers from several major drawbacks: First, like NIM, it has a large per-
iteration cost of $\mathcal{O}(d^{3})$. This stems from IGS’s Hessian updates
not being low ran, which precludes an efficient evaluation of the Hessian
inverse. Since this complexity mirrors that of NIM, it undermines the
computational advantages of incorporating QN updates. Second, the empirical
performance of IGS was not studied in Gao et al., (2020), and our experiments
(ref. Figure 2) indicate that on multiple datasets, IGS severely underperforms
compared to IQN. Finally, IGS lacked theoretical analysis to support its
lemmas and theorems. In this light, we ask the following question:

Can we devise an incremental QN method with a per-iteration complexity of
${\mathcal{O}}(d^{2})$, achieving the best-known incremental rate of
${\mathcal{O}}(e^{-t^{2}/n^{2}})$, and demonstrating superior empirical
performance?

We put forth the Sharpened Lazy IQN (SLIQN) method that meets all these
objectives. SLIQN is inspired by the recent work of Jin et al., (2022), which
showcased the superior performance of sharpened updates over greedy updates in
the non-incremental setting. We first show that a direct incorporation of
sharpened updates into the IQN framework does not work because the Hessian
update matrices corresponding to sharpened updates are not low rank, and
therefore the resulting Sharpened IQN (SIQN) method incurs a per-iteration
cost of ${\mathcal{O}}(d^{3})$. We then propose our novel Sharpened Lazy IQN
(SLIQN) algorithm that overcomes this limitation by modifying the updates of
SIQN using a clever constant multiplicative factor and incorporating a lazy
propagation strategy. The resulting algorithm incurs a per-iteration
complexity of ${\mathcal{O}}(d^{2})$ and achieves a convergence rate of
${\mathcal{O}}(\rho^{t^{2}/n^{2}})$, where $\rho:=1-\mu/dL$, which is the
best-known rate in the incremental setting. We also establish an explicit
linear rate of convergence of the Hessian approximation to the true Hessian.
Moreover, in contrast to IGS, we provide a comprehensive theoretical analysis.
Furthermore, we demonstrate the superior empirical performance of SLIQN as
compared to IQN, IGS, and other state-of-the-art incremental and stochastic QN
methods.111The code for the experiments is available on the repository:
https://github.com/aakashlahoti/sliqn. Notably, SLIQN demonstrates performance
competitive to NIM, which is a second-order algorithm, with a
${\mathcal{O}}(d^{3})$ cost, that utilizes the full Hessian information when
taking the descent step.

###  2 RELATED WORK

In recent decades, several works have developed first-order, QN, and second-
order methods for stochastic or incremental settings. Typically, the goal for
first-order methods is to achieve a linear rate at a ${\mathcal{O}}(d)$ cost,
while for QN and second-order methods, the goals are to achieve superlinear
rates at costs of ${\mathcal{O}}(d^{2})$ and ${\mathcal{O}}(d^{3})$,
respectively. These methods cater to different objectives: first-order methods
are preferred for low-precision solutions, due to their lower computational
cost, whereas higher-order methods are more effective for high-precision
solutions, due to their faster rate.

Early works like Mokhtari and Ribeiro, (2014, 2015); Byrd et al., (2016) were
only successful in developing QN methods with sub-linear convergence
guarantees. Subsequent works like Moritz et al., (2016); Chang et al., (2019);
Derezinski, (2023) employed various acceleration and variance reduction
techniques to recover a linear rate. IQN by Mokhtari et al., (2018) was the
first QN algorithm to achieve an asymptotic superlinear rate of convergence.
IGS by Gao et al., (2020) employed greedy updates, introduced in Rodomanov and
Nesterov, 2021a , within the IQN framework to derive a explicit superlinear
rate, albeit at a large $\mathcal{O}(d^{3})$ per-iteration cost. Another
recent work, Chen et al., (2022), put forth a QN style algorithm with a cost
of $\mathcal{O}(d)$ for Generalized Linear Models (GLMs). However, the method
only enjoys a linear rate of convergence, and is inefficient for general
functions with a cost of $\mathcal{O}(d^{3})$.

Other works have focused on developing first-order and second-order methods
for stochastic or incremental settings. First order methods like Defazio et
al., (2014); Johnson and Zhang, (2013) have employed variance reduction
techniques to derive methods with a linear rate of convergence. On the second-
order front, recent works include the Newton Incremental method (NIM) by
Rodomanov and Kropotov, (2016) and Stochastic Newton (SN) by Kovalev et al.,
(2019). While SN and NIM are both Newton-like methods with a per-iteration
complexity of $\mathcal{O}(d^{3})$, NIM has a fast superlinear convergence
rate while SN only enjoys a linear rate of convergence. The only setting under
which SN has been shown to converge superlinearly is with a full batch of size
$n$. We have consolidated the memory usage, computational cost, and
convergence rates of principal-related algorithms in Table 1.

In a different line of work, Newton-LESS from Derezinski et al., (2021)
utilized sketching algorithms like Newton Sketch Pilanci and Wainwright,
(2016, 2017) to attain a local linear convergence rate. Additional works such
as Gonen et al., (2016); Liu et al., (2019) have used second-order information
to accelerate SVRG, which is a first-order method. However, both methods were
only able to attain an improved linear rate of convergence.

Algorithm | Memory | Computation cost | Convergence Rate | Limit
---|---|---|---|---
SN | $\bm{\mathcal{O}(d^{2})}$ | $\mathcal{O}(d^{3})$ | Linear | Non-asymptotic
NIM | $\mathcal{O}(nd+d^{2})$ | $\mathcal{O}(d^{3})$ | Superlinear | Non-asymptotic
IQN | $\mathcal{O}(nd^{2})$ | $\bm{\mathcal{O}(d^{2})}$ | Superlinear | Asymptotic
IGS | $\mathcal{O}(nd^{2})$ | $\mathcal{O}(d^{3})$ | Superlinear | Non-asymptotic
SIQN(This work) | $\mathcal{O}(nd^{2})$ | $\mathcal{O}(d^{3})$ | Superlinear | Non-asymptotic
SLIQN(This work) | $\mathcal{O}(nd^{2})$ | $\bm{\mathcal{O}(d^{2})}$ | Superlinear | Non-asymptotic
Table 1:  Comparison of the maximum memory requirement, computation cost (per
iteration), convergence rate, and the limit of attainment of the convergence
rate for different algorithms.

###  3 NOTATION AND PRELIMINARIES

Vectors (matrices) are denoted by lowercase (uppercase) bold alphabets. The
$i$-th standard basis vector of ${\mathbb{R}}^{d}$ is denoted by
${\bm{e}}_{i}\in{\mathbb{R}}^{d}$ for $i\in[d]\coloneqq\\{1,\ldots,d\\}$. We
define the index function $i_{t}\coloneqq 1+(t-1)\bmod n$. The symbol
${\mathbf{0}}$ denotes the all-zero matrix or vector, whose size can be
inferred from the context. We use ${\bm{X}}\succeq{\mathbf{0}}$ and
${\bm{X}}\succ{\mathbf{0}}$ to denote that the symmetric matrix ${\bm{X}}$ is
positive semi-definite and positive definite, respectively. Likewise, the
notation ${\bm{X}}\succeq{\bm{Y}}$ (${\bm{X}}\succ{\bm{Y}}$) denotes
${\bm{X}}-{\bm{Y}}\succeq{\mathbf{0}}$ (${\bm{X}}-{\bm{Y}}\succ{\mathbf{0}}$).
For vectors ${\bm{u}},{\bm{v}}\in{\mathbb{R}}^{d}$, we denote the inner
product by
$\langle{\bm{u}},{\bm{v}}\rangle\coloneqq{\bm{u}}^{{\mathsf{T}}}{\bm{v}}$ and
the Euclidean norm by
$\left\|{\bm{u}}\right\|\coloneqq\sqrt{\langle{\bm{u}},{\bm{u}}\rangle}$. For
matrices ${\bm{X}},{\bm{Y}}\in{\mathbb{R}}^{d\times d}$, we define
$\langle{\bm{X}},{\bm{Y}}\rangle\coloneqq\mathrm{Tr}({\bm{X}}^{\mathsf{T}}{\bm{Y}})$,
and we let $\left\|{\bm{X}}\right\|$ denote the spectral norm of the matrix.
Given a convex function $f:{\mathbb{R}}^{d}\to{\mathbb{R}}$, we define the
norm of the vector ${\bm{y}}$ with respect to $\nabla^{2}f({\bm{x}})$ as
$\left\|{\bm{y}}\right\|_{\bm{x}}\coloneqq\sqrt{\langle{\bm{y}},\nabla^{2}f({\bm{x}}){\bm{y}}\rangle}$.
For a function $f$, we denote $\mu$ as the strong convexity parameter, $L$ as
the smoothness parameter, $\tilde{L}$ as the Hessian Lipschitz continuity
parameter, and $M=\tilde{L}\mu^{-\frac{3}{2}}$ as the strong self-concordance
parameter.

####  3.1 Quasi-Newton (QN) Methods

We introduce QN methods as an iterative algorithm to optimize the problem
($\mathcal{P}$) with $n=1$. At iteration $t\in\mathbb{Z}_{+}$, given the
current iterate ${\bm{x}}^{t}$ and the positive definite Hessian approximation
${\bm{B}}^{t}$ of $\nabla^{2}f({\bm{x}}^{t})$, the next iterate
${\bm{x}}^{t+1}$ is computed as,

$\displaystyle{\bm{x}}^{t+1}$
$\displaystyle={\bm{x}}^{t}-({\bm{B}}^{t})^{-1}\nabla f({\bm{x}}^{t}).$ (1)

The Hessian approximation for the next iteration ${\bm{B}}^{t+1}$, is obtained
by applying a constant rank update to ${\bm{B}}^{t}$. The precise update
distinguishes the exact type of QN algorithm, such as BFGS, DFP, or SR1
(Nocedal and Wright, (1999)). The efficiency of QN methods,
$\mathcal{O}(d^{2})$, over second-order methods, $\mathcal{O}(d^{3})$, stems
from the fact that the ${\bm{B}}^{t}$ update is low rank, which allows us to
use Sherman-Morrison formula (Appendix A) to efficiently evaluate
$({\bm{B}}^{t+1})^{-1}$ from $({\bm{B}}^{t})^{-1}$ in $\mathcal{O}(d^{2})$
cost.

Though in the remainder of the paper, we are primarily concerned with BFGS
updates, all the follows can also be extended to the entire restricted Broyden
class (Appendix H). Given a matrix ${\bm{K}}$ and its approximation
${\bm{B}}$, the generalized BFGS update refines this approximation along
direction ${\bm{u}}\in{\mathbb{R}}^{d}$ as,

$\displaystyle\mathrm{BFGS}({\bm{B}},{\bm{K}},{\bm{u}})={\bm{B}}_{+}\coloneqq{\bm{B}}-\frac{{\bm{B}}{\bm{u}}{\bm{u}}^{{\mathsf{T}}}{\bm{B}}}{\langle{\bm{u}},{\bm{B}}{\bm{u}}\rangle}+\frac{{\bm{K}}{\bm{u}}{\bm{u}}^{{\mathsf{T}}}{\bm{K}}}{\langle{\bm{u}},{\bm{K}}{\bm{u}}\rangle}.$
(2)

Setting
${\bm{K}}^{t}=\int_{0}^{1}\nabla^{2}f({\bm{x}}^{t}+\tau({\bm{x}}^{t+1}-{\bm{x}}^{t}))d\tau$,
and ${\bm{u}}^{t}={\bm{s}}^{t}:={\bm{x}}^{t+1}-{\bm{x}}^{t}$ yields the
classical BFGS update,

$\displaystyle{\bm{B}}^{t+1}$
$\displaystyle=\mathrm{BFGS}({\bm{B}}^{t},{\bm{K}}^{t},{\bm{u}}^{t})={\bm{B}}^{t}-\frac{{\bm{B}}^{t}{\bm{s}}^{t}({\bm{s}}^{t})^{{\mathsf{T}}}{\bm{B}}^{t}}{\langle{\bm{s}}^{t},{\bm{B}}^{t}{\bm{s}}^{t}\rangle}+\frac{{\bm{y}}^{t}({\bm{y}}^{t})^{{\mathsf{T}}}}{\langle{\bm{s}}^{t},{\bm{y}}^{t}\rangle},$
(3)

where ${\bm{y}}_{t}\coloneqq{\bm{K}}^{t}{\bm{s}}^{t}=\nabla
f({\bm{x}}^{t+1})-\nabla f({\bm{x}}^{t})$. This update seeks to approximate
the Hessian along the Newton direction ${\bm{s}}^{t}$ and has been shown by
Jin and Mokhtari, (2022); Rodomanov and Nesterov, 2021c ; Rodomanov and
Nesterov, 2021b  to achieve a superlinear convergence rate. Furthermore, since
BFGS makes a rank $2$ update to ${\bm{B}}^{t}$, we can use the Sherman-
Morrison formula twice to evaluate $({\bm{B}}^{t+1})^{-1}$ from
$({\bm{B}}^{t})^{-1}$ in $\mathcal{O}(d^{2})$ cost.

In contrast to classical BFGS update, the greedy BFGS update by Rodomanov and
Nesterov, 2021a  sets ${\bm{K}}^{t}=\nabla^{2}f({\bm{x}}^{t})$, and defines
the greedy vector,

$\displaystyle\bar{{\bm{u}}}({\bm{B}}^{t},{\bm{K}}^{t})\coloneqq\operatorname*{arg\,max}_{{\bm{u}}\in\\{{\bm{e}}_{i}\\}_{i=1}^{d}}\tfrac{\langle{\bm{u}},{\bm{B}}^{t}{\bm{u}}\rangle}{\langle{\bm{u}},{\bm{K}}^{t}{\bm{u}}\rangle},$
(4)

which results in
${\bm{B}}^{t+1}=\mathrm{BFGS}({\bm{B}}^{t},{\bm{K}}^{t},{\bar{{\bm{u}}}}^{t}({\bm{B}}^{t},{\bm{K}}^{t}))$.
Greedy BFGS, similar to classic BFGS, exhibits a superlinear rate of
convergence. However, unlike classic BFGS, it can guarantee convergence in the
$\sigma$ sense. Specifically, the Hessian approximation error,

$\displaystyle\sigma({\bm{B}}^{t},{\bm{K}}^{t})\coloneqq\langle({\bm{K}}^{t})^{-1},{\bm{B}}^{t}-{\bm{K}}^{t}\rangle,$
(5)

decays linearly with $t$. In practice, a trade-off exists between the two
updates. Greedy BFGS dedicates initial iterations to construct a reliable
Hessian approximation. In contrast, classic BFGS gains an initial advantage
because it updates along the Newton direction. Greedy BFGS only outperforms
classic BFGS if it has enough time to build an accurate approximation before
convergence. Sharpened BFGS proposed by Jin et al., (2022) incorporated both
classic and greedy BFGS updates, to achieve the best of both updates: an
explicit superlinear convergence of ${\bm{x}}^{t}$, a linear convergence of
${\bm{B}}^{t}$ and no initial “ramp-up” phase to build the approximation.

####  3.2 Incremental Quasi-Newton (IQN)

We now introduce IQN by Mokhtari et al., (2018). For each iteration $t\geq 1$,
IQN maintains tuples of the form $({\bm{z}}_{i}^{t},\nabla
f_{i}({\bm{z}}_{i}^{t}),{\bm{B}}_{i}^{t})$ for each index $i\in[n]$. Here
${\bm{z}}_{i}^{t}\in{\mathbb{R}}^{d}$ is the iterate corresponding to the
function $f_{i}$, $\nabla f_{i}({\bm{z}}_{i}^{t})$ is the gradient of $f_{i}$
at ${\bm{z}}_{i}^{t}$, and ${\bm{B}}_{i}^{t}$ is the positive definite Hessian
approximation of $\nabla^{2}f_{i}({\bm{z}}_{i}^{t})$.

IQN begins by constructing a second-order Taylor approximation $g_{i}^{t}$ of
each $f_{i}$, centered at ${\bm{z}}_{i}^{t-1}$ and using the Hessian
approximation ${\bm{B}}_{i}^{t-1}$ as,

$\displaystyle g_{i}^{t}({\bm{x}})$
$\displaystyle:=f_{i}({\bm{z}}_{i}^{t-1})+\langle\nabla{f_{i}({\bm{z}}_{i}^{t-1})},{\bm{x}}-{\bm{z}}_{i}^{t-1}\rangle+\frac{1}{2}\langle{\bm{x}}-{\bm{z}}_{i}^{t-1},{\bm{B}}_{i}^{t-1}({\bm{x}}-{\bm{z}}_{i}^{t-1})\rangle.$
(6)

The iterate ${\bm{x}}^{t}$ is then calculated as,

$\displaystyle{\bm{x}}^{t}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{x}}\frac{1}{n}\sum_{i=1}^{n}g_{i}^{t}({\bm{x}})=(\bar{{\bm{B}}}^{t-1})^{-1}\bigg{(}\sum_{i=1}^{n}{\bm{B}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}-\nabla
f_{i}({\bm{z}}_{i}^{t-1})\bigg{)},$ (7)

where $\bar{{\bm{B}}}^{t-1}=\sum_{i=1}^{n}{\bm{B}}_{i}^{t-1}$. In every
iteration $t$, IQN only updates the tuple whose index is given by
$i_{t}=1+(t-1)\bmod n$, using the following scheme:

  1. 1.

${\bm{z}}_{i_{t}}^{t}={\bm{x}}^{t}$, $\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})=\nabla f_{i_{t}}({\bm{x}}^{t})$.

  2. 2.

${\bm{z}}_{j}^{t}={\bm{z}}_{j}^{t-1}$, $\nabla f_{j}({\bm{z}}_{j}^{t})=\nabla
f_{j}({\bm{z}}_{j}^{t-1})$, and ${\bm{B}}_{j}^{t}={\bm{B}}_{j}^{t-1}$, for all
$j\neq i_{t}$.

  3. 3.

${\bm{B}}_{i_{t}}^{t}=\mathrm{BFGS}({\bm{B}}_{i_{t}}^{t-1},{\bm{K}}^{t},{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}),$
where
${\bm{K}}^{t}=\int_{0}^{1}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1}+\tau({\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}))d\tau$.

The per-iteration complexity of IQN: the cost of gradient evaluation in (1) is
$\mathcal{O}(d)$, the no-operation step in (2) incurs no cost, and the BFGS
step (3) has a $\mathcal{O}(d^{2})$ cost as it is a constant rank update to
${\bar{{\bm{B}}}}_{i_{t}}^{t}$. To compute the iterate ${\bm{x}}^{t+1}$, we
first evaluate the inverse of
${\bar{{\bm{B}}}}^{t}={\bar{{\bm{B}}}}^{t-1}+{\bm{B}}_{i_{t}}^{t}-{\bm{B}}_{i_{t}}^{t-1}$
from $(\bar{{\bm{B}}}^{t-1})^{-1}$ using the Sherman-Morrison formula at a
cost of ${\mathcal{O}}(d^{2})$. Then, we calculate
$\sum_{i=1}^{n}{\bm{B}}_{i}^{t}{\bm{z}}_{i}^{t}-\nabla
f_{i}({\bm{z}}_{i}^{t})$ from the memoized value of
$\sum_{i=1}^{n}{\bm{B}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}-\nabla
f_{i}({\bm{z}}_{i}^{t-1})$ in ${\mathcal{O}}(d^{2})$ cost. Therefore, the
overall per-iteration cost of IQN is $\mathcal{O}(d^{2})$. Please refer to
Appendix C for details.

###  4 PROPOSED ALGORITHM

We first introduce the SIQN method, which incorporates the sharpened updates
into the IQN framework. However, this direct incorporation results in an
inefficient ${\mathcal{O}}(d^{3})$ method. This is because, ensuring the
positive semi-definiteness of the Hessian approximation with respect to the
true Hessian results in Hessian updates which are not low rank. And
consequently, the Hessian inversion incurs a large $\mathcal{O}(d^{3})$ cost.
We then propose the SLIQN algorithm that overcomes this problem by modifying
the SIQN updates using a constant corrective multiplicative factor based on
the theoretical analysis of SIQN. It then utilizes a novel lazy propagation
strategy to implement this factor correction efficiently with a per-iteration
cost of $\mathcal{O}(d^{2})$.

####  4.1 Sharpened Incremental Quasi-Newton (SIQN)

Similar to IQN, SIQN also maintains tuples of the form
$({\bm{z}}_{i}^{t},\nabla f_{i}({\bm{z}}_{i}^{t}),{\bm{B}}_{i}^{t})$ for each
iteration $t\geq 1$ and each index $i\in[n]$. The iterate ${\bm{x}}^{t}$ is
calculated as,

$\displaystyle{\bm{x}}^{t}=(\bar{{\bm{B}}}^{t-1})^{-1}\bigg{(}\sum_{i=1}^{n}{\bm{B}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}-\nabla
f_{i}({\bm{z}}_{i}^{t-1})\bigg{)},$ (8)

where $\bar{{\bm{B}}}^{t-1}=\sum_{i=1}^{n}{\bm{B}}_{i}^{t-1}$. The tuples are
updated in a deterministic cyclic order as follows:

  1. 1.

${\bm{z}}_{i_{t}}^{t}={\bm{x}}^{t}$, $\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})=\nabla f_{i_{t}}({\bm{x}}^{t})$.

  2. 2.

${\bm{z}}_{j}^{t}={\bm{z}}_{j}^{t-1}$, $\nabla f_{j}({\bm{z}}_{j}^{t})=\nabla
f_{j}({\bm{z}}_{j}^{t-1})$, and ${\bm{B}}_{j}^{t}={\bm{B}}_{j}^{t-1}$, for all
$j\neq i_{t}$.

To update ${\bm{B}}_{i_{t}}^{t}$, SIQN first performs the classic BFGS update
followed by the greedy update,

  1. 3.

${{\bm{Q}}^{t}}=\mathrm{BFGS}((1+\beta_{t})^{2}{\bm{B}}_{i_{t}}^{t-1},(1+\beta_{t}){\bm{K}}^{t},{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1})$,

${\bm{B}}_{i_{t}}^{t}=\mathrm{BFGS}({{\bm{Q}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}),\bar{{\bm{u}}}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))),$

where,
${\bm{K}}^{t}\coloneqq\int_{0}^{1}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1}+\tau({\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}))d\tau$
and the scaling factor
$\beta_{t}\coloneqq\tfrac{M}{2}\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}$
ensures that before the classical BFGS update, the Hessian approximation is
positive semi-definite with respect to the corresponding Hessian,
$(1+\beta_{t})^{2}{{\bm{B}}}_{i_{t}}^{t-1}\succeq(1+\beta_{t}){\bm{K}}^{t}$.
For technical details, please refer to Lemma E.1. The pseudo-code for SIQN is
provided in Algorithm 1.

Algorithm 1 Sharpened Incremental Quasi-Newton

1: Initialize:
$\\{{\bm{z}}_{i}^{0}={\bm{x}}^{0}\\}_{i=1}^{n},\\{{\bm{B}}_{i}^{0}\\}_{i=1}^{n}$
such that for all $i\in[n]$,
${\bm{B}}_{i}^{0}\succeq\nabla^{2}f_{i}({\bm{z}}_{i})^{0}$

2: while not converged:

3: Set current index $i_{t}=(t-1)\mod n+1$;

4: Update ${\bm{x}}^{t}$ as per (8);

5: Update ${\bm{z}}_{i_{t}}^{t}={\bm{x}}^{t}$ and ${\bm{B}}_{i_{t}}^{t}$ as
per (3);

6: Update the tuples with index $j\neq i_{t}$ as
${\bm{z}}_{j}^{t}={\bm{z}}_{j}^{t-1},{\bm{B}}_{j}^{t}={\bm{B}}_{j}^{t-1}$;

7: Increment the iteration counter $t$;

8: end while

We now consider the per-iteration complexity of SIQN: the cost of gradient
evaluation is $\mathcal{O}(d)$, the no-operation step incurs no cost and the
BFGS step can be computed in $\mathcal{O}(d^{2})$ cost. To compute
${\bm{x}}^{t+1}$, we need to calculate the inverse of
${\bar{{\bm{B}}}}^{t}={\bar{{\bm{B}}}}^{t-1}+{\bm{B}}_{i_{t}}^{t}-{\bm{B}}_{i_{t}}^{t-1}$,

$\displaystyle{\bm{B}}_{i_{t}}^{t}-{\bm{B}}_{i_{t}}^{t-1}$
$\displaystyle=\mathrm{BFGS}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}),\bar{{\bm{u}}}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})))-{\bm{B}}_{i_{t}}^{t-1},$
$\displaystyle={\bm{Q}}^{t}-\frac{{\bm{Q}}^{t}\bar{{\bm{u}}}^{t}(\bar{{\bm{u}}}^{t})^{{\mathsf{T}}}{\bm{Q}}^{t}}{\langle\bar{{\bm{u}}}^{t},{\bm{Q}}^{t}\bar{{\bm{u}}}^{t}\rangle}+\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})\bar{{\bm{u}}}^{t}(\bar{{\bm{u}}}^{t})^{{\mathsf{T}}}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})}{\langle\bar{{\bm{u}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})\bar{{\bm{u}}}^{t}\rangle}-{\bm{B}}_{i_{t}}^{t-1},$
$\displaystyle=((1+\beta_{t})^{2}-1){\bm{B}}_{i_{t}}^{t-1}-(1+\beta_{t})^{2}\frac{{\bm{B}}_{i_{t}}^{t-1}{\bm{s}}_{i_{t}}^{t}{{\bm{s}}_{i_{t}}^{t}}^{{\mathsf{T}}}{\bm{B}}_{i_{t}}^{t-1}}{\langle{\bm{s}}_{i_{t}}^{t},{\bm{B}}_{i_{t}}^{t-1}{\bm{s}}_{i_{t}}^{t}\rangle}$
$\displaystyle\hskip
14.22636pt+(1+\beta_{t})\frac{{\bm{y}}_{i_{t}}^{t}{{\bm{y}}_{i_{t}}^{t}}^{{\mathsf{T}}}}{\langle{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}},{\bm{s}}_{i_{t}}^{t}\rangle}-\frac{{\bm{Q}}^{t}\bar{{\bm{u}}}^{t}(\bar{{\bm{u}}}^{t})^{{\mathsf{T}}}{\bm{Q}}^{t}}{\langle\bar{{\bm{u}}}^{t},{\bm{Q}}^{t}\bar{{\bm{u}}}^{t}\rangle}+\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})\bar{{\bm{u}}}^{t}(\bar{{\bm{u}}}^{t})^{{\mathsf{T}}}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})}{\langle\bar{{\bm{u}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})\bar{{\bm{u}}}^{t}\rangle},$
(9)

where
$\bar{{\bm{u}}}^{t}=\bar{{\bm{u}}}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$,
${\bm{s}}_{i_{t}}^{t}={\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}$ and
${\bm{y}}_{i_{t}}^{t}=\nabla f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})$. Observe that the expression (9) is
generally not a matrix of constant rank, since the rank of
$((1+\beta_{t})^{2}-1){\bm{B}}_{i_{t}}^{t-1}$ may be as large as $d$.
Therefore, it is not possible to compute the inverse of ${\bar{{\bm{B}}}}^{t}$
by using the inverse of ${\bar{{\bm{B}}}}^{t-1}$ in $\mathcal{O}(d^{2})$ cost.

####  4.2 Sharpened Lazy Incremental Quasi-Newton (SLIQN)

We now present the SLIQN method as a solution to the aforementioned issues. We
observe that the primary reason for the update of $\bar{{\bm{B}}}^{t}$ in SIQN
not being low-rank, is the presence of the scaling factor $\beta_{t}$. To
resolve this issue, we begin by noting that the convergence analysis of SIQN
(ref. Appendix E) indicates that there exists a factor $\alpha_{\lceil
t/n-1\rceil}$, such that $\beta_{t}\leq\alpha_{\lceil t/n-1\rceil}$ and using
$\alpha_{\lceil t/n-1\rceil}$ instead of $\beta_{t}$ preserves the convergence
properties of SIQN. Since $\alpha_{\lceil t/n-1\rceil}$ is constant throughout
an epoch, rather than multiplying $\alpha_{\lceil t/n-1\rceil}$ to
${\bm{B}}_{i_{t}}^{t-1}$ at iteration $t$, we instead pre-multiply
$\alpha_{\lceil t/n-1\rceil}$ to each Hessian approximation at the start of
every epoch. This enables us to compute the inverse of $\bar{{\bm{B}}}^{t}$
trivially by dividing the old inverse by $\alpha_{\lceil t/n-1\rceil}$.
However, this pre-multiplication step is a $\mathcal{O}(nd^{2})$ operation and
it undermines the utility of incremental algorithms. To address this issue, we
employ a lazy propagation strategy, wherein we scale the individual Hessian
approximations just before they are updated in their respective iterations,
but treat all memoized quantities as if the approximations are already scaled.
These key changes enable SLIQN to achieve an $\mathcal{O}(d^{2})$ per-
iteration cost along with an explicit superlinear rate.

In what follows, we will denote the Hessian approximations by
$\\{{\bm{D}}_{i}^{t}\\}_{i=1}^{n}$ instead of
$\\{{\bm{B}}_{i}^{t}\\}_{i=1}^{n}$ to distinguish SLIQN updates from SIQN
updates. We now formally present the SLIQN algorithm.

Initialization: At $t=0$, we initialize the iterates
$\\{{\bm{z}}_{i}^{0}\\}_{i=1}^{n}$ as ${\bm{z}}_{i}^{0}={\bm{x}}^{0}$ and
their corresponding hessian approximations $\\{{\bm{D}}_{i}^{0}\\}_{i=1}^{n}$
as ${\bm{D}}_{i}^{0}=(1+\alpha_{0})^{2}{\bm{I}}_{i}^{0}$, where
${\bm{x}}_{0},\alpha_{0},\\{{\bm{I}}_{i}^{0}\\}_{i=1}^{n}$ satisfy the the
premise of Lemma 2.

Iterative Updates: For each iteration $t\geq 1$, we set the iterate
${\bm{x}}^{t}$ as,

$\displaystyle{\bm{x}}^{t}=\left(\bar{{\bm{D}}}^{t-1}\right)^{-1}\bigg{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}-\nabla
f_{i}({\bm{z}}_{i}^{t-1})\bigg{)},$ (10)

where $\bar{{\bm{D}}}^{t-1}=\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}$. The scheme to
update each tuple is as follows:

  1. 1.

${\bm{z}}_{i_{t}}^{t}={\bm{x}}^{t}$, $\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})=\nabla f_{i_{t}}({\bm{x}}^{t})$.

  2. 2.

${\bm{z}}_{i}^{t}={\bm{z}}_{i}^{t-1},{\bm{D}}_{i}^{t}=\omega_{t}{\bm{D}}_{i}^{t-1}$,
$\forall i\in[n]\backslash\\{i_{t}\\}$, where $\omega_{t}:=(1+\alpha_{\lceil
t/n\rceil})^{2}$ if $t\mod n=0$ and $1$ otherwise.

  3. 3.

${\bm{Q}}^{t}=\mathrm{BFGS}({\bm{D}}_{i_{t}}^{t-1},(1+\alpha_{\lceil
t/n-1\rceil}){\bm{K}}^{t},{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1})$,

${\bm{D}}_{i_{t}}^{t}=\omega_{t}\mathrm{BFGS}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}),\bar{{\bm{u}}}^{t}({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})))$.

Algorithm 2 provides the pseudo-code for SLIQN. We now consider its per-
iteration complexity. Observe that in the $j+1^{\text{th}}$ epoch, the updates
(2) and (3) are carried out differently for $t\mod n=0$ and $t\mod n\neq 0$.
Specifically, for $t\in\\{nj+1,\dots,nj+n-1\\}$, we can carry out the update
(10) in $\mathcal{O}(d^{2})$ cost using Sherman-Morrison formula (ref.
Appendix F.1) and can compute other iterative updates in $\mathcal{O}(d^{2})$
cost as they consist of a constant number of matrix-vector multiplications.
However, for $t=nj+n$, each ${\bm{D}}_{i}^{t}$ is multiplied with the scaling
factor of $(1+\alpha_{j+1})^{2}$, which incurs a large $\mathcal{O}(nd^{2})$
overhead. Instead, we implement this step by lazily scaling the Hessian
approximations at the iteration in which they are updated while treating all
memoized quantities as if the approximations are already scaled. The details
of the lazy strategy are provided in Appendix F.2.

Algorithm 2 Sharpened Lazy IQN

1: Initialization: Initialize the iterates $\\{{\bm{z}}_{i}^{0}\\}_{i=1}^{n}$
as ${\bm{z}}_{i}^{0}={\bm{x}}^{0}$ and their corresponding hessian
approximations $\\{{\bm{D}}_{i}^{0}\\}_{i=1}^{n}$ as
${\bm{D}}_{i}^{0}=(1+\alpha_{0})^{2}{\bm{I}}_{i}^{0}$

2: while not converged:

3: Update ${\bm{x}}^{t}$ as per (10);

4: Update ${\bm{z}}_{i_{t}}^{t}$ as ${\bm{z}}_{i_{t}}^{t}={\bm{x}}^{t}$;

5: Update ${\bm{Q}}^{t}$ and ${\bm{D}}_{i_{t}}^{t}$ as per (3);

6: Update the tuples with index $i\neq i_{t}$ as per (2);

7: Increment the iteration counter $t$;

8: end while

###### Remark 1

We remark that both IQN and SLIQN exhibit a per-iteration cost of
$\mathcal{O}(d^{2})$ which is in contrast to the $\mathcal{O}(nd^{2})$ cost
for QN methods. However, this efficiency comes with an increased memory cost
of $\mathcal{O}(nd^{2})$. To address this issue we propose a pipelining scheme
in Appendix D, in which we leverage a much larger disk to augment the main
memory by prefetching the data and processing the updates in parallel.

###  5 THEORETICAL ANALYSIS OF ALGORITHM 2

####  5.1 Assumptions

We analyze SLIQN under the assumptions of smoothness, strong convexity, and
Lipschitz continuity of the Hessian. These assumptions are commonly used in
the analysis of Quasi-Newton methods.

###### A1

(Strong convexity and smoothness) The functions $\\{f_{i}\\}_{i=1}^{n}$ are
$\mu$-strongly convex and $L$-smooth, that is
$\frac{\mu}{2}\left\|{\bm{y}}-{\bm{x}}\right\|^{2}\leq
f_{i}({\bm{y}})-f_{i}({\bm{x}})-\langle\nabla
f_{i}({\bm{x}}),{\bm{y}}-{\bm{x}}\rangle$, and
$f_{i}({\bm{y}})-f_{i}({\bm{x}})-\langle\nabla
f_{i}({\bm{x}}),{\bm{y}}-{\bm{x}}\rangle\leq\frac{L}{2}\left\|{\bm{y}}-{\bm{x}}\right\|^{2}$
hold for all ${\bm{x}},{\bm{y}}\in{\mathbb{R}}^{d}$, and for all $i\in[n]$.

###### A2

(Lipschitz continuous Hessian) The Hessians $\\{\nabla^{2}f_{i}\\}_{i=1}^{n}$
are $\tilde{L}$-Lipschitz continuous, that is,
$\left\|\nabla^{2}f_{i}({\bm{y}})-\nabla^{2}f_{i}({\bm{x}})\right\|\leq\tilde{L}\left\|{\bm{y}}-{\bm{x}}\right\|$,
for all ${\bm{x}},{\bm{y}}\in{\mathbb{R}}^{d}$, and for all $i\in[n]$.

The above assumptions also imply that the functions $\\{f_{i}\\}_{i=1}^{n}$
are $M$-strongly self-concordant, which is defined as,
$\nabla^{2}f_{i}({\bm{y}})-\nabla^{2}f_{i}({\bm{x}})\preceq
M\left\|{\bm{y}}-{\bm{x}}\right\|_{\bm{z}}\nabla^{2}f_{i}({\bm{w}})$,
$\forall{\bm{x}},{\bm{y}},{\bm{z}},{\bm{w}}\in{\mathbb{R}}^{d}$ with
$M\coloneqq\tilde{L}\mu^{-\frac{3}{2}}$ (Rodomanov and Nesterov, 2021a, , Ex
4.1).

####  5.2 Convergence Lemmas

We establish the convergence guarantees in three steps: Lemma 1 establishes a
one-step inequality that bounds the residual
$\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|$ in terms of the previous
residuals $\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|$ and the norm
error in the Hessian approximation
$\left\|{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|$,
$\forall i\in[n]$. Lemma 2 uses the result of Lemma 1 to inductively show that
both the residual $\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|$ and the
Hessian approximation error in the $\sigma$ sense, i.e.,
$\sigma({\bm{D}}_{i_{t}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$,
decrease linearly with $\lceil t/n\rceil$. Using the result of Lemma 1 and
Lemma 2, Lemma 3 establishes a mean-superlinear convergence result. We finally
show in Theorem 1 that the residuals can be upper bounded by a superlinearly
convergent sequence.

###### Lemma 1

If Assumptions A1 and A2 hold, the sequence of iterates generated by Algorithm
2 satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|$
$\displaystyle\leq{\frac{\tilde{L}\Gamma^{t-1}}{2}}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|^{2}+{\Gamma^{t-1}}\sum_{i=1}^{n}$
$\displaystyle\left\|{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|,$
(11)

for all $t\geq 1$, where
$\Gamma^{t}:=\left\|\big{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t}\big{)}^{-1}\right\|$.

The proof of this result can be found in Appendix G.1. It is important to note
that our bound in (11) differs from a similar result presented in (Mokhtari et
al.,, 2018, Lemma 2), where they utilize
$\left\|{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{x}}^{\star})\right\|$ instead
of $\left\|{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|$.
This modification helps connect the approximation error in the norm sense,
i.e., $\left\|{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|$
with the error in the $\sigma$ sense, i.e.,
$\sigma({\bm{D}}_{i}^{t-1},\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1}))$. This
connection is crucial for quantifying the improvements achieved by the greedy
updates.

###### Lemma 2

If Assumptions A1 and A2 hold, for any $\rho$ such that
$0<\rho<1-\frac{\mu}{dL}$, there exist positive constants $\epsilon$ and
$\sigma_{0}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon$ and
$\sigma({\bm{I}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}$ for all
$i\in[n]$, the sequence of iterates generated by Algorithm 2 satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\rho^{\lceil\frac{t}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$
(12)

Further, it holds that

$\displaystyle\sigma(\omega_{t}^{-1}{\bm{D}}_{i_{t}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))\leq\big{(}1-\frac{\mu}{dL}\big{)}^{\lceil\frac{t}{n}\rceil}\delta,$
(13)

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

The proof of Lemma 2 can be found in Appendix G.2. Under the hood, the proof
uses induction to show that if the initialized ${\bm{x}}^{0}$ is close to
${\bm{x}}^{\star}$ and ${\bm{D}}_{i}^{0}$ is close to
$\nabla^{2}f_{i}({\bm{x}}^{0})$, then the iterate ${\bm{x}}^{t}$ converges
linearly to ${\bm{x}}^{\star}$ and ${\bm{D}}_{i_{t}}^{t}$ converges linearly
to $\nabla^{2}f_{i}({\bm{z}}_{i_{t}}^{t})$. Our result is stronger than
(Mokhtari et al.,, 2018, Lemma 3) as we establish the linear convergence of
${\bm{D}}_{i}^{t}$, whereas Mokhtari et al., (2018) were only able to
establish that
$\left\|{\bm{D}}_{i}^{t}-\nabla^{2}f_{i}({\bm{x}}^{\star})\right\|$ does not
grow with $t$. Moreover, Mokhtari et al., (2018) did not guarantee convergence
of ${\bm{D}}_{i}^{t}$. Also, there is no equivalent convergence result
presented in the analysis of IGS.

###### Lemma 3

If Assumptions A1, A2 hold, then there exist positive constants $\epsilon$ and
$\sigma_{0}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon$ and
$\sigma({\bm{I}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}$ for all
$i\in[n]$, the sequence of iterates generated by Algorithm 2 satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\big{(}1-\frac{\mu}{dL}\big{)}^{\lceil\frac{t}{n}\rceil}\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{x}}^{t-i}-{\bm{x}}^{\star}\right\|.$

The proof of Lemma 3 can be found in Appendix G.3. The main idea behind the
proof is to substitute the linear convergence results, specifically (12) and
(13) from Lemma 2, back into the result from Lemma 1. By doing so, the first
term on the right-hand side of (11) converges quadratically, while the second
term converges superlinearly, which proves the result.

Our result is markedly different from (Mokhtari et al.,, 2018, Theorem 6)
which proves asymptotic superlinear convergence of IQN. Their proof is based
on a variant of the Dennis-Moré theorem to show that superlinear convergence
kicks asymptotically. Since their proof is existential, unlike Lemma 3, it
cannot be used to derive an explicit rate of superlinear convergence.

###### Theorem 1

If Assumptions A1, A2 hold, then there exist positive constants $\epsilon$ and
$\sigma_{0}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon$ and
$\sigma({\bm{I}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}$ for all
$i\in[n]$, and for the sequence of iterates generated by Algorithm 2, there
exists a sequence $\\{\zeta^{k}\\}$, $k\geq 1$ such that
$\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\zeta^{\lfloor\frac{t-1}{n}\rfloor}$
for all $t\geq 1$ and $\\{\zeta^{k}\\}$ satisfies,

$\displaystyle\zeta^{k}\leq\epsilon\big{(}1-\frac{\mu}{dL}\big{)}^{\frac{(k+2)(k+1)}{2}}.$
(14)

The proof of Theorem 1 involves the construction of a sequence that provides
an upper bound on the residual $\|{{\bm{x}}}^{t}-{{\bm{x}}}^{\star}\|$, and
can be found in Appendix G.4.

###### Remark 2

It is instructive to compare our convergence rate with that of IGS (Gao et
al.,, 2020, Theorem 3). According to their result, the rate is given as
$\|{{\bm{x}}}^{t}-{{\bm{x}}}^{\star}\|\leq\left(1-\frac{\mu}{dL}\right)^{\frac{k(k+1)}{2}}r^{k_{0}}||{{\bm{x}}}^{0}-{{\bm{x}}}^{\star}||$,
$\forall t\geq 1$,
$r\in(0,1),k=\left\lfloor\frac{t-1}{n}\right\rfloor+1-k_{0}$ and $k_{0}$ is a
constant such that $\left(1-\frac{\mu}{dL}\right)^{k_{0}}D\leq 1$. The
parameter $D$ depends on the underlying objective function. Observe that their
superlinear rate only takes effect after
$\left\lfloor\frac{t-1}{n}\right\rfloor\geq k_{0}-1$, and $k_{0}$ could
potentially be large, though it is not possible to infer the bounds on $k_{0}$
from Gao et al., (2020). In contrast, our convergence rate guarantees
superlinear convergence right from the first iteration.

###  6 NUMERICAL EXPERIMENTS

####  6.1 Quadratic Function Minimization

We begin with a comparative analysis of the empirical performance of SLIQN,
IQN, Sharpened BFGS (SBFGS) Jin et al., (2022), and IGS on a synthetic
quadratic minimization problem.

Problem Definition: We consider the function
$f({\bm{x}})=\frac{1}{n}\sum_{i=1}^{n}\big{(}\frac{1}{2}\langle{\bm{x}},{\bm{A}}_{i}{\bm{x}}\rangle+\langle\bm{b}_{i},{\bm{x}}\rangle\big{)},$
where ${\bm{A}}_{i}\succ\mathbf{0}$, and $\bm{b}_{i}\in{\mathbb{R}}^{d}$,
$\forall i\in[n]$. The detailed generation scheme for
${\bm{A}}_{i},{\bm{b}}_{i}$ can be found in Appendix I.1.

Experiments and Inference: We study the performance of the algorithms on two
extreme cases: $d\gg n$ (Fig. 1(a)), and $n\gg d$ (Fig. 1(b)). In each case,
we plot the normalized error
$\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|/\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$
against the number of effective passes or epochs. We see that in both these
cases, SLIQN outperforms IGS, IQN, and SBFGS. We also observe that in the case
where $n\gg d$, IGS outperforms IQN, whereas in the case where $d\gg n$, IQN
surpasses IGS. This is because IGS devotes the initial $\mathcal{O}(d)$
iterations to constructing a precise Hessian approximation, after which its
fast convergence phase kicks in. On the other hand, since IQN takes the
descent step in the Newton direction, its Hessian approximation is never
precise and therefore its normalized error decreases at more or less a
“consistent” rate. SLIQN combines the strengths of both IQN and IGS: during
the initial iterations when its Hessian approximation is not accurate enough,
the classical BFGS updates are responsible for the progress made. In the later
iterations, when the Hessian has been sufficiently well approximated, its fast
convergence phase kicks in.

(a) $n=20,d=500$

(b) $n=50000,d=10$

Figure 1: Normalized error vs. number of effective passes for the quadratic
minimization problem Figure 2: Normalized error vs. number of effective passes
for regularized logistic loss minimization

####  6.2 Regularized Logistic Regression

We now compare the performance of SLIQN against IQN, IGS, NIM, and SN on the
regularized logistic regression task given by

$\displaystyle\min_{{\bm{x}}}f({\bm{x}})$
$\displaystyle\coloneqq\frac{1}{N}\sum_{i=1}^{N}\big{(}y_{i}\log(1+e^{-\langle{\bm{x}},{\bm{z}}_{i}\rangle})+(1-y_{i})\log(1+e^{\langle{\bm{x}},{\bm{z}}_{i}\rangle})\big{)}+\frac{\lambda}{2}\left\|{\bm{x}}\right\|^{p},$
(15)

where $\\{{\bm{z}}_{i}\\}_{i=1}^{N}$ are the training samples,
$\\{y_{i}\\}_{i=1}^{N}$ are their corresponding binary labels and $p$ is set
at $2.1$. It is easy to observe that $f({\bm{x}})$ is smooth, strongly convex
and has a Lipschitz continuous Hessian, thereby satisfying Assumptions A1, A2.
We compare the algorithms across $9$ datasets with a large variation in the
values of $n$ and $d$. Each algorithm is initialization with the same iterate
and the regularization parameter is set as $\lambda={1}/{N}$. For SN, we set
the mini-batch size $\tau=N$ since that is the regime it works best in. Please
refer to Appendix I for a complete experimental setup. We observe in Figure 2
that SLIQN outperforms IGS, IQN, and SN on each of the $9$ datasets from
LIBSVM by Chang and Lin, (2011). This supports our claim that SLIQN offers the
best of both, IQN and IGS. Furthermore, we observe that while NIM outperforms
SLIQN, their performance remains comparable. It is important to note that NIM
utilizes the full Hessian information for the descent step and is an
$\mathcal{O}(d^{3})$ algorithm, while SLIQN has a per-iteration complexity of
$\mathcal{O}(d^{2})$. Thus, these results underscore the superiority of SLIQN
over other incremental QN style methods.

###  7 CONCLUSION AND FUTURE WORK

We introduced the SLIQN method for minimizing finite-sum problems. SLIQN
enjoys the best known incremental rate of
$\mathcal{O}((1-\tfrac{\mu}{dL})^{t^{2}/n^{2}})$, has a $\mathcal{O}(d^{2})$
per-iteration cost, an explicit superlinear convergence rate, and exhibits a
superior empirical performance compared to several other incremental and
stochastic QN methods. The key novelty in the proposed algorithm is the
construction of modified update rules using a clever multiplicative factor and
a lazy propagation strategy. We back up our empirical results with a
comprehensive theory that explains the superior performance of SLIQN. The
convergence rate of SLIQN is locally superlinear; analyzing the global
convergence of the proposed algorithm remains as a future work.

### References

  * 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 and Lin, (2011) Chang, C.-C. and Lin, C.-J. (2011).  Libsvm: A library for support vector machines.  ACM Trans. Intell. Syst. Technol., 2(3).
  * Chang et al., (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.
  * Chen et al., (2022) Chen, J., Yuan, R., Garrigos, G., and Gower, R. M. (2022).  San: Stochastic average newton algorithm for minimizing finite sums.  In Camps-Valls, G., Ruiz, F. J. R., and Valera, I., editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 279–318. PMLR.
  * Defazio et al., (2014) Defazio, A., Bach, F. R., and Lacoste-Julien, S. (2014).  Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.  In Neural Information Processing Systems.
  * Derezinski, (2023) Derezinski, M. (2023).  Stochastic variance-reduced newton: Accelerating finite-sum minimization with large batches.  In OPT 2023: Optimization for Machine Learning.
  * Derezinski et al., (2021) Derezinski, M., Lacotte, J., Pilanci, M., and Mahoney, M. W. (2021).  Newton-less: Sparsification without trade-offs for the sketched newton update.  Advances in Neural Information Processing Systems, 34:2835–2847.
  * Gao et al., (2020) Gao, Z., Koppel, A., and Ribeiro, A. (2020).  Incremental greedy bfgs: An incremental quasi-newton method with explicit superlinear rate.  In Adv. Neural Inf. Process. Syst. 12th OPT Workshop Optim. Mach. Learn.
  * Gonen et al., (2016) Gonen, A., Orabona, F., and Shalev-Shwartz, S. (2016).  Solving ridge regression using sketched preconditioned svrg.  In International conference on machine learning, pages 1397–1405. PMLR.
  * Jin et al., (2022) Jin, Q., Koppel, A., Rajawat, K., and Mokhtari, A. (2022).  Sharpened quasi-newton methods: Faster superlinear rate and larger local convergence neighborhood.  In International Conference on Machine Learning, pages 10228–10250. PMLR.
  * Jin and Mokhtari, (2022) Jin, Q. and Mokhtari, A. (2022).  Non-asymptotic superlinear convergence of standard quasi-newton methods.  Mathematical Programming, pages 1–49.
  * Johnson and Zhang, (2013) Johnson, R. and Zhang, T. (2013).  Accelerating stochastic gradient descent using predictive variance reduction.  Advances in neural information processing systems, 26.
  * Kovalev et al., (2019) Kovalev, D., Mishchenko, K., and Richtárik, P. (2019).  Stochastic newton and cubic newton methods with simple local linear-quadratic rates.  arXiv preprint arXiv:1912.01597.
  * Li et al., (2021) Li, J., Bian, S., Zeng, A., Wang, C., Pang, B., Liu, W., and Lu, C. (2021).  Human pose regression with residual log-likelihood estimation.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 11025–11034.
  * Li et al., (2020) Li, J., Zhou, P., Xiong, C., Socher, R., and Hoi, S. C. H. (2020).  Prototypical contrastive learning of unsupervised representations.  CoRR, abs/2005.04966.
  * Liu et al., (2019) Liu, Y., Feng, F., and Yin, W. (2019).  Acceleration of svrg and katyusha x by inexact preconditioning.  In International Conference on Machine Learning, pages 4003–4012. PMLR.
  * Mokhtari et al., (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 et al., (2016) Moritz, P., Nishihara, R., and Jordan, M. (2016).  A linearly-convergent stochastic l-bfgs algorithm.  In Artificial Intelligence and Statistics, pages 249–258. PMLR.
  * Nocedal and Wright, (1999) Nocedal, J. and Wright, S. J. (1999).  Numerical optimization.  Springer-Verlag, New York, NY, 2 edition.
  * Pilanci and Wainwright, (2016) Pilanci, M. and Wainwright, M. J. (2016).  Iterative hessian sketch: Fast and accurate solution approximation for constrained least-squares.  The Journal of Machine Learning Research, 17(1):1842–1879.
  * Pilanci and Wainwright, (2017) Pilanci, M. and Wainwright, M. J. (2017).  Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence.  SIAM Journal on Optimization, 27(1):205–245.
  * 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, pages 2597–2605. PMLR.
  * (25) Rodomanov, A. and Nesterov, Y. (2021a).  Greedy quasi-newton methods with explicit superlinear convergence.  SIAM Journal on Optimization, 31(1):785–811.
  * (26) 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.
  * (27) Rodomanov, A. and Nesterov, Y. (2021c).  Rates of superlinear convergence for classical quasi-newton methods.  Mathematical Programming, pages 1–32.
  * Song and Ermon, (2020) Song, Y. and Ermon, S. (2020).  Improved techniques for training score-based generative models.  ArXiv, abs/2006.09011.
  * Wu et al., (2018) Wu, J., Huang, W., Huang, J., and Zhang, T. (2018).  Error compensated quantized sgd and its applications to large-scale distributed optimization.  In International Conference on Machine Learning, pages 5325–5333. PMLR.

## Appendix

###  Appendix A ESTABLISHED RESULTS

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

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

$\displaystyle\frac{1}{1+\left\|{\bm{A}}\right\|}<\left\|({\bm{I}}+{\bm{A}})^{-1}\right\|<\frac{1}{1-\left\|{\bm{A}}\right\|}.$

######  Proposition A.1 (Sherman-Morrison Formula)

Let ${\bm{A}}\in{\mathbb{R}}^{d\times d}$ be an invertible matrix. Then, for
all vectors ${\bm{u}},{\bm{v}}\in{\mathbb{R}}^{d}$, we have

$\displaystyle({\bm{A}}+{\bm{u}}{\bm{v}}^{\mathsf{T}})^{-1}={\bm{A}}^{-1}-\frac{{\bm{A}}^{-1}{\bm{u}}{\bm{v}}^{\mathsf{T}}{\bm{A}}^{-1}}{1+\langle{\bm{v}},{\bm{A}}^{-1}{\bm{u}}\rangle}.$
(16)

######  Lemma A.2 (Lemma 2.1, Lemma 2.2 Rodomanov and Nesterov, 2021c )

Consider positive definite matrices ${\bm{A}},{\bm{G}}\in\mathbb{R}^{d\times
d}$ and suppose ${\bm{G}}_{+}:=\mathrm{BFGS}({\bm{G}},{\bm{A}},{\bm{u}})$,
where ${\bm{u}}\neq{\mathbf{0}}$. Then, the following results hold:

  1. 1.

For any constants $\xi,\eta\geq 1$, we have

$\displaystyle\frac{1}{\xi}{\bm{A}}\preceq{\bm{G}}\preceq\eta{\bm{A}}\implies\frac{1}{\xi}\preceq{\bm{G}}_{+}\preceq\eta{\bm{A}}.$

  2. 2.

If ${\bm{A}}\preceq{\bm{G}}$, then we have

$\displaystyle\sigma({\bm{A}},{\bm{G}})\geq\sigma({\bm{A}},{\bm{G}}_{+}).$

######  Lemma A.3 (Lemma 4.2 Rodomanov and Nesterov, 2021a )

Suppose an objective function $f({\bm{x}})$ is strongly self-concordant with
constant $M>0$. Consider ${\bm{x}},{\bm{y}}\in\mathbb{R}^{d}$,
$r:=\left\|{\bm{y}}-{\bm{x}}\right\|_{{\bm{x}}}$, and
${\bm{K}}:=\int_{0}^{1}\nabla^{2}f({\bm{x}}+\tau({\bm{y}}-{\bm{x}}))d\tau$.
Then, we have that

$\displaystyle\frac{\nabla^{2}f({\bm{x}})}{1+Mr}\preceq\nabla^{2}f({\bm{y}})\preceq(1+Mr)\nabla^{2}f({\bm{x}}),$
$\displaystyle\frac{\nabla^{2}f({\bm{x}})}{1+\frac{Mr}{2}}\preceq{\bm{K}}\preceq(1+\frac{Mr}{2})\nabla^{2}f({\bm{x}}),$
$\displaystyle\frac{\nabla^{2}f({\bm{y}})}{1+\frac{Mr}{2}}\preceq{\bm{K}}\preceq(1+\frac{Mr}{2})\nabla^{2}f({\bm{y}}).$

######  Lemma A.4 (Theorem 2.5 Rodomanov and Nesterov, 2021a )

Consider positive definite matrices ${\bm{A}},{\bm{G}}\in\mathbb{R}^{d\times
d}$ such that ${\bm{A}}\preceq{\bm{G}}$ and $\mu{\bm{I}}\preceq{\bm{A}}\preceq
L{\bm{I}}$ for constants $\mu,L>0$. Suppose
${\bm{G}}_{+}:=\mathrm{BFGS}({\bm{G}},{\bm{A}},\bar{{\bm{u}}}({\bm{G}},{\bm{A}}))$,
where $\bar{{\bm{u}}}({\bm{G}},{\bm{A}})$ (4) is the greedy vector of
${\bm{G}}$ with respect to ${\bm{A}}$. Then, the following holds:

$\displaystyle\sigma({\bm{G}}_{+},{\bm{A}})\leq\big{(}1-\frac{\mu}{dL}\big{)}\sigma({\bm{G}},{\bm{A}}).$

###  Appendix B SUPPORTING LEMMAS

###### Lemma B.1

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

$\displaystyle\left\|{\bm{G}}-{\bm{A}}\right\|\leq
L\sigma({\bm{G}},{\bm{A}}).$

##### Proof:

For any positive definite matrix ${\bm{X}}\in\mathbb{R}^{d\times d}$, let
$\lambda_{\max}({\bm{X}})$ denote its maximum eigenvalue and let
$\sum_{i=1}^{d}\lambda_{i}({\bm{X}})$ denote the sum of all of its
eigenvalues. We can bound $\tfrac{1}{L}\left\|{\bm{G}}-{\bm{A}}\right\|$ as

$\displaystyle\tfrac{1}{L}\left\|{\bm{G}}-{\bm{A}}\right\|$
$\displaystyle\overset{\text{def}}{=}\tfrac{1}{L}\lambda_{\max}({\bm{G}}-{\bm{A}})\leq\tfrac{1}{L}\sum_{i=1}^{d}\lambda_{i}({\bm{G}}-{\bm{A}})\overset{\text{def}}{=}\tfrac{1}{L}\mathrm{Tr}({\bm{G}}-{\bm{A}}).$

Recall from the premise, we have ${\bm{A}}\preceq L{\bm{I}}$, which implies
that ${\bm{I}}\preceq L{\bm{A}}^{-1}$. Therefore,

$\displaystyle\sigma({\bm{G}},{\bm{A}})=\langle{\bm{A}}^{-1},{\bm{G}}-{\bm{A}}\rangle\geq\tfrac{1}{L}\langle{\bm{I}},{\bm{G}}-{\bm{A}}\rangle=\tfrac{1}{L}\mathrm{Tr}({\bm{G}}-{\bm{A}})\geq\tfrac{1}{L}\left\|{\bm{G}}-{\bm{A}}\right\|.$

This completes the proof. $\square$

###### Lemma B.2

Let $f\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a real valued function
that is $\mu$-strongly convex, $L$-smooth, and $M$-strongly self-concordant.
Let ${\bm{x}},{\bm{x}}_{+}\in\mathbb{R}^{d}\backslash\\{\mathbf{0}\\}$ and
${\bm{B}}$ be a matrix such that ${\bm{B}}\succeq\nabla^{2}f({\bm{x}})$.
Define the constant $r:=\left\|{\bm{x}}_{+}-{\bm{x}}\right\|_{\bm{x}}$ and the
matrix ${\bm{P}}:=(1+\tfrac{Mr}{2})^{2}{\bm{B}}$. Consider the following BFGS
updates:

$\displaystyle{\bm{Q}}$
$\displaystyle:=\mathrm{BFGS}({\bm{P}},(1+\tfrac{Mr}{2}){\bm{K}},{\bm{x}}_{+}-{\bm{x}}),$
$\displaystyle{\bm{B}}_{+}$
$\displaystyle:=\mathrm{BFGS}({\bm{Q}},\nabla^{2}f({\bm{x}}_{+}),\bar{{\bm{u}}}({\bm{Q}},\nabla^{2}f({\bm{x}}_{+}))).$

Here, the matrix
${\bm{K}}:=\int_{0}^{1}\nabla^{2}f({\bm{x}}+\tau({\bm{x}}_{+}-{\bm{x}}))d\tau$
and the vector $\bar{{\bm{u}}}({\bm{Q}},\nabla^{2}f({\bm{x}}_{+}))$ is the
greedy vector 4. Then, ${\bm{B}}_{+}\succeq\nabla^{2}f({\bm{x}}_{+})$ and

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))\leq\big{(}1-\tfrac{\mu}{dL}\big{)}\bigg{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{4}\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+d\big{(}1+\tfrac{Mr}{2}\big{)}^{4}-d\bigg{)}.$

##### Proof:

We begin by analyzing the first BFGS update. Since
${\bm{B}}\succeq\nabla^{2}f({\bm{x}})$, we have the following:

$\displaystyle{\bm{B}}\succeq\nabla^{2}f({\bm{x}})\overset{\text{def}}{\implies}{\bm{P}}\succeq\big{(}1+\frac{Mr}{2}\big{)}^{2}f({\bm{x}})\stackrel{{\scriptstyle{\textit{Lem}}.\ref{integral_lemma}}}{{\succeq}}\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}.$

Since ${\bm{P}}\succeq\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}$, the metric
$\sigma\big{(}{\bm{P}},(1+\frac{Mr}{2}){\bm{K}}\big{)}$ is well defined.
Applying Lemma A.2, we obtain

$\displaystyle{\bm{Q}}=\mathrm{BFGS}({\bm{P}},\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}},{\bm{x}}_{+}-{\bm{x}})\succeq\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}.$

Applying Lemma A.3 to relate ${\bm{K}}$ and $\nabla^{2}f({\bm{x}}_{+})$, we
obtain

$\displaystyle{\bm{Q}}\succeq\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}\stackrel{{\scriptstyle{\textit{Lem}}.\ref{integral_lemma}}}{{\succeq}}\nabla^{2}f({\bm{x}}_{+}).$

We now begin analyzing the second BFGS update. Since
${\bm{Q}}\succeq\nabla^{2}f({\bm{x}}_{+})$, applying Lemma A.2, we obtain

$\displaystyle{\bm{B}}_{+}=\mathrm{BFGS}({\bm{Q}},\nabla^{2}f({\bm{x}}_{+}),\bar{{\bm{u}}}({\bm{Q}},\nabla^{2}f({\bm{x}}_{+}))\succeq\nabla^{2}f({\bm{x}}_{+}),$

which completes the proof of the first part. Applying Lemma A.4, we obtain

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))\leq\left(1-\tfrac{\mu}{dL}\right)\sigma({\bm{Q}},\nabla^{2}f({\bm{x}}_{+})).$

Define $c\coloneqq\left(1-\tfrac{\mu}{dL}\right)$ for brevity. We are now
ready to show the second result. Observe that

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))$
$\displaystyle\leq(1-c)\sigma({\bm{Q}},\nabla^{2}f({\bm{x}}_{+})),$
$\displaystyle=(1-c)(\langle\nabla^{2}f({\bm{x}}_{+})^{-1},{\bm{Q}}\rangle-d),$
$\displaystyle\overset{(a)}{\leq}(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}\langle{\bm{K}}^{-1},{\bm{Q}}\rangle-d\big{)},$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{2}\langle\tilde{{\bm{K}}}^{-1},{\bm{Q}}\rangle-d\big{)},$
(17)

where $(a)$ follows since

$\big{(}1+\tfrac{Mr}{2}\big{)}\nabla^{2}f({\bm{x}}_{+})\stackrel{{\scriptstyle{\textit{Lem}}.\ref{integral_lemma}}}{{\succeq}}{\bm{K}}\iff\nabla^{2}f({\bm{x}}_{+})^{-1}\preceq\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}^{-1}.$

In $(b)$, we have defined
$\tilde{{\bm{K}}}\coloneqq\big{(}1+\frac{Mr}{2}\big{)}{\bm{K}}$. Recall that
we have already established ${\bm{P}}\succeq\tilde{{\bm{K}}}$. Applying Lemma
A.2, we obtain

$\displaystyle\sigma({\bm{Q}},\tilde{{\bm{K}}})\leq\sigma({\bm{P}},\tilde{{\bm{K}}})\iff\sigma({\bm{Q}},{\bm{K}})\leq\sigma({\bm{P}},{\bm{K}}).$

Continuing from (17), we obtain

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))$
$\displaystyle\leq(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{2}\langle\tilde{{\bm{K}}}^{-1},{\bm{P}}\rangle-d\big{)},$
$\displaystyle=(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{4}\langle\tilde{{\bm{K}}}^{-1},{\bm{B}}\rangle-d\big{)},$
$\displaystyle=(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{3}\langle{{\bm{K}}}^{-1},{\bm{B}}\rangle-d\big{)},$
$\displaystyle\overset{(a)}{\leq}(1-c)\big{(}\big{(}1+\tfrac{Mr}{2}\big{)}^{4}\langle{\nabla^{2}f({\bm{x}})}^{-1},{\bm{B}}\rangle-d\big{)},$
(18)

where $(a)$ follows from
${\bm{K}}^{-1}\preceq\big{(}1+\tfrac{Mr}{2}\big{)}\nabla^{2}f({\bm{x}})^{-1}$
($\because$ Lemma A.3). Finally, setting
$\langle{\nabla^{2}f({\bm{x}})}^{-1},{\bm{B}}\rangle=\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+d$
completes the proof. $\square$

###### Corollary B.1

Under the notation established in Lemma B.2, let $\alpha\in\mathbb{R}_{+}$ be
an upper bound on $r$. Then, we have the following:

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))\leq(1-\tfrac{\mu}{dL})e^{2M\alpha}\big{(}\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+2Md\alpha\big{)}.$

##### Proof:

From Lemma B.2, we have

$\displaystyle\sigma({\bm{B}}_{+},\nabla^{2}f({\bm{x}}_{+}))$
$\displaystyle\leq\big{(}1-\tfrac{\mu}{dL}\big{)}\big{(}1+\tfrac{M\alpha}{2}\big{)}^{4}\big{(}\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+d\big{(}1-\tfrac{1}{\big{(}1+\tfrac{M\alpha}{2}\big{)}^{4}})),$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\big{(}1-\tfrac{\mu}{dL}\big{)}e^{2M\alpha}\big{(}\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+d(1-e^{-2M\alpha})\big{)},$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\big{(}1-\tfrac{\mu}{dL}\big{)}e^{2M\alpha}\big{(}\sigma({\bm{B}},\nabla^{2}f({\bm{x}}))+2Md\alpha\big{)},$

where $(a)$ follows by the inequality $1+x\leq e^{x},\forall x\in\mathbb{R}$
and $(b)$ follows from the inequality $1-e^{-x}\leq x,\forall x>0$. This
completes the proof. $\square$

###### Lemma B.3

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

$\displaystyle\left\|{\bm{x}}^{k}-\tilde{{\bm{x}}}\right\|\leq{\gamma}^{k}\left\|{\bm{x}}^{0}-\tilde{{\bm{x}}}\right\|.$

Define the constant
$r_{k}:=\left\|{\bm{x}}^{k}-{\bm{x}}^{k-1}\right\|_{{\bm{x}}^{k-1}}$ for every
$k$. Let ${\bm{B}}^{0}$ be a matrix such that it satisfies
${\bm{B}}^{0}\succeq\nabla^{2}f({\bm{x}}^{0})$. Consider the following BFGS
updates:

$\displaystyle{\bm{Q}}^{k}$
$\displaystyle:=\mathrm{BFGS}({\bm{P}}^{k-1},(1+\tfrac{Mr_{k}}{2}){\bm{K}}^{k},{\bm{x}}^{k}-{\bm{x}}^{k-1}),$
$\displaystyle{\bm{B}}^{k}$
$\displaystyle:=\mathrm{BFGS}({\bm{Q}}^{k},\nabla^{2}f({\bm{x}}^{k}),\bar{{\bm{u}}}({\bm{Q}}^{k},\nabla^{2}f({\bm{x}}^{k}))),$

where ${\bm{P}}^{k-1}:=\big{(}1+\tfrac{Mr_{k}}{2}\big{)}^{2}{\bm{B}}^{k-1}$,
${\bm{K}}^{k}:=\int_{0}^{1}\nabla^{2}f({\bm{x}}^{k-1}+\tau({\bm{x}}^{k}-{\bm{x}}^{k-1}))d\tau$,
and $\bar{{\bm{u}}}({\bm{Q}}^{k},\nabla^{2}f({\bm{x}}^{k}))$ is the greedy
vector 4. Then, the following holds for all $k\in[T]$:

$\displaystyle\sigma({\bm{B}}^{k},\nabla^{2}f({\bm{x}}^{k}))\leq(1-\tfrac{\mu}{dL})^{k}e^{\tfrac{4M\sqrt{L}||{{\bm{x}}^{0}-\tilde{{\bm{x}}}}||}{1-\gamma}}\bigg{(}\sigma({\bm{B}}^{0},\nabla^{2}f({\bm{x}}^{0}))+\left\|{\bm{x}}^{0}-\tilde{{\bm{x}}}\right\|\tfrac{4Md\sqrt{L}}{1-(1-\tfrac{\mu}{dL})^{-1}\gamma}\bigg{)}.$
(19)

##### Proof:

From Lemma B.2, it can be shown that
${\bm{B}}^{k}\succeq\nabla^{2}f({\bm{x}}^{k})$ for $k=1,\dots,T$. Therefore,
$\sigma({\bm{B}}^{k},\nabla^{2}f({\bm{x}}^{k}))$ is well defined. We introduce
the notation $\sigma_{k}:=\sigma({\bm{B}}^{k},\nabla^{2}f({\bm{x}}^{k}))$,
$d_{k}:=\left\|{\bm{x}}^{k}-\tilde{{\bm{x}}}\right\|$, and
$c:=\tfrac{\mu}{dL}$ for simplicity.

To apply Corollary B.1, we need an upper bound on $r_{k}$. This is trivial via
the Triangle inequality

$\displaystyle
r_{k}=\left\|{\bm{x}}^{k}-{\bm{x}}^{k-1}\right\|_{{\bm{x}}^{k-1}}$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{\leq}}\left\|{\bm{x}}^{k}-\tilde{{\bm{x}}}\right\|_{{\bm{x}}^{k-1}}+\left\|{\bm{x}}^{k-1}-\tilde{{\bm{x}}}\right\|_{{\bm{x}}^{k-1}}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\sqrt{L}\big{(}\left\|{\bm{x}}^{k}-\tilde{{\bm{x}}}\right\|+\left\|{\bm{x}}^{k-1}-\tilde{{\bm{x}}}\right\|\big{)}\leq
2\sqrt{L}\gamma^{k-1}\left\|{\bm{x}}^{0}-\tilde{{\bm{x}}}\right\|,$

where $(a)$ follows since $\nabla^{2}f({\bm{x}})\preceq LI$. Therefore
$\alpha_{k}=2\sqrt{L}\gamma^{k-1}\left\|{\bm{x}}^{0}-\tilde{{\bm{x}}}\right\|$
is an upper bound on $r_{k}$. Applying Corollary B.1, we obtain

$\displaystyle\sigma_{k}$
$\displaystyle\leq(1-c)e^{2M\alpha_{k}}(\sigma_{k-1}+2Md\alpha_{k})$
$\displaystyle\leq(1-c)e^{4M\sqrt{L}d_{0}\gamma^{k-1}}(\sigma_{k-1}+4Md\sqrt{L}d_{0}\gamma^{k-1}).$
(20)

We solve the recursion (20) as follows:

$\displaystyle\sigma_{k}$
$\displaystyle\leq(1-c)e^{4M\sqrt{L}d_{0}\gamma^{k-1}}(\sigma_{k-1}+4Md\sqrt{L}d_{0}\gamma^{k-1}),$
$\displaystyle\leq(1-c)^{2}e^{4M\sqrt{L}d_{0}(\gamma^{k-1}+\gamma^{k-2})}\sigma_{k-2}+4Md\sqrt{L}d_{0}\big{(}\gamma^{k-2}(1-c)^{2}e^{4M\sqrt{L}d_{0}(\gamma^{k-1}+\gamma^{k-2})}$
$\displaystyle\hskip
227.62204pt+\gamma^{k-1}(1-c)e^{4M\sqrt{L}d_{0}\gamma^{k-1}}\big{)},$
$\displaystyle\leq(1-c)^{k}e^{4M\sqrt{L}d_{0}\sum_{j=0}^{k-1}\gamma^{j}}\sigma_{0}+4Md\sqrt{L}d_{0}\big{(}\sum_{j=0}^{k-1}\gamma^{j}(1-c)^{k-j}e^{4M\sqrt{L}d_{0}\sum_{i=1}^{k-j}\gamma^{k-i}}\big{)},$
$\displaystyle\leq(1-c)^{k}e^{4M\sqrt{L}d_{0}\sum_{j=0}^{\infty}\gamma^{j}}\sigma_{0}+4Md\sqrt{L}d_{0}e^{4M\sqrt{L}d_{0}\sum_{i=0}^{\infty}\gamma^{i}}\sum_{j=0}^{k-1}\gamma^{j}(1-c)^{k-j},$
$\displaystyle\leq(1-c)^{k}e^{\tfrac{4M\sqrt{L}d_{0}}{1-\gamma}}\big{(}\sigma_{0}+\tfrac{4Md\sqrt{L}d_{0}}{1-\tfrac{\gamma}{1-c}}\big{)}.$

This completes the proof. $\square$

###### Remark 3

It can be concluded from the proof of Lemma B.3 that redefining
$r_{k}:=2\sqrt{L}\gamma^{k-1}d_{0}$ (which is an upper bound on
$\left\|{\bm{x}}^{k}-{\bm{x}}^{k-1}\right\|_{{\bm{x}}^{k-1}}$), the results of
Lemma B.3 remain unchanged.

###  Appendix C EFFICIENT IMPLEMENTATION OF IQN

Recall that, that at time $t$, IQN updates

$\displaystyle{\bm{B}}_{i_{t}}^{t}=\mathrm{BFGS}({\bm{B}}_{i_{t}}^{t-1},{\bm{K}}^{t},{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}).$

Define
$\bm{\phi}^{t}:=\sum_{i=1}^{n}{\bm{B}}_{i}^{t}{\bm{z}}_{i}^{t}-\sum_{i=1}^{n}\nabla
f_{i}({\bm{z}}_{i}^{t})$. At time $t$, since IQN only updates the tuple with
index $i_{t}$, we have

$\displaystyle\bm{\phi}^{t}=\bm{\phi}^{t-1}+\big{(}{\bm{B}}_{i_{t}}^{t}{\bm{z}}_{i_{t}}^{t}-{\bm{B}}_{i_{t}}^{t-1}{\bm{z}}_{i_{t}}^{t-1}\big{)}-\big{(}\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})\big{)}.$ (21)

Therefore, given access to $\bm{\phi}^{t-1}$, we can compute $\bm{\phi}^{t}$
in $\mathcal{O}(d^{2})$ time. This updating scheme can be implemented
iteratively, where we only evaluate $\bm{\phi}^{0}$ explicitly and evaluate
$\bm{\phi}^{t}$, for all $t\geq 1$ by (21). It only remains to compute
$\big{(}\bar{{\bm{B}}}^{t}\big{)}^{-1}$, where
$\bar{{\bm{B}}}^{t}:=\big{(}\sum_{i=1}^{n}{\bm{B}}_{i}^{t}\big{)}$. This can
be done by applying the Sherman-Morrison formula (16) twice to the matrix on
the right (22).

$\displaystyle\bar{{\bm{B}}}^{t}=\bar{{\bm{B}}}^{t-1}+{\bm{B}}_{i_{t}}^{t}-{\bm{B}}_{i_{t}}^{t-1}=\bar{{\bm{B}}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}{{\bm{s}}_{i_{t}}^{t}}^{{\mathsf{T}}}{\bm{B}}_{i_{t}}^{t-1}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle},$
(22)

where
${\bm{y}}_{i_{t}}^{t}={\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1},{\bm{s}}_{i_{t}}^{t}=\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})$.
Applying(16) twice, we get (23) and (24)

$\displaystyle(\bar{{\bm{B}}}^{t})^{-1}={{\bm{Z}}}^{t}+\frac{{{\bm{Z}}}^{t}({\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})({\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{{\mathsf{T}}}{{\bm{Z}}}^{t}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle-\langle({\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}),{\bm{Z}}^{t}({\bm{B}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})\rangle},$
(23)

where ${\bm{Z}}^{t}$ is given by

$\displaystyle{\bm{Z}}^{t}=(\bar{{\bm{B}}}^{t-1})^{-1}-\frac{(\bar{{\bm{B}}}^{t-1})^{-1}{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}(\bar{{\bm{B}}}^{t-1})^{-1}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle+\langle{{\bm{y}}_{i_{t}}^{t}},(\bar{{\bm{B}}}^{t-1})^{-1}{{\bm{y}}_{i_{t}}^{t}}\rangle}.$
(24)

By iteratively implementing this scheme, we only need to compute
$(\bar{{\bm{B}}}^{0})^{-1}$ explicitly and $(\bar{{\bm{B}}}^{t})^{-1}$, for
all $t\geq 1$ by the Sherman-Morrison formula.

###  Appendix D LOW MEMORY IMPLEMENTATION OF INCREMENTAL METHODS

In this section, we illustrate how incremental methods can be effectively
implemented for large-scale real-world scenarios which are characterized by
substantial memory requirements of $\mathcal{O}(nd^{2})$. Our solution
leverages the disk, which offers significantly larger storage capacity
compared to main memory but comes with an increased load-store latency. To
mitigate this latency issue, we employ a pipelining scheme. In this scheme, we
partition the data into blocks and simultaneously run compute operations on
one block while performing load-store operations on the blocks adjacent to it.
This parallelization effectively extends the main memory capacity to the
available the disk size, all the while avoiding its larger latency.

Formally, let the available main memory capacity be $g$ GB, the number of data
samples be $n$, the dimensionality of the data be $d$, and the space
requirement for each sample be $s$. We assume that the disk is sufficiently
large to store the data for all samples, that it the size of the disk is
greater than $ns$. We divide the data into $m=\frac{2ns}{g}$ blocks, denoted
as $b_{i}$ for $i\in[m]$. This choice of $m$ ensures that two blocks can be
accommodated in memory simultaneously. The processing proceeds as follows:

  1. 1.

We assume that the memory holds blocks $b_{1}$ and $b_{2}$, along with the
global memoized quantities for SLIQN. All data blocks are also stored on the
disk.

  2. 2.

At iteration $t=1$, we process the block $b_{1}$ by executing the
corresponding algorithm updates on it.

  3. 3.

At any iteration $t>1$, we execute the algorithm updates on $b_{i\%n+1}$ while
concurrently storing the already processed block $b_{i}$ back into the disk
and loading the block $b_{(i+1)\%n+1}$ into memory to be processed next.

In practice, modern disks have access speeds of around 500 MBps, making
processing the bottleneck in this parallel architecture, rather than disk
access. For example, in our implementation of SLIQN with $g=1200$ MB,
$n=20,000$, $d=123$, $s=0.1117$ MB, and $m=4$, we observed that processing one
block took 7.8 seconds, while the load-store operation required only 3.8
seconds.

###  Appendix E CONVERGENCE ANALYSIS OF SIQN

In this section, we provide the convergence analysis of SIQN, as the
convergence analysis motivates the replacement of $\beta_{t}$ with
$\alpha_{\lceil\frac{t}{n}-1\rceil}$. We begin by showing that at each time
$t$, the matrix ${\bm{Q}}^{t}$ obtained after the first BFGS update satisfies
${\bm{Q}}^{t}\succeq\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})$. Using this, we
show that the updated Hessian approximation satisfies
${\bm{B}}_{i_{t}}^{t}\succeq\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t})$. These
observations are essential in order ensure that
$\sigma({\bm{Q}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$ and
$\sigma({\bm{B}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$ are well
defined.

###### Lemma E.1

For all $t\geq 1$ the following hold true:

$\displaystyle{\bm{Q}}^{t}\succeq\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}),$
(25) $\displaystyle{\bm{B}}_{i}^{t}\succeq\nabla^{2}f_{i}({\bm{z}}_{i}^{t}),$
(26)

for all $i\in[n]$.

##### Proof:

The proof follows by induction on $t$.

Base case: At $t=0$, (26) holds due to the initialization.

Induction Hypothesis (IH): Assume that (25) and (26) hold for $t=m-1$, for
some $m\geq 1$. We prove that (25) and (26) hold for $t=m$ in the Induction
step.

Induction step: Since
${\bm{B}}_{i_{m}}^{m-1}\succeq\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m-1})$, we
have the following:

$\displaystyle\big{(}1+\beta_{m}\big{)}^{2}{\bm{B}}_{i_{m}}^{m-1}\succeq\big{(}1+\beta_{m}\big{)}^{2}\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m-1})\stackrel{{\scriptstyle{\textit{Lem}}.\ref{integral_lemma}}}{{\succeq}}\big{(}1+\beta_{m}\big{)}{\bm{K}}^{m},$

where recall that
$\beta_{m}=\tfrac{M}{2}\left\|{\bm{z}}_{i_{m}}^{m}-{\bm{z}}_{i_{m}}^{m-1}\right\|_{{\bm{z}}_{i_{m}}^{m-1}}$.
Applying Lemma A.2, we obtain

$\displaystyle{\bm{Q}}^{m}=\mathrm{BFGS}(\big{(}1+\beta_{m}\big{)}^{2}{\bm{B}}_{i_{m}}^{m-1},\big{(}1+\beta_{m}\big{)}{\bm{K}}^{m},{\bm{z}}_{i_{m}}^{m}-{\bm{z}}_{i_{m}}^{m-1})\stackrel{{\scriptstyle{\textit{Lem}}.\ref{boundedness}}}{{\succeq}}\big{(}1+\beta_{m}\big{)}{\bm{K}}^{m}.$

Applying Lemma A.3 to relate ${\bm{K}}^{m}$ and
$\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m})$, we obtain

$\displaystyle{\bm{Q}}^{m}\succeq\big{(}1+\beta_{m}\big{)}{\bm{K}}^{m}\stackrel{{\scriptstyle{\textit{Lem}}.\ref{integral_lemma}}}{{\succeq}}\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m}).$

Therefore, (25) holds for $t=m$. Since
${\bm{Q}}^{m}\succeq\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m})$, applying Lemma
A.2, we obtain

$\displaystyle{\bm{B}}_{i_{m}}^{m}=\mathrm{BFGS}({\bm{Q}}^{m},\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m}),\bar{{\bm{u}}}({\bm{Q}}^{m},\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m}))\stackrel{{\scriptstyle{\textit{Lem}}.\ref{boundedness}}}{{\succeq}}\nabla^{2}f_{i_{m}}({\bm{z}}_{i_{m}}^{m}).$

Therefore, (26) holds for $t=m$. This completes the induction step. The proof
is hence complete by induction. $\square$

Key steps: We establish the convergence guarantees in three steps: Lemma E.2
establishes a one-step inequality that bounds the residual
$\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|$ in terms of the previous
residuals $\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|$ and the norm
error in the Hessian approximation
$\left\|{\bm{B}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|$, for
all $i\in[n]$. Lemma E.3 uses the result of Lemma E.2 to inductively show that
both the residual $\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|$ and the
Hessian approximation error
$\sigma({\bm{B}}_{i_{t}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$,
decrease linearly with $\lceil t/n\rceil$. Using the result of Lemma E.2 and
Lemma E.3, Lemma E.4 establishes a mean-superlinear convergence result. We
finally show in Theorem 1 that the residuals can be upper bounded by a
superlinearly convergent sequence.

We now present our one-step inequality. A similar inequality with
${\bm{B}}_{i}^{t-1}$ replaced with ${\bm{D}}_{i}^{t-1}$ also appears in Lemma
1 (for SLIQN). Since the proofs are identical, we refer Appendix G.1 directly
for the proof of Lemma E.2.

###### Lemma E.2

If Assumptions A1 and A2 hold, the sequence of iterates generated by SIQN
satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq{\tfrac{\tilde{L}\Gamma^{t-1}}{2}}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|^{2}+{\Gamma^{t-1}}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|,$
(27)

for all $t\geq 1$, where
$\Gamma^{t}:=\left\|\big{(}\sum_{i=1}^{n}{\bm{B}}_{i}^{t}\big{)}^{-1}\right\|$.

##### Proof:

Refer Appendix G.1. $\square$

###### Lemma E.3

If Assumptions A1 and A2 hold, for any $\rho$ such that
$0<\rho<1-\tfrac{\mu}{dL}$, there exist positive constants
$\epsilon^{{\textit{siqn}}}$ and $\sigma_{0}^{{\textit{siqn}}}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon^{{\textit{siqn}}}$
and
$\sigma({\bm{B}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}^{{\textit{siqn}}}$
for all $i\in[n]$, the sequence of iterates generated by SIQN satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\rho^{\lceil\tfrac{t}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$
(28)

Further, it holds that

$\displaystyle\sigma({\bm{B}}_{i_{t}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))\leq\big{(}1-\tfrac{\mu}{dL}\big{)}^{\lceil\tfrac{t}{n}\rceil}{\delta^{{\textit{siqn}}}},$
(29)

where
${\delta^{{\textit{siqn}}}}:=e^{\tfrac{4M\sqrt{L}\epsilon^{{\textit{siqn}}}}{1-\rho}}\big{(}\sigma_{0}^{{\textit{siqn}}}+\epsilon^{{\textit{siqn}}}\tfrac{4Md\sqrt{L}}{1-\tfrac{\rho}{1-\tfrac{\mu}{dL}}}\big{)},M=\tilde{L}/{\mu^{\tfrac{3}{2}}}$.

##### Proof:

For a given $\rho$ that satisfies $0<\rho<1-\frac{\mu}{dL}$, let the variables
$\epsilon^{{\textit{siqn}}},\sigma_{0}^{{\textit{siqn}}}$ be chosen to satisfy

$\displaystyle\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L{\delta^{{\textit{siqn}}}}}{\mu}\leq\frac{\rho}{1+\rho}<1.$
(30)

###### Remark 4

Indeed there exists positive constants
$\epsilon^{\textit{siqn}},\sigma_{0}^{\textit{siqn}}$ that satisfy (30).
Recall from the premise that $\delta$ is a function of
$\epsilon^{\textit{siqn}},\sigma_{0}^{\textit{siqn}}$, and we can define the
left-hand-side of (30) as a function
$g(\epsilon^{{\textit{siqn}}},\sigma^{{\textit{siqn}}}_{0})$ as

$g(\epsilon^{{\textit{siqn}}},\sigma^{{\textit{siqn}}}_{0}):=\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L{\delta^{{\textit{siqn}}}}}{\mu}=\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+Le^{\frac{4M\sqrt{L}\epsilon^{{\textit{siqn}}}}{1-\rho}}\big{(}\sigma^{{\textit{siqn}}}_{0}+\epsilon^{{\textit{siqn}}}\frac{4Md\sqrt{L}}{1-\frac{\rho}{1-\frac{\mu}{dL}}}\big{)}}{\mu}.$

Fix ${\sigma_{0}^{\textit{siqn}}=\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho})>0}$. It
is easy to see that $g(\epsilon^{\textit{siqn}},\sigma_{0}^{\textit{siqn}})$
is continuous and monotonically increasing in $\epsilon^{\textit{siqn}}$.
Also, note that
$g(0,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))=\tfrac{\rho}{2(1+\rho)}$ and
$\lim_{\epsilon\rightarrow\infty}g(\epsilon,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))=\infty$.
We can therefore apply the Intermediate Value Theorem (IVT) to guarantee that
there exists $\epsilon>0$ such that
$g(\epsilon,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))\leq\tfrac{\rho}{1+\rho}$.

Base case: At $t=1$, applying Lemma E.2, we have

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|\leq\frac{{\tilde{L}}\Gamma^{0}}{2}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{0}-{\bm{x}}^{\star}\right\|^{2}+\Gamma^{0}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\left\|{\bm{z}}_{i}^{0}-{\bm{x}}^{\star}\right\|.$

Since ${\bm{B}}_{i}^{0}\succeq\nabla^{2}f_{i}({\bm{z}}_{i}^{0})$ and
$\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\preceq L{\bm{I}}$, applying Lemma B.1, we
have
$\sigma({\bm{B}}_{i}^{0},\nabla^{2}f_{i}({\bm{z}}_{i}^{0}))\geq\frac{1}{L}\left\|{\bm{B}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|$.
This gives

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|\leq
n\Gamma^{0}\big{(}\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L\sigma^{{\textit{siqn}}}_{0}\big{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$

where we have used ${\bm{z}}_{i}^{0}={\bm{x}}^{0}$ and
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon^{{\textit{siqn}}},\sigma({\bm{B}}_{i}^{0},\nabla^{2}f_{i}({\bm{z}}_{i}^{0}))\leq\sigma^{{\textit{siqn}}}_{0}$.

We now bound $\Gamma^{0}$. Define
${\bm{X}}^{0}:=\frac{1}{n}\sum_{i=1}^{n}{\bm{B}}_{i}^{0},{\bm{Y}}^{0}:=\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{0})$.
We have the following:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\stackrel{{\scriptstyle{\Delta}}}{{\geq}}\left\|{\bm{X}}^{0}-{\bm{Y}}^{0}\right\|=\left\|({\bm{Y}}^{0})(({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}})\right\|\stackrel{{\scriptstyle(a)}}{{\geq}}\mu\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|,$

where $(a)$ follows since
${\bm{Y}}^{0}=\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\succeq\mu{\bm{I}}$
($\because$ Assumption A1). This gives us

$\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|\leq\frac{L\sigma^{{\textit{siqn}}}_{0}}{\mu}\leq\frac{L{\delta^{{\textit{siqn}}}}}{\mu}\stackrel{{\scriptstyle\eqref{choice_variables}}}{{<}}\frac{\rho}{1+\rho}.$

We can now upper bound $\Gamma^{0}$ using Banach’s Lemma A.1. Consider the
matrix $({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}$ and note that it satisfies
the requirement of A.1 from the above result. Therefore, we have

$\displaystyle\left\|({\bm{X}}^{0})^{-1}{\bm{Y}}^{0}\right\|=\left\|\big{(}{\bm{I}}+(({\bm{Y}}^{0})^{-1}{\bm{X}}^{0})-{\bm{I}})\big{)}^{-1}\right\|\stackrel{{\scriptstyle{\textit{Lem}}.\ref{lem:banach}}}{{\leq}}\frac{1}{1-\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|}\leq
1+\rho.$

Recall that $\mu{\bm{I}}\preceq{\bm{Y}}^{0}$. Using this and the previous
result, we can upper bound $\left\|({\bm{X}}^{0})^{-1}\right\|$ as
$\mu\left\|({\bm{X}}^{0})^{-1}\right\|\leq\left\|({\bm{X}}^{0})^{-1}{\bm{Y}}^{0}\right\|\leq
1+\rho.$ This gives us,
$\Gamma^{0}=\frac{1}{n}\left\|({\bm{X}}^{0})^{-1}\right\|\leq\frac{1+\rho}{n\mu}$.
Therefore, we have the following bound on
$\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|$:

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|\leq\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L\sigma^{{\textit{siqn}}}_{0}}{\mu}(1+\rho)\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L{\delta^{{\textit{siqn}}}}}{\mu}(1+\rho)\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\stackrel{{\scriptstyle\eqref{choice_variables}}}{{\leq}}\rho\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$

By the updates performed by Algorithm 1, we have
${\bm{z}}_{1}^{1}={\bm{x}}^{1}$ and ${\bm{z}}_{i}^{1}={\bm{x}}^{0}$ for $i\neq
1$. Applying Lemma B.3, we obtain

$\displaystyle\sigma({\bm{B}}^{1}_{1},\nabla^{2}f_{1}({\bm{z}}_{1}^{1}))\leq\big{(}1-c\big{)}e^{\frac{4M\sqrt{L}\epsilon^{{\textit{siqn}}}}{1-\rho}}\big{(}\underbrace{\sigma({\bm{B}}_{1}^{0},\nabla^{2}f_{1}({\bm{z}}_{1}^{0}))}_{\leq\sigma^{{\textit{siqn}}}_{0}}+\epsilon^{{\textit{siqn}}}\frac{4Md\sqrt{L}}{1-\frac{\rho}{1-c}}\big{)}\leq(1-c){\delta^{{\textit{siqn}}}}.$

This completes the proof for $t=1$.

Induction Hypothesis (IH): Let (28) and (29) hold for $t\in[jn+m]$, where
$j\geq 0,0\leq m<n$.

Induction step: We then prove that (28) and (29) also hold for $t=jn+m+1$.
Recall that the tuples are updated in a deterministic cyclic order, and at the
current time $t$, we are in the $j^{\text{th}}$ cycle and have updated the
$m^{\text{th}}$ tuple. Therefore, it is easy to note that
${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn+i}$, for all $i\in[m]$, which refer to the
tuples updated in this cycle, and ${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn-n+i}$ for
all $i\in[n]\backslash[m]$, which refer to the tuples updated in the previous
cycle. From the induction hypothesis, we have

$\displaystyle\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\begin{cases}\rho^{\lceil\frac{jn+i}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[m],\\\
\rho^{\lceil\frac{(j-1)n+i)}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[n]\backslash[m].\end{cases}$
(31)

Step 1

Applying Lemma E.2 on updating ${\bm{z}}_{m+1}$, we have

$\displaystyle\left\|{\bm{z}}_{m+1}^{jn+m+1}-{\bm{x}}^{\star}\right\|$
$\displaystyle\leq\frac{\tilde{L}\Gamma^{jn+m}}{2}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|^{2}+$
$\displaystyle\hskip
56.9055pt\Gamma^{jn+m}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|,$
$\displaystyle=\frac{\tilde{L}\Gamma^{jn+m}}{2}\big{(}\sum_{i=1}^{m}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|^{2}+\sum_{i=m+1}^{n}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|^{2}\big{)}+$
$\displaystyle\hskip
56.9055pt\Gamma^{jn+m}\big{(}\sum_{i=1}^{m}\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\big{)}+$
$\displaystyle\hskip
56.9055pt\Gamma^{jn+m}\big{(}\sum_{i=m+1}^{n}\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\big{)}.$

From the induction hypothesis, we have

$\displaystyle\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|{\leq}\begin{cases}\rho^{\lceil\frac{jn+i}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[m],\\\
\rho^{\lceil\frac{(j-1)n+i)}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[n]\backslash[m],\end{cases}$
(32)

Since ${\bm{B}}_{i}^{jn+m}\succeq\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})$ (Lemma
E.1), applying Lemma B.1 and the induction hypothesis for $\sigma$, we have
the following bound on
$\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|$:

$\displaystyle\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\stackrel{{\scriptstyle{\textit{Lem}}.\ref{lem:simple_lemma}}}{{\leq}}L\sigma^{{\textit{siqn}}}({\bm{B}}_{i}^{jn+m},\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m}))\leq\begin{cases}L{\delta^{{\textit{siqn}}}}(1-c)^{\lceil\frac{jn+i}{n}\rceil}&i\in[m],\\\
L{\delta^{{\textit{siqn}}}}(1-c)^{\lceil\frac{(j-1)n+i}{n}\rceil}&i\in[n]\backslash[m].\end{cases}$
(33)

Therefore, the bound on
$\left\|{\bm{z}}_{m+1}^{jn+m+1}-{\bm{x}}^{\star}\right\|$ simplifies to

$\displaystyle\left\|{\bm{z}}_{m+1}^{jn+m+1}-{\bm{x}}^{\star}\right\|$
$\displaystyle\stackrel{{\scriptstyle\eqref{induction_hyp_1}\eqref{induction_hyp_2}}}{{\leq}}\Gamma^{jn+m}\big{(}m\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}\rho^{2j+2}+(n-m)\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}\rho^{2j}+mL{\delta^{{\textit{siqn}}}}(1-c)^{j+1}\rho^{j+1}$
(34) $\displaystyle\hskip
142.26378pt+(n-m)L{\delta^{{\textit{siqn}}}}(1-c)^{j}\rho^{j}\big{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
(35)
$\displaystyle\leq\Gamma^{jn+m}\rho^{j}(n\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}\underbrace{\rho^{j}}_{\leq
1}+mL{\delta^{{\textit{siqn}}}}\underbrace{(1-c)^{j+1}\rho}_{\leq
1}+(n-m)L{\delta^{{\textit{siqn}}}}\underbrace{(1-c)^{j})}_{\leq
1}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\leq\Gamma^{jn+m}\rho^{j}\big{(}n\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+nL{\delta^{{\textit{siqn}}}}\big{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$
(36)

We now bound $\Gamma^{jn+m}$. Define
${\bm{X}}^{jn+m}\coloneqq\tfrac{1}{n}\sum_{i=1}^{n}{\bm{B}}_{i}^{jn+m}$ and
${\bm{Y}}^{jn+m}\coloneqq\tfrac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})$.
We follow the same recipe to bound $\Gamma^{0}$ in the base case. Observe that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\stackrel{{\scriptstyle\Delta}}{{\geq}}\left\|{\bm{X}}^{jn+m}-{\bm{Y}}^{jn+m}\right\|\overset{(a)}{\geq}\mu\left\|({\bm{Y}}^{jn+m})^{-1}{\bm{X}}^{jn+m}-{\bm{I}}\right\|.$

The inequality $(a)$ follows from Assumption A1 which implies that
$\mu{\bm{I}}\preceq\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})={\bm{Y}}^{jn+m}$.
By tracking steps (35)-(36), we can establish that

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\leq
L\delta^{\textit{siqn}}.$

From the above two chain of inequalities, we deduce

$\displaystyle\|({\bm{Y}}^{jn+m})^{-1}{\bm{X}}^{jn+m}-{\bm{I}}\|\leq\tfrac{L\delta^{\textit{siqn}}}{\mu}\stackrel{{\scriptstyle\eqref{choice_variables}}}{{<}}\frac{\rho}{1+\rho}.$

We can now upper bound $\Gamma^{jn+m}$ using Banach’s Lemma by exactly
following the procedure laid out in the base case. We get that
$\Gamma^{jn+m}=\|\sum_{i=1}^{n}{\bm{B}}_{i}^{jn+m})^{-1}\|\leq\frac{1+\rho}{n\mu}$.
Therefore, we obtain

$\displaystyle\left\|{\bm{z}}_{m+1}^{jn+m+1}-{\bm{x}}^{\star}\right\|\leq\frac{\frac{\tilde{L}\epsilon^{{\textit{siqn}}}}{2}+L{\delta^{{\textit{siqn}}}}}{\mu}(1+\rho)\rho^{j}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\stackrel{{\scriptstyle\eqref{choice_variables}}}{{\leq}}\rho^{j+1}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$
(37)

Since ${\bm{z}}_{m+1}^{jn+m+1}={\bm{x}}^{jn+m+1}$, (28) holds for $t=jn+m+1$.

Step 2

Next, we prove that
$\sigma({\bm{B}}_{m+1}^{jn+m+1},\nabla^{2}f_{m+1}({\bm{z}}_{m+1}^{jn+m+1}))\leq(1-c)^{j+1}{\delta^{{\textit{siqn}}}}$.
We define the sequence $\\{{\bm{y}}_{k}\\}$, for $k=0,\dots,j+1$, such that
$\\{{\bm{y}}_{k}\\}=\\{{\bm{x}}^{0},{\bm{z}}_{m+1}^{m+1},\dots,{\bm{z}}_{m+1}^{jn+m+1}\\}$.
The sequence $\\{{\bm{y}}_{k}\\}_{k=1}^{j+1}$ comprises of the updated value
of ${\bm{z}}_{m+1}$ till the current cycle $j$. Since $\\{{\bm{y}}_{k}\\}$
comes about from the application of BFGS updates as described in the statement
of Lemma B.3, therefore $\\{{\bm{y}}_{k}\\}$ satisfies the conditions of Lemma
B.3. This implies that

$\displaystyle\left\|{\bm{y}}_{k}-{{\bm{x}}^{\star}}\right\|\leq\rho^{\lceil\frac{(k-1)n+m+1}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|=\rho^{k}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$

for $k\in[j+1]$. Since ${\bm{y}}_{j+1}={\bm{z}}_{m+1}^{jn+m+1}$, we have
proved (29) for $t=jn+m+1$. The proof is hence complete via induction.
$\square$

###### Corollary E.1

If Assumptions A1 and A2 hold, the following holds true for all $t\geq 1$:

$\displaystyle\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}\leq
2U_{t},$ (38)

where
$U_{t}:=\sqrt{L}\rho^{\lceil\frac{t}{n}-1\rceil}\epsilon^{{\textit{siqn}}}$.

##### Proof:

We can bound
$\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}$
in the following manner:

$\displaystyle\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}\stackrel{{\scriptstyle(a)}}{{\leq}}\sqrt{L}\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|$
$\displaystyle\stackrel{{\scriptstyle{\Delta}}}{{\leq}}\sqrt{L}(\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{x}}^{\star}\right\|+\left\|{\bm{z}}_{i_{t}}^{t-1}-{\bm{x}}^{\star}\right\|),$
$\displaystyle\leq\sqrt{L}\big{(}\rho^{\lceil\frac{t}{n}\rceil}+\rho^{\lceil\frac{t-n}{n}\rceil}\big{)}\epsilon^{{\textit{siqn}}}\leq
2\sqrt{L}\rho^{\lceil\frac{t}{n}-1\rceil}\epsilon^{{\textit{siqn}}},$

where $(a)$ follows since $\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})\preceq
L{\bm{I}}$ ($\because$ Assumption A1). Therefore, the correction term
$\beta_{t}=\frac{M}{2}\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}\leq
M\sqrt{L}\rho^{\lceil\frac{t}{n}-1\rceil}\epsilon^{{\textit{siqn}}}=2U_{t}$,
which establishes (38). $\square$

######  Remark 5 ($\beta_{t}$ can be bounded by a quantity that remains
constant in a cycle)

Recall that the correction factor
$\beta_{t}=\frac{M}{2}\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}$
in SIQN was introduced to ensure that
$(1+\beta_{t})^{2}{{\bm{B}}_{i_{t}}}^{t-1}\succeq(1+\beta_{t}){\bm{K}}^{t}$
(we formalized this in Lemma E.1). Intuitively, executing SIQN with a higher
correction factor
$\beta_{t}^{\textit{new}}=MU_{t}\geq\frac{M}{2}\left\|{\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}\right\|_{{\bm{z}}_{i_{t}}^{t-1}}$
(follows from Corollary E.1) which remains constant in a cycle, it should
continue to hold that
$(1+\beta_{t}^{\textit{new}})^{2}{\bm{B}}_{i_{t}}^{t-1}\succeq(1+\beta_{t}^{\textit{new}}){\bm{K}}^{t}$.
We skip the proof for the sake of brevity as it is similar to the analysis
SIQN.

Next, we present our mean superlinear convergence result for the iterates of
SIQN. The main idea behind the proof is to substitute the linear convergence
results, specifically (28) and (29) from Lemma E.3, back into the result from
Lemma E.2. By doing so, the first term on the right-hand side of (27)
converges quadratically, while the second term converges superlinearly. This
combination leads to the desired result.

###### Lemma E.4

If Assumptions A1 and A2 hold, for any $\rho$ such that
$0<\rho<1-\tfrac{\mu}{dL}$, there exist positive constants
$\epsilon^{{\textit{siqn}}}$ and $\sigma_{0}^{{\textit{siqn}}}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon^{{\textit{siqn}}}$
and
$\sigma({\bm{B}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}^{{\textit{siqn}}}$
for all $i\in[n]$, the sequence of iterates produced by SIQN satisfy

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\big{(}1-\tfrac{\mu}{dL}\big{)}^{\lceil\tfrac{t}{n}\rceil}\tfrac{1}{n}\sum_{i=1}^{n}\left\|{\bm{x}}^{t-i}-{\bm{x}}^{\star}\right\|.$

##### Proof:

Let $t=jn+m+1$, where $0\leq j$ and $0\leq m<n$. From the proof of Lemma E.3,
we have an uniform upper bound on
$\Gamma^{jn+m}=\left\|\big{(}\sum_{i=1}^{n}{\bm{B}}_{i}^{jn+m}\big{)}^{-1}\right\|$,
given by

$\displaystyle\Gamma^{jn+m}\leq\frac{1+\rho}{n\mu},$

and the following upper bound on
$\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|$:

$\displaystyle\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\stackrel{{\scriptstyle{\textit{Lem}}.\ref{lem:simple_lemma}}}{{\leq}}L\sigma({\bm{B}}_{i}^{jn+m},\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m}))\leq\begin{cases}L\delta^{{\textit{siqn}}}(1-c)^{\lceil\frac{jn+i}{n}\rceil}&i\in[m],\\\
L\delta^{{\textit{siqn}}}(1-c)^{\lceil\frac{(j-1)n+i}{n}\rceil}&i\in[n]\backslash[m].\end{cases}$

This gives
$\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\leq
L\delta^{{\textit{siqn}}}(1-c)^{j}$, for all $i\in[n]$.

Further, from LemmaE.3, we also have
$\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\rho^{j+1}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$,
for all $i\in[n]$ and
$\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\rho^{j}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$,
for all $i\in[n]\backslash[m]$. This clearly implies
$\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\rho^{j}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$,
for $i=1,\dots,n$.

Applying Lemma E.2 at $t=jn+m+1$ and upper bounding $\Gamma^{jn+m}$, we obtain

$\displaystyle\left\|{\bm{x}}^{jn+m+1}-{\bm{x}}^{\star}\right\|$
$\displaystyle\leq\frac{1+\rho}{n\mu}\sum_{i=1}^{n}\big{(}\frac{\tilde{L}}{2}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|+\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\big{)}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\frac{1+\rho}{n\mu}\big{(}\frac{\tilde{L}\epsilon}{2}+L\delta^{{\textit{siqn}}}\big{)}(1-c)^{j}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\leq\rho(1-c)^{j}\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\stackrel{{\scriptstyle(b)}}{{\leq}}(1-c)^{\lceil\frac{t}{n}\rceil}\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{x}}^{t-i}-{\bm{x}}^{\star}\right\|,$

where $(a)$ follows from the bounds
$\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\rho^{j}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$
and
$\left\|{\bm{B}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\leq
L\delta^{{\textit{siqn}}}(1-c)^{j}$ discussed above and $\rho<1-c$. Also $(b)$
follows since ${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn+i}$, for all $i\in[m]$ and
${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn-n+i}$, for all $i\in[n]\backslash[m]$, which
implies
$\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|=\sum_{i=1}^{n}\left\|{\bm{x}}^{t-i}-{\bm{x}}^{\star}\right\|$.
This completes the proof. $\square$

The mean superlinear convergence result of Lemma E.4 ultimately gives a
superlinear rate for SIQN. Note that an identical result as Lemma E.4 is given
by Theorem 1(SLIQN). Therefore, we directly provide the proof of Lemma E.4 in
Appendix G.4 while proving Theorem 1 for SLIQN.

###### Lemma E.5

If Assumptions A1 and A2 hold, for any $\rho$ such that
$0<\rho<1-\tfrac{\mu}{dL}$, there exist positive constants
$\epsilon^{{\textit{siqn}}}$ and $\sigma_{0}^{{\textit{siqn}}}$ such that if
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon^{{\textit{siqn}}}$
and
$\sigma({\bm{B}}_{i}^{0},\nabla^{2}f_{i}({\bm{x}}^{0}))\leq\sigma_{0}^{{\textit{siqn}}}$
for all $i\in[n]$, for the sequence of iterates $\\{{\bm{x}}^{t}\\}$ generated
by SIQN, there exists a sequence $\\{\zeta^{k}\\}$ such that
$\left\|{\bm{x}}^{t}-{\bm{x}}^{\star}\right\|\leq\zeta^{\lfloor\tfrac{t-1}{n}\rfloor}$
for all $t\geq 1$ and the sequence $\\{\zeta^{k}\\}$ satisfies

$\displaystyle\zeta^{k}\leq\epsilon\big{(}1-\tfrac{\mu}{dL}\big{)}^{\tfrac{(k+2)(k+1)}{2}},$
(39)

forall $k\geq 0$.

##### Proof:

Refer Appendix G.4. $\square$

###  Appendix F EFFICIENT IMPLEMENTATION OF SLIQN

We begin by showing that the update of ${\bm{x}}^{t}$ in the SLIQN algorithm
(Algorithm 2) can be carried out in $\mathcal{O}(d^{2})$ cost.

####  F.1 Carrying out (10) in Algorithm 2 in $\mathcal{O}(d^{2})$ cost

We begin by defining the following variables that track the summands in (10)

$\displaystyle\bar{{\bm{D}}}^{t}:=\sum_{i=1}^{n}{\bm{D}}_{i}^{t},\bm{\phi}^{t}:=\sum_{i=1}^{n}{\bm{D}}_{i}^{t}{\bm{z}}_{i}^{t},{\bm{g}}^{t}:=\sum_{i=1}^{n}\nabla
f_{i}({\bm{z}}_{i}^{t}).$ (40)

The update (10) (for time $t+1$) can be expressed in terms of the defined
variables as

$\displaystyle{\bm{x}}^{t+1}=\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}\big{(}\bm{\phi}^{t}-{\bm{g}}^{t}\big{)}.$
(41)

From the updates performed by Algorithm 2 at time $t$, we have

$\displaystyle{\bm{g}}^{t}={\bm{g}}^{t-1}+\big{(}\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})\big{)}.$ (42)

Further, we can express $\bm{\phi}^{t}$ in terms of $\bm{\phi}^{t-1}$ as
follows:

$\displaystyle\bm{\phi}^{t}=\sum_{i=1,i\neq
i_{t}}^{n}\omega_{t}{\bm{D}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}+{\bm{D}}_{i_{t}}^{t}{\bm{z}}_{i_{t}}^{t}=\omega_{t}\bm{\phi}^{t-1}-\omega_{t}{\bm{D}}_{i_{t}}^{t-1}{\bm{z}}_{i_{t}}^{t-1}+{\bm{D}}_{i_{t}}^{t}{\bm{z}}_{i_{t}}^{t}.$
(43)

This updating scheme can be implemented iteratively, where we only evaluate
$\bm{\phi}^{0},{\bm{g}}^{0}$ explicitly and evaluate $\bm{\phi}^{t}$ in
$\mathcal{O}(d^{2})$ cost, for all $t\geq 1$ by (43) and ${\bm{g}}^{t}$ in
$\mathcal{O}(d)$ cost, for all $t\geq 1$ by (42).

Next, we demonstrate the method for updating $(\bar{{\bm{D}}}^{t})^{-1}$. We
begin by expressing $\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}$ in terms of
$\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}$ in the following manner:

$\displaystyle(\bar{{\bm{D}}}^{t})^{-1}$
$\displaystyle=\bigg{(}\sum_{i=1,i\neq
i_{t}}^{n-1}\omega_{t}{\bm{D}}_{i}^{t-1}+{\bm{D}}_{i_{t}}^{t}\bigg{)}^{-1}=\big{(}\omega_{t}\big{(}\bar{{\bm{D}}}^{t-1}-{\bm{D}}_{i_{t}}^{t-1}\big{)}+{\bm{D}}_{i_{t}}^{t}\big{)}^{-1},$
$\displaystyle=\omega_{t}^{-1}\big{(}\bar{{\bm{D}}}^{t-1}+\omega_{t}^{-1}{\bm{D}}_{i_{t}}^{t}-{\bm{D}}_{i_{t}}^{t-1}\big{)}^{-1}.$
(44)

Expanding the BFGS update (3), we can express
$\omega_{t}^{-1}{\bm{D}}_{i_{t}}^{t}$ as

$\displaystyle\omega_{t}^{-1}{\bm{D}}_{i_{t}}^{t}\stackrel{{\scriptstyle\eqref{eff_second_modified_update}}}{{=}}{\bm{Q}}^{t}-\frac{{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}+\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}(\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle},$
(45)

where we have used the shorthand ${\bar{{\bm{u}}}}^{t}$ for
${{\bar{{\bm{u}}}}}({\bm{Q}}_{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}))$.
Further, ${\bm{Q}}^{t}$ can be expressed as

$\displaystyle{\bm{Q}}^{t}\stackrel{{\scriptstyle\eqref{eff_first_modified_update_}}}{{=}}{\bm{D}}_{i_{t}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle},$
(46)

where
${{\bm{s}}_{i_{t}}^{t}}={\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1},{{\bm{y}}_{i_{t}}^{t}}=(1+\alpha_{\lceil\nicefrac{{t}}{{n}}-1\rceil}){\bm{K}}^{t}{{\bm{s}}_{i_{t}}^{t}}=(1+\alpha_{\lceil\nicefrac{{t}}{{n}}-1\rceil})(\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1}))$.

Adding (45) and (46), we obtain

$\displaystyle\omega_{t}^{-1}{\bm{D}}_{i_{t}}^{t}-{\bm{D}}_{i_{t}}^{t-1}$
$\displaystyle=\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}+$
$\displaystyle\hskip
113.81102pt\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}(\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle}.$
(47)

Continuing from (44), we obtain

$\displaystyle\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}=\omega_{t}^{-1}\bigg{(}\bar{{\bm{D}}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}$
$\displaystyle-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}$
$\displaystyle+\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}(\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle}\bigg{)}^{-1}.$

Next, we define the following matrix:

$\displaystyle{\bm{\psi}}_{1}:=\bar{{\bm{D}}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{{\bm{y}}_{i_{t}}^{t}}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}.$

Expressing $\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}$ in terms of
${\bm{\psi}}_{1}^{-1}$ we get the following:

$\displaystyle\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}$
$\displaystyle=\omega_{t}^{-1}\bigg{(}{\bm{\psi}}_{1}+\frac{\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}(\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle}\bigg{)}^{-1},$
$\displaystyle\stackrel{{\scriptstyle\eqref{shermon}}}{{=}}\omega_{t}^{-1}\bigg{(}{\bm{\psi}}^{-1}-\frac{{\bm{\psi}}_{1}^{-1}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}(\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t})^{\mathsf{T}}{\bm{\psi}}_{1}^{-1}}{\langle{\bar{{\bm{u}}}}^{t},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle+{\langle\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t},{\bm{\psi}}_{1}^{-1}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}^{t}){\bar{{\bm{u}}}}^{t}\rangle}}\bigg{)}.$
(48)

Define the following matrix:

$\displaystyle{\bm{\psi}}_{2}:=\bar{{\bm{D}}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}\implies{\bm{\psi}}_{1}={\bm{\psi}}_{2}-\frac{{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}.$

Expressing ${\bm{\psi}}_{1}^{-1}$ in terms of ${\bm{\psi}}_{2}^{-1}$ we get

$\displaystyle{\bm{\psi}}_{1}^{-1}\stackrel{{\scriptstyle\eqref{shermon}}}{{=}}{\bm{\psi}}_{2}^{-1}+\frac{{\bm{\psi}}_{2}^{-1}{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}({\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t})^{\mathsf{T}}{\bm{\psi}}_{2}^{-1}}{\langle{\bar{{\bm{u}}}}^{t},{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle-\langle{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t},{\bm{\psi}}_{2}^{-1}{\bm{Q}}^{t}{\bar{{\bm{u}}}}^{t}\rangle}.$
(49)

Define the following matrix:

$\displaystyle{\bm{\psi}}_{3}:=\bar{{\bm{D}}}^{t-1}+\frac{{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle}\implies{\bm{\psi}}_{2}={\bm{\psi}}_{3}-\frac{{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}.$

Expressing ${\bm{\psi}}_{2}^{-1}$ in terms of ${\bm{\psi}}_{3}^{-1}$ we get
the following:

$\displaystyle{\bm{\psi}}_{2}^{-1}\stackrel{{\scriptstyle\eqref{shermon}}}{{=}}{\bm{\psi}}_{3}^{-1}+\frac{{\bm{\psi}}_{3}^{-1}{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}({\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}})^{\mathsf{T}}{\bm{\psi}}_{3}^{-1}}{\langle{{\bm{s}}_{i_{t}}^{t}},{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle-\langle{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}},{\bm{\psi}}_{3}^{-1}{\bm{D}}_{i_{t}}^{t-1}{{\bm{s}}_{i_{t}}^{t}}\rangle}.$
(50)

Finally, using (16) to evaluate ${\bm{\psi}}_{3}^{-1}$, we obtain

$\displaystyle{\bm{\psi}}_{3}^{-1}=\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}-\frac{\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}{{\bm{y}}_{i_{t}}^{t}}{{\bm{y}}_{i_{t}}^{t}}^{\mathsf{T}}\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}}{\langle{{\bm{y}}_{i_{t}}^{t}},{{\bm{s}}_{i_{t}}^{t}}\rangle+\langle{{\bm{y}}_{i_{t}}^{t}},\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}{{\bm{y}}_{i_{t}}^{t}}\rangle}.$
(51)

Therefore, given access to $\big{(}\bar{{\bm{D}}}^{t-1}\big{)}^{-1}$, we can
evaluate $\psi_{3}^{-1}$ in $\mathcal{O}(d^{2})$ time. Similarly, given access
to ${\bm{\psi}}_{3}^{-1}$, we can evaluate ${\bm{\psi}}_{2}^{-1}$ in
$\mathcal{O}(d^{2})$ time. Continuing similarly, we can evaluate
${\bm{\psi}}_{1}^{-1}$ and $\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1}$ in
$\mathcal{O}(d^{2})$ time. This scheme can be enumerated iteratively where we
only compute $\big{(}\bar{{\bm{D}}}^{0}\big{)}^{-1}$ explicitly and evaluate
$\big{(}\bar{{\bm{D}}}^{t}\big{)}^{-1},\forall t\geq 1$ by the steps (51),
(50), (49), and (48). Therefore, the update (10) in Algorithm 2 can be
performed in $\mathcal{O}(d^{2})$ time.

####  F.2 Efficient Implementation of Algorithm 2 in $\mathcal{O}(d^{2})$ cost

By lazily carrying out the scaling of the ${\bm{D}}$’s, i.e., only scaling
when they are used, we can improve the per-iteration complexity of SLIQN to
$\mathcal{O}(d^{2})$. The resulting algorithm is specified by the following
pseudo code:

Algorithm 3 Sharpened Lazy Incremental Quasi-Newton (SLIQN)

1: Function {Sherman-Morrison} {${\bm{A}}^{-1},{{{\bm{u}}}},{\bm{v}}$}

2: return
${\bm{A}}^{-1}-\tfrac{{\bm{A}}^{-1}{{{\bm{u}}}}{\bm{v}}^{\mathsf{T}}{\bm{A}}^{-1}}{1+\langle{\bm{v}},{\bm{A}}^{-1}{{{\bm{u}}}}\rangle}$

3: EndFunction

4: Initialize: Initialize $\\{{\bm{z}}_{i},{\bm{D}}_{i}\\}_{i=1}^{n}$ similar
to Algorithm F.2;

5: Evaluate
$\bar{{\bm{D}}}:=\big{(}\sum_{i=1}^{n}{\bm{D}}_{i}\big{)}^{-1},\bm{\phi}:=\sum_{i=1}^{n}{\bm{D}}_{i}{\bm{z}}_{i}$,
and ${\bm{g}}:=\sum_{i=1}^{n}\nabla f_{i}({\bm{z}}_{i})$;

6: Maintain running auxiliary variables
${\bm{x}},{\bar{{\bm{u}}}},{\bm{y}},{\bm{s}},{\bm{Q}},{\bm{D}}^{\text{old}},{\bm{K}}$;
// ${\bm{x}}$ keeps a track of ${\bm{x}}^{t}$, $\bar{{\bm{u}}}$ keeps a track
of the greedy vector, ${\bm{y}}$ keeps a track of $\nabla
f_{i_{t}}({\bm{z}}_{i_{t}}^{t})-\nabla f_{i_{t}}({\bm{z}}_{i_{t}}^{t-1})$,
${\bm{s}}$ keeps a track of ${\bm{z}}_{i_{t}}^{t}-{\bm{z}}_{i_{t}}^{t-1}$,
whereas ${\bm{Q}},{\bm{D}}^{\text{old}},{\bm{K}}$ keep track of the
intermediate matrices

7: while not converged:

8: Current index to be updated is $i_{t}\leftarrow(t-1)\mod n+1;$

9: Update ${\bm{x}}$ as
${\bm{x}}\leftarrow\big{(}\bar{{\bm{D}}}\big{)}\big{(}\bm{\phi}-{\bm{g}}\big{)};$

10: Update ${\bm{s}}\leftarrow{\bm{x}}-{\bm{z}}_{i_{t}};$

11: Update ${\bm{y}}\leftarrow\nabla f_{i_{t}}({\bm{x}}^{t})-\nabla
f_{i_{t}}({\bm{z}}_{i_{t}});$

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

13: Update ${\bm{Q}}$ as ${\bm{Q}}\leftarrow\mathrm{BFGS}((1+\alpha_{\lceil
t/n-1\rceil})^{2}{\bm{D}}_{i_{t}},(1+\alpha_{\lceil
t/n-1\rceil}){\bm{K}},{\bm{s}})$, where // Lazy Step
${\bm{K}}\leftarrow\int_{0}^{1}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}+\tau({\bm{x}}^{t}-{\bm{z}}_{i_{t}}))d\tau.$

14: Update ${\bar{{\bm{u}}}}$ as
${\bar{{\bm{u}}}}\leftarrow\operatorname*{arg\,max}_{{{{\bm{u}}}}\in\\{{\bm{e}}_{i}\\}_{i=1}^{d}}\frac{\langle{{{\bm{u}}}},{\bm{Q}}{{{\bm{u}}}}\rangle}{\langle{{{\bm{u}}}},\nabla^{2}f_{i_{t}}({\bm{x}}^{t}){{{\bm{u}}}}\rangle}$;

15: Update ${\bm{D}}^{\text{old}}$ as
${\bm{D}}^{\text{old}}\leftarrow{\bm{D}}_{i_{t}}$;

16: Update ${\bm{D}}_{i_{t}}$ as
${\bm{D}}_{i_{t}}\leftarrow\mathrm{BFGS}({\bm{Q}},\nabla^{2}f_{i_{t}}({\bm{x}}^{t}),{\bar{{\bm{u}}}})$;

17: Update $\bm{\phi}$ as
$\bm{\phi}\leftarrow\omega_{t}(\bm{\phi}-{\bm{D}}^{\text{old}}{\bm{z}}_{i_{t}})+{\bm{D}}_{i_{t}}{\bm{x}}^{t}$;

18: Update ${\bm{g}}$ as ${\bm{g}}\leftarrow{\bm{g}}+\big{(}\nabla
f_{i_{t}}({\bm{x}}^{t})-\nabla f_{i_{t}}({\bm{z}}_{i_{t}})\big{)}$;

19: Update $\bar{{\bm{D}}}$ as $\bar{{\bm{D}}}\leftarrow\text{Sherman-
Morrison}(\bar{{\bm{D}}},{\bm{y}},\frac{1}{\langle{\bm{y}},{\bm{s}}\rangle}{\bm{y}})$;

20: Update $\bar{{\bm{D}}}$ as $\bar{{\bm{D}}}\leftarrow\text{Sherman-
Morrison}(\bar{{\bm{D}}},-{\bm{D}}^{\text{old}}{\bm{s}},\frac{1}{\langle{\bm{s}},{\bm{D}}^{\text{old}}{\bm{s}}\rangle}{\bm{D}}^{\text{old}}{\bm{s}})$;

21: Update $\bar{{\bm{D}}}$ as $\bar{{\bm{D}}}\leftarrow\text{Sherman-
Morrison}(\bar{{\bm{D}}},-{\bm{Q}}{\bar{{\bm{u}}}},\frac{1}{\langle{\bar{{\bm{u}}}},{\bm{Q}}{\bar{{\bm{u}}}}\rangle}{\bm{Q}}{\bar{{\bm{u}}}})$;

22: Update $\bar{{\bm{D}}}$ as
$\bar{{\bm{D}}}\leftarrow\omega_{t}^{-1}\text{Sherman-
Morrison}(\bar{{\bm{D}}},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}){\bar{{\bm{u}}}},\frac{1}{\langle{\bar{{\bm{u}}}},\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}){\bar{{\bm{u}}}}\rangle}\nabla^{2}f_{i_{t}}({\bm{z}}_{i_{t}}){\bar{{\bm{u}}}})$;

23: Update ${\bm{z}}_{i_{t}}$ as ${\bm{z}}_{i_{t}}\leftarrow{\bm{x}}$;

24: Increment the iteration counter $t$;

25: end while

###  Appendix G CONVERGENCE ANALYSIS OF SLIQN

####  G.1 Proof of Lemma 1

For all $t\geq 0$, we define
${\bm{H}}^{t}:=\big{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t}\big{)}^{-1}$. From the
update for ${\bm{x}}^{t}$ (8), we have

$\displaystyle{\bm{x}}^{t}-{\bm{x}}^{*}$
$\displaystyle={\bm{H}}^{t-1}\bigg{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}{\bm{z}}_{i}^{t-1}-\sum_{i=1}^{n}\nabla
f_{i}({\bm{z}}_{i}^{t-1})\bigg{)}-{\bm{x}}^{*},$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}{\bm{H}}^{t-1}\bigg{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})-\sum_{i=1}^{n}\nabla
f_{i}({\bm{z}}_{i}^{t-1})\bigg{)},$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{=}}{\bm{H}}^{t-1}\bigg{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})-\sum_{i=1}^{n}(\nabla
f_{i}({\bm{z}}_{i}^{t-1})-\nabla f_{i}({\bm{x}}^{\star}))\bigg{)},$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}{\bm{H}}^{t-1}\bigg{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})-\sum_{i=1}^{n}\int_{0}^{1}\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v)({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})dv\bigg{)},$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}{\bm{H}}^{t-1}\bigg{(}\sum_{i=1}^{n}\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})+$
$\displaystyle\hskip
85.35826pt\sum_{i=1}^{n}\int_{0}^{1}(\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})-\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v))({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})dv\bigg{)}.$

The equality $(a)$ follows from the definition
${\bm{H}}^{t-1}=\big{(}\sum_{i=1}^{n}{\bm{D}}_{i}^{t-1}\big{)}^{-1}$. The
equality $(b)$ uses the fact that $\nabla
f({\bm{x}}^{*})=\frac{1}{n}\sum_{i=1}^{n}\nabla
f_{i}({\bm{x}}^{\star})=\mathbf{0}$. The equality (c) follows from the
Fundamental Theorem of Calculus, and the equality $(d)$ follows by adding and
subtracting
$\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})$.
Taking norm on both sides and applying the Triangle inequality, we obtain

$\displaystyle\left\|{\bm{x}}^{t}-{\bm{x}}^{*}\right\|$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{\leq}}\left\|{\bm{H}}^{t-1}\right\|\bigg{(}\sum_{i=1}^{n}\left\|\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|+$
$\displaystyle\hskip
56.9055pt\sum_{i=1}^{n}\left\|\int_{0}^{1}(\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})-\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v))({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})dv\right\|\bigg{)},$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\left\|{\bm{H}}^{t-1}\right\|\bigg{(}\sum_{i=1}^{n}\left\|\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|+$
$\displaystyle\sum_{i=1}^{n}\int_{0}^{1}\left\|(\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})-\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v))({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})\right\|dv\bigg{)},$
$\displaystyle\leq\left\|{\bm{H}}^{t-1}\right\|\bigg{(}\sum_{i=1}^{n}\left\|\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|+$
$\displaystyle\sum_{i=1}^{n}\int_{0}^{1}\left\|(\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})-\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v))\right\|\left\|({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})\right\|dv\bigg{)},$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\left\|{\bm{H}}^{t-1}\right\|\bigg{(}\sum_{i=1}^{n}\left\|\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|+\tilde{L}\int_{0}^{1}(1-v)dv\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{*}\right\|^{2}\bigg{)},$
$\displaystyle\leq\left\|{\bm{H}}^{t-1}\right\|\bigg{(}\sum_{i=1}^{n}\left\|\big{(}{\bm{D}}_{i}^{t-1}-\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})\big{)}\right\|\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|+\frac{\tilde{L}}{2}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{*}\right\|^{2}\bigg{)}.$

The inequality $(a)$ follows from the known result that if
$g\colon\mathbb{R}\rightarrow\mathbb{R}^{d}$ is a continuous function, then
$\left\|\int_{0}^{1}g(v)\,dv\right\|\leq\int_{0}^{1}\left\|g(v)\right\|\,dv$,
and the inequality $(b)$ follows from the assumption that the Hessian of
$f_{i}$ is $\tilde{L}-$ Lipschitz (A2), i.e.
$\left\|\nabla^{2}f_{i}({\bm{z}}_{i}^{t-1})-\nabla^{2}f({\bm{x}}^{*}+({\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star})v)\right\|\leq\tilde{L}(1-v)\left\|{\bm{z}}_{i}^{t-1}-{\bm{x}}^{\star}\right\|$.
This completes the proof.

####  G.2 Proof of Lemma 2

For a given $\rho$ that satisfies $0<\rho<1-\frac{\mu}{dL}$, we choose
$\epsilon,\delta>0$ such that they satisfy

$\displaystyle\tfrac{1}{\mu}\left(\tfrac{\tilde{L}\epsilon}{2}+L\delta(1+M\sqrt{L}\epsilon)^{2}+L^{\tfrac{3}{2}}M\epsilon(2+M\sqrt{L}\epsilon)\right)\leq\frac{\rho}{1+\rho}<1.$
(52)

###### Remark 6

Indeed, there exists positive constants $\epsilon,\delta>0$ that satisfy 52.
Recall from the premise that $\delta$ is a function of $\epsilon,\sigma_{0}$,
and we can define the left-hand-side of 52 as a function
$h(\epsilon,\sigma_{0})$ as

$\displaystyle
h(\epsilon,\sigma_{0})\coloneqq\tfrac{1}{\mu}\left(\tfrac{\tilde{L}\epsilon}{2}+Le^{\frac{4M\sqrt{L}\epsilon}{1-\rho}}\big{(}\sigma_{0}+\epsilon\tfrac{4Md\sqrt{L}}{1-\rho(1-\tfrac{\mu}{dL})^{-1}}\big{)}(1+M\sqrt{L}\epsilon)^{2}+L^{\tfrac{3}{2}}M\epsilon(2+M\sqrt{L}\epsilon)\right).$

Fix ${\sigma_{0}=\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho})>0}$. It is easy to see
that $h(\epsilon,\sigma_{0})$ is continuous and monotonically increasing in
$\epsilon$. Also, note that
$h(0,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))=\tfrac{\rho}{2(1+\rho)}$ and
$\lim_{\epsilon\rightarrow\infty}h(\epsilon,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))=\infty$.
We can therefore apply the Intermediate Value Theorem (IVT) to guarantee that
there exists $\epsilon>0$ such that
$h(\epsilon,\tfrac{\mu}{2L}(\tfrac{\rho}{1+\rho}))\leq\tfrac{\rho}{1+\rho}$.

Similar to Lemma E.3, we use Induction on $t$ to prove the result.

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

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|\leq\frac{\tilde{L}\Gamma^{0}}{2}\sum_{i=1}^{n}\left\|{\bm{z}}_{i}^{0}-{\bm{x}}^{\star}\right\|^{2}+\Gamma^{0}\sum_{i=1}^{n}\left\|{\bm{D}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\left\|{\bm{z}}_{i}^{0}-{\bm{x}}^{\star}\right\|.$

From the initialization, we have that
${\bm{D}}_{i}^{0}=(1+\alpha_{0})^{2}{\bm{I}}_{i}^{0}$ and
${\bm{z}}_{i}^{0}={\bm{x}}^{0}$, for all $i\in[n]$, and
$\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\leq\epsilon$. Substituting these
in the above expression, we obtain

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|$
$\displaystyle\leq\Gamma^{0}\bigg{(}n\tfrac{\tilde{L}\epsilon}{2}+\sum_{i=1}^{n}\left\|(1+\alpha_{0})^{2}{\bm{I}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\bigg{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
(53)
$\displaystyle\overset{(a)}{\leq}\Gamma^{0}\bigg{(}n\tfrac{\tilde{L}\epsilon}{2}+(1+\alpha_{0})^{2}\sum_{i=1}^{n}\left\|{\bm{I}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|+\alpha_{0}(\alpha_{0}+2)\sum_{i=1}^{n}\underbrace{\left\|\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|}_{\leq
L}\bigg{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\overset{(b)}{\leq}\Gamma^{0}\
\bigg{(}n\tfrac{\tilde{L}\epsilon}{2}+nL(1+\alpha_{0})^{2}\sigma_{0}+nL\alpha_{0}(\alpha_{0}+2)\bigg{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\leq\Gamma^{0}\bigg{(}n\tfrac{\tilde{L}\epsilon}{2}+nL(1+M\sqrt{L}\epsilon)^{2}\sigma_{0}+nL^{\frac{3}{2}}M\epsilon(M\sqrt{L}\epsilon+2)\bigg{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
$\displaystyle\leq\Gamma^{0}\bigg{(}n\tfrac{\tilde{L}\epsilon}{2}+nL(1+M\sqrt{L}\epsilon)^{2}\delta+nL^{\frac{3}{2}}M\epsilon(M\sqrt{L}\epsilon+2)\bigg{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|,$
(54)

where $(a)$ follows by adding and subtracting
$(1+\alpha_{0})^{2}\nabla^{2}f_{i}({\bm{z}}_{i}^{0})$ to
$(1+\alpha_{0})^{2}{\bm{I}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})$ and
applying the Triangle inequality. To see why inequality $(b)$ is true, first
recall from the initialization that
$\sigma({\bm{I}}_{i}^{0},\nabla^{2}f_{i}({\bm{z}}_{i}^{0}))\leq\sigma_{0}$ and
${\bm{I}}_{i}^{0}\succeq\nabla^{2}f_{i}({\bm{z}}_{i}^{0})$. Applying Lemma
B.1, we have
$\left\|{\bm{I}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\leq
L\sigma({\bm{I}}_{i}^{0},\nabla^{2}f_{i}({\bm{z}}_{i}^{0}))\leq L\sigma_{0}$.

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

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\|{\bm{D}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\stackrel{{\scriptstyle\Delta}}{{\geq}}\left\|{\bm{X}}^{0}-{\bm{Y}}^{0}\right\|=\left\|({\bm{Y}}^{0})(({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}})\right\|\stackrel{{\scriptstyle(a)}}{{\geq}}\mu\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|,$

where the inequality $(a)$ follows from Assumption A1 which implies
$\mu{\bm{I}}\preceq\frac{1}{n}\sum_{i=1}^{n}\nabla^{2}f_{i}({\bm{z}}_{i}^{0})={\bm{Y}}^{0}$.
By tracking the parts of steps (53)-(54) which bound
$\left\|{\bm{D}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|$, we get

$\displaystyle\sum_{i=1}^{n}\left\|{\bm{D}}_{i}^{0}-\nabla^{2}f_{i}({\bm{z}}_{i}^{0})\right\|\leq
nL(1+M\sqrt{L}\epsilon)^{2}\delta+nL^{\frac{3}{2}}M\epsilon(M\sqrt{L}\epsilon+2).$

From the above two chain of inequalities, we obtain

$\displaystyle\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|\leq\tfrac{1}{\mu}\bigg{(}L(1+M\sqrt{L}\epsilon)^{2}\delta+L^{\frac{3}{2}}M\epsilon(M\sqrt{L}\epsilon+2)\bigg{)}\stackrel{{\scriptstyle\eqref{choice_variables_II}}}{{<}}\frac{\rho}{1+\rho}.$

We can now upper bound $\Gamma^{0}$ using Banach’s Lemma A.1. Consider the
matrix $({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}$ and note that it satisfies
the requirement of A.1 from the above result. Therefore, we have

$\displaystyle\left\|({\bm{X}}^{0})^{-1}{\bm{Y}}^{0}\right\|=\left\|\big{(}{\bm{I}}+(({\bm{Y}}^{0})^{-1}{\bm{X}}^{0})-{\bm{I}})\big{)}^{-1}\right\|\stackrel{{\scriptstyle{\textit{Lem}}.\ref{lem:banach}}}{{\leq}}\frac{1}{1-\left\|({\bm{Y}}^{0})^{-1}{\bm{X}}^{0}-{\bm{I}}\right\|}\leq
1+\rho.$

Recall that $\mu{\bm{I}}\preceq{\bm{Y}}^{0}$. Using this and the previous
result, we can upper bound $\left\|({\bm{X}}^{0})^{-1}\right\|$ as
$\mu\left\|({\bm{X}}^{0})^{-1}\right\|\leq\left\|({\bm{X}}^{0})^{-1}{\bm{Y}}^{0}\right\|\leq
1+\rho.$ This gives us,
$\Gamma^{0}=\frac{1}{n}\left\|({\bm{X}}^{0})^{-1}\right\|\leq\frac{1+\rho}{n\mu}$.

Substituting this bound on $\Gamma^{0}$ in 54, we obtain

$\displaystyle\left\|{\bm{x}}^{1}-{\bm{x}}^{\star}\right\|\leq\tfrac{1+\rho}{\mu}\big{(}\tfrac{\tilde{L}\epsilon}{2}+L(1+M\sqrt{L}\epsilon)^{2}\delta+L^{\frac{3}{2}}M\epsilon(M\sqrt{L}\epsilon+2)\big{)}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|\stackrel{{\scriptstyle\eqref{choice_variables_II}}}{{\leq}}\rho\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|.$

To complete the base step, we now upper bound
$\sigma(\omega_{1}^{-1}{\bm{D}}^{1}_{1},\nabla^{2}f_{1}({\bm{z}}_{1}^{1}))$,
where $\omega_{1}=1$ . Applying Lemma B.3 with parameters as
$T=1,\tilde{{\bm{x}}}={\bm{x}}^{\star},{\bm{P}}^{0}=(1+\alpha_{0})^{2}{\bm{I}}_{1}^{0}={\bm{D}}_{1}^{0}$
(refer to Remark 3 we made for Lemma B.3, where we stated that the results of
Lemma B.3 remain unchanged on redefining
$r_{k}:=2\sqrt{L}\rho^{k-1}\epsilon\geq
2\sqrt{L}\rho^{k-1}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|$), we get

$\displaystyle\sigma({\bm{D}}^{1}_{1},\nabla^{2}f_{1}({\bm{z}}_{1}^{1}))\leq\big{(}1-c\big{)}e^{\frac{4M\sqrt{L}\epsilon}{1-\rho}}\big{(}\underbrace{\sigma({\bm{I}}_{1}^{0},\nabla^{2}f_{1}({\bm{z}}_{1}^{0}))}_{\leq\sigma_{0}}+\epsilon\frac{4Md\sqrt{L}}{1-\frac{\rho}{1-c}}\big{)}\leq(1-c)\delta.$

Finally, as a technical remark, the proof of Lemma B.3 already shows that
${\bm{D}}^{1}_{1}\succeq\nabla^{2}f_{1}({\bm{z}}_{1}^{1})$, and therefore
$\sigma({\bm{D}}^{1}_{1},\nabla^{2}f_{1}({\bm{z}}_{1}^{1}))$ is well defined.
This completes the proof for for the base case.

We now prove that (12) and (13) hold for any $t>1$ by induction.

Induction hypothesis (IH): Let (12) and (13) hold for all $t\in[jn+m]$ for
some $j\geq 0$ and $0\leq m<n$.

Induction step: We then prove that (12) and (13) also hold for $t=jn+m+1$.
Recall that the tuples are updated in a deterministic cyclic order, and at the
current time $t$, we are in the $j^{\text{th}}$ cycle and have updated the
$m^{\text{th}}$ tuple. Therefore, it is easy to note that
${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn+i}$, for all $i\in[m]$, which refer to the
tuples updated in this cycle, and ${\bm{z}}_{i}^{jn+m}={\bm{x}}^{jn-n+i}$ for
all $i\in[n]\backslash[m]$, which refer to the tuples updated in the previous
cycle. From the induction hypothesis, we have

$\displaystyle\left\|{\bm{z}}_{i}^{jn+m}-{\bm{x}}^{\star}\right\|\leq\begin{cases}\rho^{\lceil\frac{jn+i}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[m],\\\
\rho^{\lceil\frac{(j-1)n+i)}{n}\rceil}\left\|{\bm{x}}^{0}-{\bm{x}}^{\star}\right\|&i\in[n]\backslash[m].\end{cases}$
(55)

We will execute the induction step in three distinct stages. In the first
stage we will establish an upper bound on
$\sum_{i=1}^{n}\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\|$.
In the second stage, we will use the previous result and Lemma 1 to bound
$\|{\bm{z}}_{m+1}^{jn+m+1}-{\bm{x}}^{\star}\|$. In the final stage, we will
prove the linear convergence of the updated Hessian approximation, i.e.,
$\sigma(\omega_{jn+m+1}^{-1}{\bm{D}}_{m+1}^{jn+m+1},\nabla^{2}f_{m+1}({\bm{z}}_{m+1}^{jn+m+1}))\leq(1-c)^{j+1}\delta$.

Since ${\bm{D}}_{i_{t}}^{t}$ is updated in a different manner for $t\bmod
n\neq 0$ and $t\bmod n=0$, we split the first stage into two cases
corresponding to $t\bmod n\neq 0$ and $t\bmod n=0$.

Stage 1, Case 1: ${t\bmod 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 the pre-multiplication of the scaling
factor, it is easy to note that ${\bm{D}}_{i}^{jn+m}={\bm{D}}_{i}^{jn+i}$, for
all $i\in[m]$,
${\bm{D}}_{i}^{jn+m}=(1+M\sqrt{L}\epsilon\rho^{j})^{2}{\bm{D}}_{i}^{jn-n+i}$,
for all $i\in[n-1]\backslash[m]$, and ${\bm{D}}_{i}^{jn+m}={\bm{D}}_{i}^{jn}$,
for $i=n$.

We want to bound
$\sum_{i=1}^{n}\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\|$.
For all $i\in[m]$, from the induction hypothesis, we have

$\displaystyle\sigma({\bm{D}}_{i}^{jn+m},\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m}))\leq(1-c)^{\lceil\frac{jn+i}{n}\rceil}\delta,$
$\displaystyle\overset{{\textit{Lem}}.{\ref{lem:simple_lemma}}}{\implies}$
$\displaystyle\left\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|=\left\|{\bm{D}}_{i}^{jn+i}-\nabla^{2}f_{i}({\bm{x}}^{jn+i})\right\|\leq
L(1-c)^{\lceil\frac{jn+i}{n}\rceil}\delta.$ (56)

For all $i\in[n-1]\backslash[m]$, we follow in the footsteps of 53-54 from the
base case to get

$\displaystyle\left\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|$
$\displaystyle=\left\|(1+M\sqrt{L}\epsilon\rho^{j})^{2}{\bm{D}}_{i}^{jn-n+i}-\nabla^{2}f_{i}({\bm{x}}^{jn-n+i})\right\|,$
$\hskip
22.76219pt\leq(1+M\sqrt{L}\epsilon\rho^{j})^{2}\left\|{\bm{D}}_{i}^{jn-n+i}-\nabla^{2}f_{i}({\bm{x}}^{jn-n+i})\right\|+\\\
M\sqrt{L}\epsilon\rho^{j}(2+M\sqrt{L}\epsilon\rho^{j})\underbrace{\left\|\nabla^{2}f_{i}({\bm{x}}^{jn-n+i})\right\|}_{\leq
L},$
$\displaystyle\leq(1+M\sqrt{L}\epsilon\rho^{j})^{2}\left\|{\bm{D}}_{i}^{jn-n+i}-\nabla^{2}f_{i}({\bm{x}}^{jn-n+i})\right\|+ML^{\frac{3}{2}}\epsilon\rho^{j}(2+M\sqrt{L}\epsilon\rho^{j}),$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}(1+M\sqrt{L}\epsilon\rho^{j})^{2}L\delta(1-c)^{j}+ML^{\frac{3}{2}}\epsilon\rho^{j}(2+M\sqrt{L}\epsilon\rho^{j}).$
(57)

The inequality $(a)$ follows from
$\|{\bm{D}}_{i}^{jn-n+i}-\nabla^{2}f_{i}({\bm{x}}^{jn-n+i})\|\leq(1-c)^{j}\delta$,
which can be established from the induction hypothesis in a similar way as we
did for the case with $i\in[m]$. Next, for $i=n$, we have

$\displaystyle\left\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|$
$\displaystyle=\left\|{\bm{D}}_{i}^{jn}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\right\|,$
$\displaystyle=\left\|\omega_{jn}(\omega_{jn}^{-1}{\bm{D}}_{i}^{jn}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn}))+(\omega_{jn}-1)\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|,$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{\leq}}\omega_{jn}\|\omega_{jn}^{-1}{\bm{D}}_{i}^{jn}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\|+\lvert\omega_{jn}-1\rvert\|\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\|,$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}(1+M\sqrt{L}\epsilon\rho^{j})L\delta(1-c)^{j}+\lvert\omega_{jn}-1\rvert\|\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\|,$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}(1+M\sqrt{L}\epsilon\rho^{j})L\delta(1-c)^{j}+ML^{\frac{3}{2}}\epsilon\rho^{j}(2+M\sqrt{L}\epsilon\rho^{j}).$
(58)

To see the deduction $(a)$, we follow in the footsteps of the case with
$i\in[m]$. Concretely, from induction and Lemma B.1,
$\|\omega_{jn}^{-1}{\bm{D}}_{i}^{jn}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\|\leq
L\sigma(\omega_{jn}^{-1}{\bm{D}}_{i}^{jn},\nabla^{2}f_{i}({\bm{z}}_{i}^{jn}))\leq
L\delta(1-c)^{j}$. Inequality $(b)$ uses the fact that
$\omega_{jn}=1+M\sqrt{L}\epsilon\rho^{j}$ and
$\|\nabla^{2}f_{i}({\bm{z}}_{i}^{jn})\|\leq L$ ($\because$ Assumption A1).

We can now bound the quantity
$\sum_{i=1}^{n}\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\|$
using(56), (57) and (58), as follows:

$\sum_{i=1}^{n}\left\|{\bm{D}}_{i}^{jn+m}-\nabla^{2}f_{i}({\bm{z}}_{i}^{jn+m})\right\|\leq
mL\delta(1-c)^{j+1}+(n-m)\big{(}(1+M\sqrt{L}\epsilon\rho^{j})^{2}L\delta(1-c)^{j}\\\
+ML^{\frac{3}{2}}\epsilon\rho^{j}(2+M\sqrt{L}\epsilon\rho^{j})\big{)},$