Stochastic L-BFGS: Improved Convergence Rates and Practical Acceleration Strategies

TL;DR

This paper introduces a new convergence analysis framework for stochastic L-BFGS, demonstrating improved theoretical rates and practical acceleration strategies that significantly enhance empirical performance on large-scale problems.

Contribution

The paper provides a novel convergence analysis framework for stochastic L-BFGS and proposes practical acceleration strategies with theoretical validation.

Findings

Improved convergence rates and computational complexities for stochastic L-BFGS.

Practical acceleration strategies significantly speed up empirical performance.

Experimental results show superior performance on large-scale logistic and ridge regression tasks.

Abstract

We revisit the stochastic limited-memory BFGS (L-BFGS) algorithm. By proposing a new framework for the convergence analysis, we prove improved convergence rates and computational complexities of the stochastic L-BFGS algorithms compared to previous works. In addition, we propose several practical acceleration strategies to speed up the empirical performance of such algorithms. We also provide theoretical analyses for most of the strategies. Experiments on large-scale logistic and ridge regression problems demonstrate that our proposed strategies yield significant improvements vis-\`a-vis competing state-of-the-art algorithms.

Tables

Table 1. TABLE I: Statistics of rcv1 and E2006-tf datasets.

Datasets	$d$	$n$	$d/n$	$\varrho$	$\kappa$
rcv1	47236	20242	2.33	0.0154	113.17
E2006-tf	15036	16087	0.93	0.0404	1.70

Equations

$$\underset{\mathbf{x}\in {\mathbb{R}}^d}{\min }\left[f(\mathbf{x})\triangleq \frac{1}{n}\sum _{i=1}^n{f}_i(\mathbf{x})\right],$$

$$f_i(\mathbf{x})\triangleq \ell ({\mathbf{a}}_i^T\mathbf{x},b_i)+\lambda R(\mathbf{x}).$$

$$\nabla f_{{\mathcal{B}}_{s,t}}(\mathbf{x})\triangleq \frac{1}{b}\sum _{i\in {\mathcal{B}}_{s,t}}\frac{1}{n{p}_i}\nabla f_i(\mathbf{x}),\forall \mathbf{x}\in {\mathbb{R}}^d.$$

$${\nabla }^2{f}_{{\mathcal{T}}_r}({\overline{\mathbf{x}}}_r)\triangleq \frac{1}{b_{\mathrm{H}}}\sum _{i\in {\mathcal{T}}_r}{\nabla }^2{f}_i({\overline{\mathbf{x}}}_r),$$

$${\mathbf{H}}_r^{(r-M^{\prime })}\triangleq \frac{{\mathbf{s}}_r^T{\mathbf{y}}_r}{{\Vert {\mathbf{y}}_r\Vert }^2}\mathbf{I}$$

$${\mathbf{H}}_r^{(k)}=\left(\mathbf{I}-\frac{{\mathbf{y}}_k{{\mathbf{s}}_k}^T}{{{\mathbf{y}}_k}^T{\mathbf{s}}_k}\right){\mathbf{H}}_r^{(k-1)}\left(\mathbf{I}-\frac{{\mathbf{s}}_k{{\mathbf{y}}_k}^T}{{{\mathbf{y}}_k}^T{\mathbf{s}}_k}\right)+\frac{{\mathbf{s}}_k{{\mathbf{s}}_k}^T}{{{\mathbf{y}}_k}^T{\mathbf{s}}_k}.$$

$${\mathbf{v}}_{s,t}:=\nabla f_{{\mathcal{B}}_{s,t}}({\mathbf{x}}_{s,t})-\nabla f_{{\mathcal{B}}_{s,t}}({\mathbf{x}}^s)+{\mathbf{g}}_s$$

$${\overline{\mu }}_{\tilde{n}}\triangleq \frac{1}{\tilde{n}}\sum _{j=1}^{\tilde{n}}{\mu }_{{\underline{i}}_j}\text{mbox}and{\overline{L}}_{\tilde{n}}\triangleq \frac{1}{\tilde{n}}\sum _{j=n-\tilde{n}+1}^n{L}_{{\overline{i}}_j}.$$

$${\kappa }_{\max }\triangleq \frac{L_{\max }}{{\mu }_{\min }}\text{mbox}and{\kappa }_{\tilde{n}}\triangleq \frac{{\overline{L}}_{\tilde{n}}}{{\overline{\mu }}_{\tilde{n}}}.$$

$$\underset{n\to \infty }{\mathrm{limsup}}\frac{\Vert {\mathbf{x}}_{n+1}-\overline{\mathbf{x}}\Vert }{\Vert {\mathbf{x}}_n-\overline{\mathbf{x}}\Vert }\le \iota .$$

$${\tilde{f}}_{i,r}({\tilde{\mathbf{x}}}_{s,t,r})$$

$$\nabla {\tilde{f}}_{i,r}({\tilde{\mathbf{x}}}_{s,t,r})$$

$${\nabla }^2{\tilde{f}}_{i,r}({\tilde{\mathbf{x}}}_{s,t,r})$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}\left[{\Vert {\mathbf{v}}_{s,t}-\nabla f({\mathbf{x}}_{s,t})\Vert }^2|{\mathcal{F}}_{s,t}\right]$$

$$\le \frac{4\overline{L}}{b}\left(f({\mathbf{x}}_{s,t})-f({\mathbf{x}}^{\ast })+f({\mathbf{x}}^s)-f({\mathbf{x}}^{\ast })\right).$$

$$4L_{\max }\left(f({\mathbf{x}}_{s,t})-f({\mathbf{x}}^{\ast })+f({\mathbf{x}}^s)-f({\mathbf{x}}^{\ast })\right).$$

$$4L_{\max }((f({\mathbf{x}}_{s,t})-f({\mathbf{x}}^{\ast }))$$

$$+\left(1-1/{\kappa }_{\max }\right)(f({\mathbf{x}}^s)-f({\mathbf{x}}^{\ast }))).$$

$$\gamma \triangleq \frac{1}{(M+1){\overline{L}}_{b_{\mathrm{H}}}}\text{mbox}and\Gamma \triangleq \frac{{\kappa }_{b_{\mathrm{H}}}^{M+1}}{{\overline{\mu }}_{b_{\mathrm{H}}}({\kappa }_{b_{\mathrm{H}}}-1)}.$$

$${\mathbf{B}}_r^{(k)}={\mathbf{B}}_r^{(k-1)}-\frac{{\mathbf{B}}_r^{(k-1)}{\mathbf{s}}_k{{\mathbf{s}}_k}^T{\mathbf{B}}_r^{(k-1)}}{{{\mathbf{s}}_k}^T{\mathbf{B}}_r^{(k-1)}{\mathbf{s}}_k}+\frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{s}}_k^T{\mathbf{y}}_k},$$

$$\tilde{\gamma }=\frac{1}{(d+M)L_{\max }}\text{mbox}and\tilde{\Gamma }$$

$${\kappa }_{\mathrm{H}}$$

$${\tilde{\kappa }}_{\mathrm{H}}$$

$$\mathbb{E}$$

$$\rho$$

$${\tilde{\mathbf{x}}}_{s,t+1,r}={\tilde{\mathbf{x}}}_{s,t,r}-\eta {\tilde{\mathbf{v}}}_{s,t,r}.$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}\left[{\tilde{\mathbf{v}}}_{s,t,r}\right|{\mathcal{F}}_{s,t}]={\mathbf{H}}_r^{1/2}\nabla f({\mathbf{x}}_{s,t})=\nabla {\tilde{f}}_r({\tilde{\mathbf{x}}}_{s,t,r})\text{mbox}and$$

$${\mathbb{E}}_{{\mathcal{B}}_{s,t}}\left[\Vert {\tilde{\mathbf{\delta }}}_{s,t,r}{\Vert }^2\right|{\mathcal{F}}_{s,t}]\le {\Vert {\mathbf{H}}^{1/2}\Vert }^2\frac{4\overline{L}}{b}(f({\mathbf{x}}_{s,t})-f({\mathbf{x}}^{\ast })$$

$$+f({\mathbf{x}}^s)-f({\mathbf{x}}^{\ast }))$$

$$\le \frac{4\Gamma \overline{L}}{b}\left({\tilde{f}}_r({\tilde{\mathbf{x}}}_{s,t,r})-{\tilde{f}}_r({\tilde{\mathbf{x}}}_r^{\ast })+{\tilde{f}}_{r^{\prime }}({\tilde{\mathbf{x}}}^{s,r^{\prime }})-{\tilde{f}}_{r^{\prime }}({\tilde{\mathbf{x}}}_{r^{\prime }}^{\ast })\right),$$

$${\Vert {\tilde{\mathbf{x}}}_{s,t+1,r}-{\tilde{\mathbf{x}}}_r^{\ast }\Vert }^2={\Vert {\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }\Vert }^2$$

$$+2\eta \left(\frac{\eta }{2}{\Vert {\tilde{\mathbf{v}}}_{s,t,r}\Vert }^2-⟨{\tilde{\mathbf{v}}}_{s,t,r},{\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }⟩\right).$$

$$\frac{\eta }{2}{\Vert {\tilde{\mathbf{v}}}_{s,t,r}\Vert }^2-⟨{\tilde{\mathbf{v}}}_{s,t,r},{\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }⟩\le -\left({\tilde{f}}_r({\tilde{\mathbf{x}}}_{s,t+1,r})-{\tilde{f}}_r({\tilde{\mathbf{x}}}_r^{\ast })\right)$$

$$-⟨{\tilde{\mathbf{\delta }}}_{s,t,r},{\tilde{\mathbf{x}}}_{s,t+1,r}-{\tilde{\mathbf{x}}}_r^{\ast }⟩-\frac{\gamma \overline{\mu }}{2}{\Vert {\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }\Vert }^2$$

$${\tilde{f}}_r({\tilde{\mathbf{x}}}_{s,t+1,r})+\frac{\eta }{2}{\Vert {\tilde{\mathbf{v}}}_{s,t,r}\Vert }^2-⟨{\tilde{\mathbf{v}}}_{s,t,r},{\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }⟩+⟨{\tilde{\mathbf{\delta }}}_{s,t,r},{\tilde{\mathbf{x}}}_{s,t+1,r}-{\tilde{\mathbf{x}}}_r^{\ast }⟩+\frac{\gamma \overline{\mu }}{2}{\Vert {\tilde{\mathbf{x}}}_{s,t,r}-{\tilde{\mathbf{x}}}_r^{\ast }\Vert }^2$$

$$=$$

Full text

I Introduction

We are interested in the following (unconstrained) convex finite-sum minimization problem

[TABLE]

where $d$ and $n$ denote the ambient dimension of the decision vector and the number of component functions respectively. Problems in the form of (1) play important roles in machine learning and signal processing. One important class of such problems is the empirical risk minimization (ERM) problem, where each $f_{i}$ assumes the form

[TABLE]

In (2), $\ell:\mathbb{R}\times\mathbb{R}\to\mathbb{R}_{+}$ denotes a smooth loss function, $R:\mathbb{R}^{d}\to\mathbb{R}_{+}$ a smooth convex regularizer (e.g., Tikhonov), $\lambda\geq 0$ the regularization weight and $\{(\mathbf{a}_{i},b_{i})\}_{i=1}^{n}\subseteq\mathbb{R}^{d+1}$ the set of feature-response pairs. Depending on the form of $\ell$ and $R$, many important machine learning problems—such as logistic regression, ridge regression and soft-margin support vector machines—are special cases of ERM.

We focus on the case where both $n$ and $d$ are large, and $f$ is ill-conditioned (i.e., the condition number of $f$ is large).111In this work, the condition number of a (strongly) convex function refers to that of its Hessian. In the context of ERM, this means the dataset $\{(\mathbf{a}_{i},b_{i})\}_{i=1}^{n}$ that defines (1) is large and the feature vectors $\mathbf{a}_{i}$ have high ambient dimension. However, the points typically belong to a low-dimensional manifold. Such a setting is particularly relevant in the big-data era, due to unprecedented data acquisition abilities.

When $n$ is large, the computational costs incurred by the batch optimization methods (both first- and second-order) are prohibitive, since in such methods the gradients of all the component functions $\{f_{i}\}_{i=1}^{n}$ need to be computed at each iteration. Therefore, stochastic (randomized) optimization methods have become very popular. At each iteration, only a subset of component functions, rather than all of them, are processed. In this way, for a given time budget, much more progress can be made towards global optima compared to a single-step taken for batch methods. When $d$ is large, Newton- or quasi-Newton-based methods  (both batch and stochastic) incur both high computational and storage complexities. Consequently only first-order and limited-memory quasi-Newton methods (e.g., L-BFGS [2]) are practical in this setting.

I-A Related Works

When both $n$ and $d$ are large, as in our setting, most of research efforts have been devoted to stochastic first-order methods, which include stochastic gradient descent (SGD) [3, 4] and its variance-reduced modifications [5, 6, 7, 8]. However, these methods do not make use of the curvature information. This limits their abilities to find highly accurate solutions for ill-conditioned problems. In order to incorporate the curvature information in the limited-memory setting, recently much progress have been made toward developing stochastic L-BFGS algorithm. A partial list of such works includes [9, 10, 11, 12, 13, 14, 15, 16, 17]. In particular, the first convergent algorithm was proposed by Mokhtari and Ribeiro [12]. Later, the algorithm in [15] makes use of the subsampled Hessian-vector products to form the correction pairs (as opposed to using difference of stochastic gradients) and achieves better results than previous methods. However, the convergence rate is sublinear (in the strongly-convex case), similar to that of SGD. Later, Moritz et al. [16] combines this method with stochastic variance-reduced gradient (SVRG) and proves linear convergence of the resulting algorithm. The algorithm in [17] maintains the structure of this algorithm but incorporates the block BFGS update to collect more curvature information in the optimization process. Although the convergence rate of this new method is similar to that in [16], experimental results demonstrate practical speedups introduced by the block BFGS update. Finally, there also exist a large volume of works on decentralized second-order methods [18, 19, 20, 21, 22, 23, 24] that aim to coordinate multiple distributed agents (with computational and storage abilities) in the optimization task. Since we are not concerned with decentralized optimization algorithms in this paper, we do not discuss these works in details here.

I-B Motivations and Main Contributions

Our work can be motivated from both theory and practice. In terms of theory, although linear convergence (in expectation) has been shown for both algorithms in [16] and [17], the convergence rates (and hence computational complexities) therein can be potentially further improved. (The analysis method in [17] mainly follows that in [16], so we treat the analyses in both works in a unified manner.) In addition, these results may be strengthened in a probabilistic sense, e.g., from convergence in expectation to convergence in probability or almost surely. In terms of practice, in addition to block BFGS update, there may exist several other practical strategies that can potentially further accelerate222In this work, we refer “acceleration” to general strategies that speed up the algorithm, not necessarily the ones based on momentum methods. the algorithm in [16]. Based on these two aspects, our work consists of the following main contributions.

1. We propose a coordinate transformation framework to analyze stochastic L-BFGS-type algorithms in [16] and [17]. Our analysis framework yields a much improved (linear) convergence rate (both in expectation and almost surely) and computational complexity. The essential idea of our method is to unify the analysis of stochastic first-order and second-order methods; as a result, it opens new avenues of designing and analyzing other variants of stochastic second-order algorithms based on their first-order counterparts.

2. We conduct a computational complexity analysis for the stochastic L-BFGS algorithms, which is the first of its kind.

3. We propose several practical acceleration strategies to speed up the convergence of the stochastic L-BFGS algorithm in [16]. We demonstrate the efficacy of our strategies through numerical experiments on logistic and ridge regression problems. We also prove linear convergence for most of these strategies.

II Preliminaries

II-A Notations

We use lowercase, bold lowercase and bold uppercase letters to denote scalars, vectors and matrices respectively. For a matrix $\mathbf{U}\in\mathbb{R}^{m_{1}\times m_{0}}$, we denotes its $(p,q)$-th entry as $u_{pq}$. For a function $f:\mathbb{R}^{m_{1}}\to\mathbb{R}^{m_{2}}$, define the function $f\circ\mathbf{U}$ as the composition $(f\circ\mathbf{U})(\mathbf{z})\triangleq f(\mathbf{U}\mathbf{z})$, for any $\mathbf{z}\in\mathbb{R}^{m_{0}}$. A continuously differentiable function $g:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-smooth ($L\!>\!0$) if and only if $\nabla g$ is $L$-Lipschitz on $\mathbb{R}^{d}$. We use $\mathbb{N}$ to denote the set of natural numbers. For any $n\in\mathbb{N}$, we define $[n]\triangleq\{1,\ldots,n\}$ and $(n]\triangleq\{0,1,\ldots,n\}$. Accordingly, for a sequence of sets $\{\mathcal{A}_{n}\}_{n\geq 0}$, define $\mathcal{A}_{(n]}\triangleq\{\mathcal{A}_{0},\ldots,\mathcal{A}_{n}\}$, for any $n\in\mathbb{N}$. As usual, $\liminf_{n\to\infty}\mathcal{A}_{n}\triangleq\cup_{n\geq 0}\cap_{j\geq n}\mathcal{A}_{j}$ and $\limsup_{n\to\infty}\mathcal{A}_{n}\triangleq\cap_{n\geq 0}\cup_{j\geq n}\mathcal{A}_{j}$. For a set $\mathcal{A}$, denote its complement as $\mathcal{A}^{\mathrm{c}}$. For any sequence $\{x_{i}\}_{i\geq 0}$, we define $\sum_{i=p}^{q}x_{i}\triangleq 0$ if $p>q$. We use $\left\|\cdot\right\|$ to denote both the $\ell_{2}$ norm of a vector and the spectral norm of a matrix. We use $\mathbf{B}$ and $\mathbf{H}$ (with subscripts and superscripts) to denote the approximate Hessian and approximate inverse Hessian in L-BFGS algorithms respectively, following the convention in [25]. $\mathbf{H}$ is also known as the metric matrix [26]. In this work, technical lemmas (whose indices begin with ‘T’) will appear in Appendix E.

II-B Assumptions on Component Functions $f_{i}$

Assumption 1**.**

For each $i\in[n]$, $f_{i}$ is convex and twice differentiable on $\mathbb{R}^{d}$. For ERM problems (2), we assume these two properties are satisfied by the loss function $\ell$ in its first argument on $\mathbb{R}$ and by the regularizer $R$ on $\mathbb{R}^{d}$.

Assumption 2**.**

For each $i\in[n]$, $f_{i}$ is $\mu_{i}$-strongly convex and $L_{i}$-smooth on $\mathbb{R}^{d}$, where $0<\mu_{i}\leq L_{i}$.

Remark 1*.*

Assumptions 1 and 2 are standard in the analysis of both deterministic and stochastic second-order optimization methods. The strong convexity of $f_{i}$ in Assumption 2 ensures positive curvature at any point in $\mathbb{R}^{d}$, which in turn guarantees the well-definedness of the BFGS update. As a common practice in the literature [15, 16], this condition can typically be enforced by adding a strongly convex regularizer (e.g., Tikhonov) to $f_{i}$. Due to the strong convexity, (1) has a unique solution, denoted as $\mathbf{x}^{*}$.

III Algorithm

We will provide a refined analysis of the optimization algorithm (with some modifications) suggested in [16] and so we recapitulate it in Algorithm 1. This algorithm can be regarded as a judicious combination of SVRG and L-BFGS algorithms. We use $s$ and $t$ to denote the outer and inner iteration indices respectively and $r$ to denote the index of metric matrices $\{\mathbf{H}_{r}\}_{r\geq 0}$. We also use $\mathbf{x}_{s,t}$ and $\mathbf{x}^{s}$ to denote an inner iterate and outer iterate respectively.

Each outer iteration $s$ consists of $m$ inner iterations. Before the inner iterations, we first compute a full gradient $\mathbf{g}_{s}\triangleq\nabla f(\mathbf{x}^{s})$. In each inner iteration $(s,t)$, the only modification that we make with respect to (w.r.t.) the original algorithm in [16] is that in computing the stochastic gradient $\mathbf{v}_{s,t}$, the index set $\mathcal{B}_{s,t}\subseteq[n]$ of size $b$ is sampled with replacement nonuniformly [27, 28]. Specifically, the elements in $\mathcal{B}_{s,t}$ are sampled i.i.d. from a discrete distribution $P\triangleq(p_{1},\ldots,p_{n})$, such that for any $i\in[n]$, $p_{i}={L_{i}}/{\sum_{j=1}^{n}L_{j}}$. As will be seen in Lemma 3, compared to uniform sampling, nonuniform sampling leads to a better variance bound on the stochastic gradient $\mathbf{v}_{s,t}$. Using $\mathcal{B}_{s,t}$ and $\nabla f(\mathbf{x}^{s})$, we then compute $\mathbf{v}_{s,t}$ according to (7) in Algorithm 1, where

[TABLE]

This specific way to construct $\mathbf{v}_{s,t}$ reduces the variance of $\mathbf{v}_{s,t}$ to zero as $s\to\infty$ (see Lemma 3), and serves as a crucial step in the SVRG framework.

Then we compute the search direction $\mathbf{H}_{r}\mathbf{v}_{s,t}$. The metric matrix $\mathbf{H}_{r}$ serves as an approximate of the inverse Hessian matrix and therefore contains the local curvature information at the recent iterates. Consequently, $\mathbf{H}_{r}\mathbf{v}_{s,t}$ may act as a better descent direction than $\mathbf{v}_{s,t}$. Since storing $\mathbf{H}_{r}$ may incur high storage cost (indeed, $\Theta(d^{2})$ space) for large $d$, (stochastic) L-BFGS-type methods compute $\mathbf{H}_{r}\mathbf{v}_{s,t}$ each time from a set of recent correction pairs $\mathcal{H}_{r}$ (that only occupies $\Theta(d)$ space) and $\mathbf{v}_{s,t}$. In this way, the limitation on memory can be overcome.

Denote $M\in\mathbb{N}$ as the memory parameter. We next describe the construction of the set of recent correction pairs $\mathcal{H}_{r}\triangleq\{(\mathbf{s}_{j},\mathbf{y}_{j})\}_{j=r-M^{\prime}+1}^{r}$, where $M^{\prime}\!\triangleq\!\min\{r,M\}$. Together with $\mathbf{v}_{s,t}$, this set will be used to compute the matrix-vector product $\mathbf{H}_{r}\mathbf{v}_{s,t}$. Before doing so, in line 13, we first compute the averaged past iterates $\{\overline{\bf x}_{r}\}_{r\geq 0}$ from $\{\mathbf{x}_{s,t}\}_{s\geq 0,t\in(m]}$ for every $\Upsilon$ inner iterates, where $\Upsilon\!\in\![m]$. Based on $\overline{\bf x}_{r-1}$ and $\overline{\bf x}_{r}$, we compute the most recent correction pair $(\mathbf{s}_{r},\mathbf{y}_{r})$ in line 16. Following [15], in computing $\mathbf{y}_{r}$, we first sample an index set $\mathcal{T}_{r}\subseteq[n]$ of size $b_{\mathrm{H}}$ uniformly without replacement, and then let $\mathbf{y}_{r}=\nabla^{2}f_{\mathcal{T}_{r}}(\overline{\bf x}_{r})\mathbf{s}_{r}$, where

[TABLE]

denotes the sub-sampled Hessian at $\overline{\bf x}_{r}$. Finally, we update $\mathcal{H}_{r-1}$ to $\mathcal{H}_{r}$ by inserting $(\mathbf{s}_{r},\mathbf{y}_{r})$ into $\mathcal{H}_{r-1}$ and deleting $(\mathbf{s}_{r-M^{\prime}},\mathbf{y}_{r-M^{\prime}})$ from it.

Based on $\mathcal{H}_{r}$, a direct approach to compute $\mathbf{H}_{r}\mathbf{v}_{s,t}$ would be computing $\mathbf{H}_{r}$ first and then forming the product with $\mathbf{v}_{s,t}$. Computing $\mathbf{H}_{r}$ involves applying $M^{\prime}$ BFGS updates to

[TABLE]

using $\{(\mathbf{s}_{j},\mathbf{y}_{j})\}_{j=r-M^{\prime}+1}^{r}$. For each $k\in\{r-M^{\prime}+1,\ldots,r\}$, the update is

[TABLE]

Finally we set $\mathbf{H}_{r}=\mathbf{H}_{r}^{(r)}$. Instead of using this direct approach, we adopt the two-loop recursion algorithm [25, Algorithm 7.4] to compute $\mathbf{H}_{r}\mathbf{v}_{s,t}$ as a whole. This method serves the same purpose as the direct one, but, as we shall see, with much reduced computation.

At the end of each outer iteration $s$, the starting point of next outer iteration, $\mathbf{x}^{s+1}$ is either uniformly sampled (option I) or averaged (option II) from all the past inner iterates $\{\mathbf{x}_{s,t}\}_{t\in[m]}$. As shown in Theorem 1, these two options can be analyzed in a unified manner.

Remark 2*.*

Under many scenarios (e.g., the ridge and logistic regression problems in Section VII), the smoothness parameters $\{L_{i}\}_{i=1}^{n}$ can be accurately estimated. (For ERM problems, these parameters are typically data-dependent.) If in some cases, accurate estimates of these parameters are not available, we can simply employ uniform sampling, which is a special case of our weighted sampling technique.

IV Convergence Analysis

IV-A Definitions

Let $\{\underline{i}_{j}\}_{j\in[n]}$ and $\{\overline{i}_{j}\}_{j\in[n]}$ be permutations of $[n]$ such that $\mu_{\min}\triangleq\mu_{\underline{i}_{1}}\leq\cdots\leq\mu_{\underline{i}_{n}}$ and $L_{\overline{i}_{1}}\leq\cdots\leq L_{\overline{i}_{n}}\triangleq L_{\max}$. Given any $\widetilde{n}\in[n]$, define

[TABLE]

Accordingly, define

[TABLE]

In particular, define $\overline{\mu}\triangleq\overline{\mu}_{n}$, $\overline{L}\triangleq\overline{L}_{n}$ and $\kappa\triangleq\overline{L}/\overline{\mu}$. Denote the probability space on which the sequence of (random) iterates $\{\mathbf{x}_{s,t}\}_{s\geq 0,t\in(m]}$ in Algorithm 1 is defined as $(\Omega,\Sigma,\mathbb{P})$, where $\Sigma$ is the Borel $\sigma$-algebra of $\Omega$. We also define a filtration $\{\mathcal{F}_{s,t}\}_{s\geq 0,t\in(m-1]}$ such that $\mathcal{F}_{s,t}$ contains all the information up to the time $(s,t)$. Formally, $\mathcal{F}_{s,t}\triangleq\sigma(\{\tau_{j}\}_{j=0}^{s-1}\cup\{\mathcal{B}_{i,j}\}_{i\in(s-1],j\in(m-1]}\cup\{\mathcal{B}_{s,j}\}_{j=0}^{t-1}\cup\{\mathcal{T}_{j}\}_{j=0}^{\lfloor{(sm+t)/L}\rfloor})$, where $\sigma\left(\{x_{j}\}_{j=1}^{n}\right)$ denotes the $\sigma$-algebra generated by random variables $\{x_{j}\}_{j=1}^{n}$. Define $\mathcal{F}_{s}\triangleq\mathcal{F}_{s,0}$.

To introduce our coordinate transformation framework, we define some transforms of variables appearing in Algorithm 1. Specifically, for any $s,t,r\geq 0$, define $\widetilde{\mathbf{x}}_{s,t,r}\triangleq\mathbf{H}_{r}^{-1/2}\mathbf{x}_{s,t}$, $\widetilde{\mathbf{x}}^{s,r}\triangleq\mathbf{H}_{r}^{-1/2}\mathbf{x}^{s}$, $\widetilde{\mathbf{x}}_{r}^{*}\triangleq\mathbf{H}_{r}^{-1/2}\mathbf{x}^{*}$ and $\widetilde{\mathbf{v}}_{s,t,r}\triangleq\mathbf{H}_{r}^{1/2}\mathbf{v}_{s,t}$.333As will be shown in Lemma 4, for any $r\geq 0$, $\mathbf{H}_{r}\succeq\gamma\mathbf{I}$ for some $\gamma>0$. Therefore, $\mathbf{H}_{r}^{1/2}$ (and $\mathbf{H}_{r}^{-1/2}$) are well-defined. We also define transformed functions $\widetilde{f}_{i,r}\triangleq f_{i}\circ\mathbf{H}_{r}^{1/2}$ and $\widetilde{f}_{r}\triangleq\frac{1}{n}\sum_{i=1}^{n}\widetilde{f}_{i,r}$, for any $i\in[n]$ and $r\geq 0$.

To state our convergence results, we define the notions of linear convergence and R-linear convergence.

Definition 1**.**

A sequence $\{\mathbf{x}_{n}\}_{n\geq 0}\!\subseteq\!\mathbb{R}^{d}$ is said to converge to $\overline{\bf x}\!\in\!\mathbb{R}^{d}$ linearly (or more precisely, Q-linearly) with rate $\iota\in(0,1)$ if

[TABLE]

We say $\mathbf{x}_{n}\to\overline{\bf x}$ R-linearly with rate $\iota^{\prime}\in(0,1)$ if there exists a nonnegative sequence $\{\varepsilon_{n}\}_{n\geq 0}$ such that $\|\mathbf{x}_{n}-\overline{\bf x}\|\leq\varepsilon_{n}$ for sufficiently large $n$ and $\varepsilon_{n}\to 0$ linearly with rate $\iota^{\prime}$.

IV-B Preliminary Lemmas

From the definitions of transformed variables and functions in Section IV-A, we immediately have the following lemmas.

Lemma 1**.**

For any $s,t,r\geq 0$ and $i\in[n]$, we have

[TABLE]

Lemma 2**.**

If there exist $0<\gamma^{\prime}\leq\Gamma^{\prime}$ such that $\gamma^{\prime}\mathbf{I}\preceq\mathbf{H}_{r}\preceq\Gamma^{\prime}\mathbf{I}$ for all $r\geq 0$, then for any $i\in[n]$ and $r\geq 0$, $\widetilde{f}_{i,r}$ is twice differentiable, $(\mu_{i}\gamma^{\prime})$-strongly convex and $(L_{i}\Gamma^{\prime})$-smooth on $\mathbb{R}^{d}$. Consequently, $\widetilde{f}_{r}$ is twice differentiable, $(\gamma^{\prime}\overline{\mu})$-strongly convex and $(\Gamma^{\prime}\overline{L})$-smooth on $\mathbb{R}^{d}$.

Next we derive two other lemmas that will not only be used in the analysis later, but have the potential to be applied to more general problem settings. Specifically, Lemma 3 can be applied to any stochastic optimization algorithms based on SVRG and Lemma 4 can be applied to any finite-sum minimization algorithms based on L-BFGS methods (not necessarily stochastic in nature). The proofs of Lemmas 3 and 4 are deferred to Appendices A and B respectively.

Lemma 3** (Variance bound of $\mathbf{v}_{s,t}$).**

In Algorithm 1, we have $\mathbb{E}_{\mathcal{B}_{s,t}}\left[\mathbf{v}_{s,t}|\mathcal{F}_{s,t}\right]=\nabla f(\mathbf{x}_{s,t})$ and

[TABLE]

Remark 3*.*

In previous works [16, 17], a uniform mini-batch sampling of $\mathcal{B}_{s,t}$ was employed, and different variance bounds of $\mathbf{v}_{s,t}$ were derived. In [16, Lemma 6], the bound was

[TABLE]

In [17, Lemma 2], this bound was slightly improved to

[TABLE]

However, both of these bounds fail to capture the dependence on the mini-batch size $b$. In contrast, in this work we consider a nonuniform mini-batch sampling (with replacement). Due to division by $b$ and $\overline{L}\leq L_{\max}$ (indeed in many cases, $\overline{L}\ll L_{\max}$), our bound in (14) is superior to (15) and (16). As will be seen in Theorem 1, our better bound (14) leads to a faster (linear) convergence rate of Algorithm 1.

Lemma 4** (Uniform Spectral Bound of $\{\mathbf{H}_{r}\}_{r\geq 0}$).**

The spectra of $\{\mathbf{H}_{r}\}_{r\geq 0}$ are uniformly bounded, i.e., for each $r\geq 0$, $\gamma\mathbf{I}\preceq\mathbf{H}_{r}\preceq\Gamma\mathbf{I}$, where444We assume $\kappa_{b_{\mathrm{H}}}>1$ for any ${b_{\mathrm{H}}}\in[n]$ since we focus on the setting where $f$ is ill-conditioned, i.e., $\kappa_{b_{\mathrm{H}}}\geq\kappa\gg 1$. If $\kappa_{b_{\mathrm{H}}}=1$ for some ${b_{\mathrm{H}}}\in[n]$, then $\Gamma=(\kappa_{b_{\mathrm{H}}}^{M}+M)/\overline{\mu}_{b_{\mathrm{H}}}$ and $\gamma$ remains the same. The proof for this case can be straightforwardly adapted from that in Section B.

[TABLE]

Remark 4*.*

In [13, 15, 16], the authors make use of a classical technique in [2] to derive a different uniform spectral bound of $\{\mathbf{H}_{r}\}_{r\geq 0}$. Their technique involves applying $\operatorname{tr}(\cdot)$ and $\det(\cdot)$ recursively to the BFGS update rule

[TABLE]

where $\mathbf{B}_{r}^{(k)}\triangleq(\mathbf{H}_{r}^{(k)})^{-1}$ denotes the approximate Hessian matrix at step $k$ in the reconstruction of $\mathbf{B}_{r}\triangleq(\mathbf{H}_{r})^{-1}$. The lower and upper bounds derived by this technique are

[TABLE]

respectively. As will be seen in Proposition 1, the overall computational complexity of Algorithm 1 heavily depends on the estimated (uniform) condition number of $\{\mathbf{H}_{r}\}_{r\geq 0}$. Therefore, it is instructive to compute this quantity for both $(\gamma,\Gamma)$ and $(\widetilde{\gamma},\widetilde{\Gamma})$ as

[TABLE]

where the approximation in (19) follows from $\kappa_{b_{\mathrm{H}}}\gg 1$ (see Footnote 4). By comparing (19) and (20), we notice our estimate for the condition number of $\{\mathbf{H}_{r}\}_{r\geq 0}$, namely $\kappa_{\mathrm{H}}$, is smaller than those in [15] and [16], namely $\widetilde{\kappa}_{\mathrm{H}}$, in several aspects. First, $\kappa_{\mathrm{H}}$ does not grow (exponentially) with the data dimension $d$. Second, $\kappa_{\mathrm{H}}$ depends on $\kappa_{b_{\mathrm{H}}}$, which is usually much smaller than $\kappa_{\max}$. Third, even if we set $d=1$ in (20), the factor $M+1$ in (19) is much smaller than the factor $(M+1)^{M+1}$ in (20). As a result, our improved estimate of the condition number of $\{\mathbf{H}_{r}\}_{r\geq 0}$ will lead to a much better computational complexity estimate (see Proposition 1).

IV-C Main Results

Our main convergence results consist of Theorem 1 and Corollary 1, which provide linear convergence guarantees of $f(\mathbf{x}^{s})$ to $f(\mathbf{x}^{*})$ in expectation and almost surely, respectively.

Theorem 1**.**

In Algorithm 1, choose $\eta\!<\!\min\{b/12,1\}/(\Gamma\overline{L})$ and $m$ sufficiently large. With either option I or II, we have

[TABLE]

Proof.

Fix an outer iteration $s$ and consider an inner iteration $t$. Define $r\triangleq\lfloor{(sm+t)/L}\rfloor$. For brevity, we omit the dependence of $r$ on $s$ and $t$. The iteration in line 10 of Algorithm 1 becomes

[TABLE]

Define ​$\widetilde{\bm{\delta}}_{s,t,r}\!\triangleq\!\widetilde{\mathbf{v}}_{s,t,r}\!-\!\nabla\widetilde{f}_{r}(\widetilde{\mathbf{x}}_{s,t,r})$. From Lemmas 1 and 3,

[TABLE]

where $r^{\prime}\triangleq\lfloor{sm/L}\rfloor$. Using (23), we can express the distance between $\widetilde{\mathbf{x}}_{s,t+1,r}$ and $\widetilde{\mathbf{x}}_{r}^{*}$ as

[TABLE]

We can show

[TABLE]

from steps (28) to (32) on the next page.

In (30), we use the condition $\eta\!\leq\!1/(\Gamma\overline{L})$. In (32), we use the $(\Gamma\overline{L})$-smoothness of $\widetilde{f}_{r}$ and the $(\gamma\overline{\mu})$-strong convexity of $\widetilde{f}_{r}$ in Lemma 2 respectively.

Now, substituting (27) into (26), we have

[TABLE]

Taking expectation w.r.t. $\mathcal{B}_{s,t}$ and using (24) and (25), we have

[TABLE]

By bounding the factor $1-\eta\gamma\overline{\mu}$ by $1$, we can telescope (34) over $t=0,\ldots,m-1$ and obtain

[TABLE]

If we use option II to choose $\mathbf{x}^{s+1}$ (in line 21), we have $\widetilde{\mathbf{x}}^{s+1,r^{\prime\prime}}=1/m\sum_{i=1}^{m}\widetilde{\mathbf{x}}_{s,t,r^{\prime\prime}}$, where $r^{\prime\prime}\triangleq\lfloor{(s+1)m/L}\rfloor$. Using (11) and Jensen’s inequality, we have

[TABLE]

Alternatively, if we use option I to determine $\mathbf{x}^{s+1}$ (in line 20), we still have (36) (with inequality replaced by equality). If we further use Lemma T-1 to upper bound the term $\|\widetilde{\mathbf{x}}^{s,r^{\prime}}-\widetilde{\mathbf{x}}_{r^{\prime}}^{*}\|^{2}$ in (35), we have

[TABLE]

Using (11) again and rearranging, we have

[TABLE]

We take expectation on both sides to complete the proof. ∎

Remark 5*.*

We compare our linear convergence rate $\rho$ in (22) with those in [16] and [17]. Since the convergence rates in these two works are almost the same, we use the rate in [16] for comparison. The linear rate $\widetilde{\rho}$ in [16] equals

[TABLE]

For simplicity, if we let $b=1$, $\mu_{\min}=\overline{\mu}$, $L_{\max}=\overline{L}$, $\widetilde{\gamma}=\gamma$, $\widetilde{\Gamma}=\Gamma$ and ignore other constant factors,555However, bear in mind that by doing all these substitutions, $\widetilde{\rho}$ has already been much improved. we notice that there is an additional multiplicative factor $\kappa\kappa_{\mathrm{H}}$ associated with $\eta\Gamma\overline{L}$ in $\widetilde{\rho}$. As a result $\widetilde{\rho}>\rho$. A more direct way to observe the detrimental effects of this additional $\kappa\kappa_{\mathrm{H}}$ is to compare the computational complexities resulting from $\rho$ and $\widetilde{\rho}$. See Remark 9 for details. The reason that we manage to avoid this factor in our rate $\rho$ is precisely because we adopt the coordinate transformation framework in our analysis (see proof above). Specifically, by absorbing the sequence of metric matrices $\{\mathbf{H}_{r}\}_{r\geq 0}$ into decision vectors and functions, we are able to proceed through bounding the (expected squared Euclidean) distance between $\widetilde{\mathbf{x}}_{s,t+1,r}$ and $\widetilde{\mathbf{x}}_{r}^{*}$, instead of directly bounding $f(\mathbf{x}_{s,t+1})$ via the smoothness property of $f$ (cf. proof of Theorem 7 in [16]). Thus in our analysis, we avoid the additional appearance of $\overline{L}$ and $\Gamma$ (which leads to the additional factor $\kappa\kappa_{\mathrm{H}}$).

Corollary 1**.**

In Algorithm 1, $\{f(\mathbf{x}^{s})\}_{s\geq 0}$ converges to $f(\mathbf{x}^{*})$ R-linearly almost surely with rate $\rho$.

Proof.

Our proof is inspired by [29, Corollary 2] and mainly leverages the Borel-Cantelli lemma [30]. For any $\epsilon^{\prime}>0$ and $0\!<\!\delta\!<\!1-\rho$, consider the sequence of events $\{\mathcal{E}_{s}\}_{s\in\mathbb{N}}$ such that

[TABLE]

Therefore,

[TABLE]

where in (a) we use Markov’s inequality, in (b) we use (21) in Theorem 1, in (c) we use

[TABLE]

and in (d) we use the fact that $s\geq 1$. Thus by the Borel-Cantelli lemma, ${\mathbb{P}}(\limsup_{s\to\infty}\!\mathcal{E}_{s})\!=\!0$, or equivalently,

[TABLE]

The definition of $\mathcal{E}_{s}$ in (39) implies that

[TABLE]

Since $\epsilon^{\prime}>0$ is arbitrary, we have

[TABLE]

or equivalently, $f(\mathbf{x}^{s})-f(\mathbf{x}^{*})=o((\rho+1/\sqrt{s})^{s})$ almost surely, for any $s\geq 0$. For convenience, define a sequence $\{\varpi_{s}\}_{s\geq 0}$ such that $\varpi_{s}\triangleq(\rho+1/\sqrt{s})^{s}$, for any $s\geq 0$. By applying (41) to $\varpi_{s}$, we have $\varpi_{s}=\rho^{s-\Theta(\sqrt{s})}$ (the implied constant in the $\Theta(\cdot)$ notation is positive). Since $\rho^{-\Theta(\sqrt{s+1})}/\rho^{-\Theta(\sqrt{s})}\to 1$ as $s\to\infty$, $\lim_{s\to\infty}\varpi_{s+1}/\varpi_{s}=\rho$ as desired. ∎

Remark 6*.*

Note that the analysis techniques in Corollary 1 can be applied to any stochastic algorithm with linear convergence in expectation, therefore are of independent interest. Although a similar result was proved in [29, Corollary 2], it is weaker than Corollary 1. Specifically, the almost sure linear convergence rate therein is strictly worse than the corresponding linear rate in expectation. In contrast, by leveraging refined analysis techniques, we show that such a degradation can indeed be avoided.

Remark 7*.*

By the $\overline{\mu}$-strong convexity of $f$ (see Assumption 2), we have $\left\|\mathbf{x}-\mathbf{x}^{*}\right\|^{2}\!\leq\!(2/\overline{\mu})(f(\mathbf{x})-f(\mathbf{x}^{*}))$, for any $\mathbf{x}\!\in\!\mathbb{R}^{d}$. Therefore, the linear convergence of $\{f(\mathbf{x}^{s})\}_{s\geq 0}$ to $f(\mathbf{x}^{*})$ in expectation (in Theorem 1) also implies the R-linear convergences of $\{\mathbf{x}^{s}\}_{s\geq 0}$ to $\mathbf{x}^{*}$ in expectation. Similarly, we can also derive the almost sure R-linear convergence of $\{\mathbf{x}^{s}\}_{s\geq 0}$ to $\mathbf{x}^{*}$ from Corollary 1.

V Complexity Analysis

In this section we provide a framework for analyzing the computational complexity of the stochastic L-BFGS algorithm in Algorithm 1. Our framework can be easily generalized to other stochastic second-order algorithms, e.g., SQN algorithm in [15]. To begin with, we make two additional assumptions.

Assumption 3**.**

For any $\mathbf{x}\in\mathbb{R}^{d}$ and $i\in[n]$, the gradient $\nabla f_{i}(\mathbf{x})$ can be computed in $O(d)$ operations.666An operation refers to evaluation of an elementary function, such as addition, multiplication and logarithm.

Assumption 4**.**

For any $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$ and $i\in[n]$, the Hessian-vector product $\nabla^{2}f_{i}(\mathbf{x})\mathbf{y}$ can be computed in $O(d)$ operations.

Remark 8*.*

These two assumptions are naturally satisfied for ERM problems (2) with Tikhonov regularization. For these problems, $R(\mathbf{x})=\frac{1}{2}\left\|\mathbf{x}\right\|^{2}$ and

[TABLE]

where $\ell^{\prime}(\cdot,\cdot)$ and $\ell^{\prime\prime}(\cdot,\cdot)$ are first and second derivatives of $\ell(\cdot,\cdot)$ w.r.t. the first argument. We easily see that the right-hand sides of both (44) and (45) can be computed in $O(d)$ operations.

From Algorithm 1, we observe that its total computational cost $C$ can be split into three parts. The first part $C_{1}$ involves computing the variance-reduced gradient $\mathbf{v}_{s,t}$ in (7), the second part $C_{2}$ involves computing $\mathbf{H}_{r}\mathbf{v}_{s,t}$ (via two-loop recursion) in line 9, and the third part $C_{3}$ involves computing the correction pair $(\mathbf{s}_{r},\mathbf{y}_{r})$ in line 16.

Proposition 1**.**

Let Assumptions 1 to 4 hold. In Algorithm 1,

[TABLE]

Thus the total computational cost $C\triangleq\sum_{i=1}^{3}C_{i}$ equals

[TABLE]

Proof.

We leverage techniques that have become standard in the SVRG literature (e.g., [27]). In (22), if we choose $\eta=\theta b/(\Gamma\overline{L})$ for some $0<\theta<1/12$, $m=\theta^{\prime}\kappa\kappa_{\mathrm{H}}/b$ for some large enough positive constant $\theta^{\prime}$, and use $1+1/m\leq 2$, then

[TABLE]

As a result, the required number of outer iterations to achieve $\epsilon$-suboptimality is $O(\log(1/\epsilon))$. Thus (46) follows from Assumption 3 and that $2mb$ gradients (of component functions) are computed in each inner iteration. If we further choose $M=\Theta(b)$, then (47) follows from the fact that two-loop recursion can be done in $O(Md)$ time [25, Chapter 7]. Lastly, if we choose $b_{\mathrm{H}}=\Theta(\Upsilon)$, then we obtain (48) using Assumption 4. ∎

Remark 9*.*

Following a similar argument, we can deduce the total complexity estimate $\widetilde{C}$ based on the linear rate $\widetilde{\rho}$ (see Remark 5) derived in [16] as

[TABLE]

Compared with $\widetilde{C}$, we observe that our complexity estimate $C$ in (49) is much better, in several aspects. First, the dependence of $C$ on the condition number $\kappa\kappa_{\mathrm{H}}$ is linear, rather than quadratic. The quadratic dependence of $\kappa_{\max}\widetilde{\kappa}_{\mathrm{H}}$ in $\widetilde{C}$ is precisely caused by the additional $\kappa_{\max}\widetilde{\kappa}_{\mathrm{H}}$ in $\widetilde{\rho}$ (see Remark 5). Second, $C$ is independent of the mini-batch size $b$. The appearance of $b$ in $\widetilde{C}$ is a result of the loose bound on variance of $\mathbf{v}_{s,t}$ (cf. (15) and (16)). Third, the condition number $\kappa\kappa_{\mathrm{H}}$ in $C$ is much more smaller than $\kappa_{\max}\widetilde{\kappa}_{\mathrm{H}}$ in $\widetilde{C}$ for ill-conditioned problems. This is a result of the non-uniform sampling of $\mathcal{B}_{s,t}$ and our improved bound on the spectra of $\{\mathbf{H}_{r}\}_{r\geq 0}$.

Remark 10*.*

As our coordinate transformation framework unifies the design and analysis of stochastic first- and second-order algorithms, we believe that momentum-based acceleration techniques for stochastic first-order methods [31, 32] can be applied to Algorithm 1 as well. (Details are left to future work.) In this case, the dependence on $\kappa\kappa_{\mathrm{H}}$ in $C$ may be further improved to $\sqrt{\kappa\kappa_{\mathrm{H}}}$ [32].

VI Acceleration Strategies

In this section, we propose three practical acceleration strategies. We follow the notations in Section III and Algorithm 1. As will be shown in Section VII-B, all of these strategies result in faster convergence in practice. For the first and second strategies, we also provide their theoretical analyses in Propositions 2 and 3, respectively. See Appendices C and D for the proofs of these two propositions.

VI-A Geometric Sampling/Averaging Scheme

Instead of choosing $\mathbf{x}^{s+1}$ according to option I or II in Algorithm 1, inspired by [33], we can introduce a “forgetting” effect by considering two other schemes:

option III: Sample $\tau_{s}$ randomly from $[m]$ from the distribution $Q\triangleq(\beta^{m-1}/c,\beta^{m-2}/c,\ldots,1/c)$ and set $\mathbf{x}^{s+1}:=\mathbf{x}_{s,\tau_{s}}$,

option IV: $\mathbf{x}^{s+1}:=\frac{1}{c}\sum_{t=1}^{m}{\beta^{m-t}}\mathbf{x}_{s,t}$,

where $0<\beta\leq 1-\eta\gamma\overline{\mu}<1$ and the normalization constant $c\triangleq\sum_{t=1}^{m}{\beta^{m-t}}$. Since $\beta\in(0,1)$, we observe that in both options III and IV, more recent iterates (i.e., iterates $\mathbf{x}_{s,t}$ with larger indices $t$) will have larger contributions to $\mathbf{x}^{s+1}$. Theoretically, these two schemes can be analyzed in a unified manner, as shown in the following proposition.

Proposition 2**.**

In Algorithm 1, choose $\eta\!<\!\min\{b/12,\!1\}/(\Gamma\overline{L})$ and $m$ sufficiently large. With either option III or IV, we have

[TABLE]

[TABLE]

and $c^{\prime}\triangleq c/(1-\eta\gamma\overline{\mu})^{m}$.

Remark 11*.*

In the literature [16, 17], usually option I or II (in Algorithm 1) is analyzed to prove that the stochastic L-BFGS algorithms therein converge linearly. However, for faster convergence in practice, $\mathbf{x}^{s+1}$ is chosen to be the last inner iterate $\mathbf{x}_{s,m}$. However, the latter strategy is not amenable to linear convergence analysis. This gap between theory and practice is filled in by our geometric sampling or averaging scheme, i.e., option III or IV. Specifically, as shown in Figure 2, our scheme not only yields linear convergence in theory, but also performs as well as the “last inner iterate” scheme in practice.

VI-B Subsampled Gradient Stabilization

In Algorithm 1, at the beginning of each outer iteration indexed by $s$, we compute a full gradient $\mathbf{g}_{s}$ to stabilize the subsequent (inner) iterations. Inspired by [7], we propose a strategy that only computes a subsampled gradient $\widetilde{\mathbf{g}}_{s}$ at the start of each outer iteration $s$. Specifically, we uniformly sample a subset $\widetilde{\mathcal{B}}_{s}$ of $[n]$ with size $\widetilde{b}_{s}$ without replacement and then form $\widetilde{\mathbf{g}}_{s}\triangleq({1}/{\widetilde{b}_{s}})\sum_{i\in\widetilde{\mathcal{B}}_{s}}\nabla f_{i}(\mathbf{x}^{s})$. The size of $\widetilde{\mathcal{B}}_{s}$, namely $\widetilde{b}_{s}$, increases with the index $s$ until it reaches $n$. By judiciously choosing $\widetilde{b}_{s}$, we can show that the resulting algorithm still enjoys linear convergence with rate $\overline{\rho}$, when integrated with the geometric sampling/averaging scheme in Section VI-A. Before we formally state this result in Proposition 3, we first make an assumption in addition to Assumptions 1 and 2.

Assumption 5**.**

The inner iterates $\{\mathbf{x}_{s,t}\}_{s\geq 0,t\in[m]}$ generated by the modified algorithm in Section VI-B are bounded almost surely, i.e., there exists $B<\infty$ such that for any $s\geq 0$ and $t\in[m]$, $\left\|\mathbf{x}_{s,t}-\mathbf{x}^{*}\right\|\leq B$.

Proposition 3**.**

Let Assumptions 1, 2 and 5 hold. For any $\xi>0$ and $S\in\mathbb{N}$, and for any $s\in(S]$, if we choose $\widetilde{b}_{s}\geq{\underline{\widetilde{b}}}_{s}\triangleq nS^{2}\alpha_{s}/(S^{2}\alpha_{s}+(n-1)\xi^{2}\overline{\rho}^{2s})$, where $\alpha_{s}\triangleq{1}/{n}\sum_{i=1}^{n}\left\|\nabla f_{i}(\mathbf{x}^{s})\right\|^{2}$, we have

[TABLE]

Remark 12*.*

Several remarks are in order. First, we remark that assumptions involving almost sure boundedness of iterates (e.g., Assumption 5) are commonplace in the stochastic optimization literature [34, 35, 7] and are always observed to hold in experiments. Second, under this assumption, we can show that $\{\alpha_{s}\}_{s\geq 0}$ are bounded almost surely using the Lipschitz continuity of $\nabla f_{i}$, for any $i\in[n]$ in Assumption 2. Consequently, there exists $B^{\prime}<\infty$ such that $\alpha_{s}\leq B^{\prime}$ for any $s\geq 0$ and hence

[TABLE]

As a sanity check, we observe that (55) increases to $n$ as $s\to\infty$. By further upper bounding (55) by (56), we obtain a practical rule to select $\widetilde{b}_{s}$. Namely, it suffices to choose $\widetilde{b}_{s}=\min\{\zeta\upsilon^{s},n\}$, for some constants $\zeta>0$ and $\upsilon>1$. As shown in Section VII, this rule works well in practice. Third, for any $\epsilon>0$, we can choose $S\in\mathbb{N}$ such that our algorithm achieves $\epsilon$-suboptimality, i.e., $\mathbb{E}\left[f(\mathbf{x}^{S})-f(\mathbf{x}^{*})\right]<\epsilon$.

VI-C Low-dimensional Approximate Hessians

In additional to high dimensionality and large size, sparsity is also a typical attribute for modern data, i.e., many feature vectors only have a few nonzero entries. For ERM problems (2) (with Tikhonov regularization), this implies that the Hessian of $f$ in (1),

[TABLE]

tends to be sparse. Based on this observation, we propose a strategy that aims to approximate $\nabla^{2}f(\mathbf{x})$ by several smaller Hessian matrices and update them efficiently. For sparse data, collecting curvature information via smaller dense Hessians can be more effective than directly manipulating the high-dimensional sparse Hessian [25]. As a result, the algorithm converges faster in practice (see Figure 4).

Before describing our strategy, we first introduce some notations. We partition $[n]$ into $K$ groups, and denote the set of partitions as $\mathcal{P}\triangleq\{\mathcal{P}_{1},\ldots,\mathcal{P}_{K}\}$. For any $i\in[K]$, we define $\mathcal{S}_{i}\triangleq\cup_{j\in\mathcal{P}_{i}}\operatorname{supp}(\mathbf{a}_{j})$, where $\operatorname{supp}(\mathbf{a}_{j})$ denotes the support of the vector $\mathbf{a}_{j}$. We define $d_{i}\triangleq\left\lvert\mathcal{S}_{i}\right\rvert$ and denote the elements in $\mathcal{S}_{i}$ as $\{s_{i,1},\ldots,s_{i,d_{i}}\}$. We also define $F_{i}=\sum_{j\in\mathcal{P}_{i}}f_{j}$ so that $f=\frac{1}{n}\sum_{i=1}^{K}F_{i}$. We define a projection matrix $\mathbf{U}_{i}\in\mathbb{R}^{d_{i}\times d}$ such that for any $p\in[d_{i}]$ and $q\in[d]$, $u_{pq}=1$ if $q=s_{i,p}$ and [math] otherwise. Accordingly, for any $l\in\mathcal{P}_{i}$, define a function $\phi_{i,l}:\mathbb{R}^{d_{i}}\to\mathbb{R}$ such that $f_{l}\triangleq\phi_{i,l}\circ\mathbf{U}_{i}$. Note that $\phi_{i,l}$ is uniquely defined by the definition $\mathbf{U}_{i}$. Also define $\phi_{i}\triangleq\sum_{l\in\mathcal{P}_{i}}\phi_{i,l}$. Therefore,

[TABLE]

We now describe our strategy. In Algorithm 1, for any $i\in[K]$ and any $j\in\{r-M^{\prime}+1,\ldots,r\}$, define correction pairs $\mathbf{s}_{j,i}\!\!\triangleq\!\!\mathbf{U}_{i}(\overline{\bf x}_{j}-\overline{\bf x}_{j-1})$ and $\mathbf{y}_{j,i}\!\!\triangleq\!\!\sum_{l\in\mathcal{T}_{r,i}}\!\!\nabla^{2}\phi_{i,l}(\mathbf{U}_{i}\overline{\bf x}_{j})\mathbf{s}_{j,i}$, where $\mathcal{T}_{r,i}$ with size $b_{\mathrm{H}}/K$ is uniformly sampled from $\mathcal{P}_{i}$. Accordingly, define $\mathbf{S}_{r,i}\triangleq[\mathbf{s}_{r-M^{\prime}+1,i},\ldots,\mathbf{s}_{r,i}]$, $\mathbf{Y}_{r,i}\triangleq[\mathbf{y}_{r-M^{\prime}+1,i},\ldots,\mathbf{y}_{r,i}]$. Instead of storing $\mathcal{H}_{r}$, we only store matrices $\{(\mathbf{S}_{r,i},\mathbf{Y}_{r,i})\}_{i=1}^{K}$. To reconstruct approximation $\mathbf{B}_{r,i}$ for each $\nabla^{2}\phi_{i}$ at $\mathbf{U}_{i}\overline{\bf x}_{r}$, as usual, we apply $M^{\prime}$ BFGS updates (18) (using the correction pairs stored in $\mathbf{S}_{r,i}$ and $\mathbf{Y}_{r,i}$) to $\mathbf{B}_{r,i}^{(0)}\triangleq\delta_{r,i}\mathbf{I}$, where $\delta_{r,i}\!\triangleq\!\|\mathbf{y}_{r,i}\|^{2}/{\mathbf{s}_{r,i}^{T}\mathbf{y}_{r,i}}$. This procedure can be implemented efficiently via the method of compact representation [25], i.e.,

[TABLE]

where $\mathbf{W}_{r,i}\!\triangleq\![\delta_{r}^{(i)}{\mathbf{S}_{r}^{(i)}},\mathbf{Y}_{r}^{(i)}]$ and $\mathbf{L}_{r,i}$ and $\mathbf{D}_{r,i}$ are the lower triangular matrix (excluding diagonal) and diagonal matrix of $\mathbf{S}_{r,i}^{T}\mathbf{Y}_{r,i}$ respectively. Analogous to (58), we define the approximation of $\nabla^{2}f$ at $\mathbf{U}_{i}\overline{\bf x}_{r}$, denoted as $\mathbf{B}_{r}$, as

[TABLE]

We remark that the strong convexity of each function $f_{i}$ (see Assumption 2), together with the full-row-rank property of $\mathbf{U}_{i}$, ensures $\nabla^{2}\phi_{i}\!\succ\!0$ on $\mathbb{R}^{d_{i}}$. This implies positive curvature $\mathbf{s}_{r,i}^{T}\mathbf{y}_{r,i}>0$ and hence the positive definiteness of $\mathbf{B}_{r,i}$, for any $r\!\geq\!0$ and $i\!\in\![K]$. As a result, $\mathbf{B}_{r}\!\succ\!0$ on $\mathbb{R}^{d}$. This suggests the usage of the conjugate gradient (CG) method to compute the search direction at time $(s,t)$, namely $\mathbf{p}_{s,t}\triangleq-\mathbf{B}_{r}^{-1}\mathbf{v}_{s,t}$, via solving the positive definite system $\mathbf{B}_{r}\mathbf{p}_{s,t}=-\mathbf{v}_{s,t}$. In particular, for any $\mathbf{z}\in\mathbb{R}^{d}$, $\mathbf{B}_{r}\mathbf{z}$ and $\mathbf{z}^{T}\mathbf{B}_{r}\mathbf{z}$ in CG can be computed very efficiently using (59) and (61). For example,

[TABLE]

where $\mathbf{z}_{i}\triangleq\mathbf{U}_{i}\mathbf{z}$. We observe that the total computational cost in (62) is $O(M^{\prime}({M^{\prime}}^{2}+d^{\prime}))$, where $d^{\prime}\triangleq\sum_{i=1}^{K}d_{i}$. For sparse data, usually $d^{\prime}=O(d)$, so this cost is still linear in $d$. In addition, we can compute (62) in parallel across $i\in[K]$. (Intuitively, this amounts to collecting curvature from each function $\phi_{i}$ in parallel.) In this case, the computational time will be greatly reduced to $O(M^{\prime}({M^{\prime}}^{2}+\max_{i}d_{i}))$. Since typically $\max_{i}d_{i}\ll d$, the computational savings from parallel curvature collection can be significant.777The memory parameter $M$ (note that $M^{\prime}\leq M$) is usually set to a small constant, e.g., 5 or 10. Thus it has less effect on the computational complexity compared to $d$ or $\max_{i}d_{i}$.

Remark 13*.*

Note that if we interpret the matrices $\{\mathbf{U}_{i}^{T}\}_{i\in[K]}$ as sketching matrices [36], then the acceleration technique in Section VI-C can be regarded as a way of performing (approximate) Hessian sketching. However, most existing methods in the literature [37, 38, 17, 39] either use random sketching matrices or the (deterministic) frequent directions approach [40], which is based on the singular value decomposition. A certain amount of information contained in the Hessian is lost or modified in these sketching processes. In contrast, by using the sparse binary matrices $\{\mathbf{U}_{i}^{T}\}_{i\in[K]}$, our approach merely (deterministically) compresses the large sparse Hessian matrix into $K$ small dense Hessians, without changing any information contained therein.

Remark 14*.*

In [17], the authors proposed another acceleration strategy called the block BFGS update [41, 42]. We remark that this strategy can be straightforwardly combined with all the other acceleration strategies proposed above, and may result in further acceleration of the convergence of Algorithm 1.

VII Numerical Experiments

VII-A Experimental Setup

We consider two ERM problems, including logistic regression (with Tikhonov regularization) and ridge regression. For logistic regression, $b_{i}\in\{-1,1\}$ and

[TABLE]

For ridge regression, $b_{i}\in\mathbb{R}$ and

[TABLE]

Accordingly, define $f^{\log}\!\!\triangleq\!\frac{1}{n}\!\sum_{i=1}^{n}f_{i}^{\log}$ and $f^{{\rm rid}}\!\!\triangleq\!\frac{1}{n}\!\sum_{i=1}^{n}f_{i}^{{\rm rid}}$. Simple calculations reveal that the smoothness parameters $L_{i}$ of $f_{i}^{\log}$ and $f_{i}^{{\rm rid}}$ are given by $\left\|\mathbf{a}_{i}\right\|^{2}\!\!\!/4+\!\lambda$ and $2\left\|\mathbf{a}_{i}\right\|^{2}\!+\!\lambda$, respectively. Define the data matrix $\mathbf{A}\triangleq[\mathbf{a}_{1},\ldots,\mathbf{a}_{n}]$. The condition numbers $\kappa$ of $f^{\log}$ and $f^{{\rm rid}}$ are given by $(\frac{1}{4n}\sigma_{\max}^{2}(\mathbf{A})+\lambda)/\lambda$ and $({2}\sigma_{\max}^{2}(\mathbf{A})\!+\!n\lambda)/({2}\sigma_{\min}^{2}(\mathbf{A})+n\lambda)$ respectively, where $\sigma_{\max}(\mathbf{A})$ and $\sigma_{\min}(\mathbf{A})$ denote the largest and smallest singular values of $\mathbf{A}$ respectively. In both (63) and (64), we choose $\lambda=1/n$, following the convention in the literature (e.g., [17]).

We tested logistic and ridge regression problems on rcv1.binary and E2006-tfidf datasets respectively [43]. (In the sequel we abbreviate them as rcv1 and E2006-tf.) From (57), we defined a sparsity estimate of $\nabla^{2}f$ at any $\mathbf{x}\in\mathbb{R}^{d}$ as $\varrho\triangleq\left\lvert\operatorname{supp}(\mathbf{A}\mathbf{A}^{T}+\mathbf{I})\right\rvert/d^{2}$. The statistics of both datasets, including the (ambient) data dimension $d$, number of data samples $n$, sparsity parameter $\varrho$ and condition number $\kappa$ (of $f^{\log}$ or $f^{{\rm rid}}$ defined by the datasets), are summarized in Table I.888The data dimension $d$ for the original E2006-tf dataset is 150360. Due to memory issues, we randomly subsampled its features so that $d=15036$. From it, we observe that both datasets are large-scale and sparse, but with different $d$-to-$n$ ratios and condition numbers. Through these differences, we are able to infer the reasons for the different performances of some acceleration strategies on different ERM problems (shown in Section VII-B).

For both datasets, the norms of all feature vectors $\{\mathbf{a}_{i}\}_{i=1}^{n}$ have been normalized to unit. Since the smoothness parameters $L_{i}$ for both ERM problems are only dependent on $\left\|\mathbf{a}_{i}\right\|$ and $\lambda$, we have $L_{i}=L_{j}$ for any $i,j\in[n]$. Therefore the nonuniform distribution $P$ in Section III becomes uniform, and the merit of nonuniform sampling of $\mathcal{B}_{s,t}$ cannot be observed.

To estimate the global optimum $\mathbf{x}^{*}$ as ground truth, we used batch L-BFGS-B algorithm [44]. We randomly initialized $\mathbf{x}^{0}$ according to a scaled standard normal distribution. (The performance of our algorithms were observed to be insensitive to the initialization of $\mathbf{x}^{0}$.) We used the number of data passes (i.e., number of data points accessed divided by $n$), rather than running time, to measure the convergence rates of all the algorithms under comparison. This has been a well-established convention in the literature on both stochastic first-order [5, 6, 7] and second-order methods [15, 16, 17] to make convergence results agnostic to the actual implementation of the algorithms, e.g., programming languages.999Specifically, our algorithm (Algorithm 1) was implemented in Matlab®. However, some benchmarking algorithms (see Section VII-C) were implemented in other languages, e.g., the SVRG algorithm [5] and the stochastic L-BFGS algorithm [16] were implemented in C++ and Julia respectively. The differences in programming language typically have a significant impact on the actual running time of algorithms.

Finally we describe the parameter setting. We set the mini-batch size $b=\sqrt{n}$, Hessian update period $\Upsilon=10$ and the memory parameter $M=10$. We set $b_{\mathrm{H}}=b\Upsilon$ so that the computation for $\mathbf{y}_{r}$ can be amortized to each inner iteration. We set the number of inner iterations $m=n/b$, so that each outer iteration will access $2n$ data points. Lastly, we set $\eta=1\times 10^{-2}$. From Figure 1, we observe that when $\eta$ is too large, e.g., $\eta=0.1$, Algorithm 1 only converges sublinearly; whereas when $\eta$ is too small, e.g., $\eta=1\times 10^{-3}$, Algorithm 1 converges linearly but slowly. This corresponds well to our theoretical analysis in Theorem 1, which indicates that when $\eta$ falls below a threshold, $\rho$ increases as $\eta$ decreases. For both ERM problems, we see that $\eta=1\times 10^{-2}$ yields fast linear convergence.

VII-B Performance of Acceleration Strategies

We first examine the performance of Algorithm 1 with different schemes of choosing $\mathbf{x}^{s+1}$. We consider five schemes in total, including (a) uniform sampling (option I), (b) uniform averaging (option II), (c) geometric sampling (option III), (d) geometric averaging (option IV) and (e) last inner iterate (in Remark 11). For options III and IV, we set $\beta=1/2$. From Figure 2, we observe that options III and IV perform as well as the “last inner iterate” scheme, on both ERM problems, and outperform the schemes based on uniform sampling/averaging significantly. For all the subsequent experiments, we use option IV to select $\mathbf{x}^{s+1}$.

We next compare the performance of Algorithm 1 with and without using the subsampled gradient stabilization strategy in Section VI-B. As suggested by Remark 12, we chose $\widetilde{b}_{s}=\min\{\zeta\upsilon^{s},n\}$, where $\zeta=n/\upsilon^{q}$, $\upsilon=3$ and $q=8$. That is, the number of component gradients in $\widetilde{\mathbf{g}}_{s}$ exponentially increases in the first $p=8$ outer iterations and then remains at $n$. From Figure 3, we observe that this simple parameter selection method works well on both ERM problems, especially in the initial phase (when $s$ is small). In addition, we also observe when $s$ is large, both algorithms have almost the same (linear) convergence rates. This corroborates our analysis in Proposition 3.

We finally compare the performance of Algorithm 1 with and without using the low-dimensional approximate Hessian strategy in Section VI-C. We set the number of partitions $K=5$ and partition $[n]$ evenly and randomly. From Figure 4, we observe that our strategy leads to improvements of convergence on both logistic and ridge regression problems, and the improvement on the latter is very significant. It could be possible that the smaller condition number of the E2006-tf dataset enables more effective curvature collection by low-dimensional Hessians. Nevertheless, for the rcv1 dataset, which has a large condition number, our strategy is still efficacious. Additionally, we observe that our strategy preserves the linear convergence of Algorithm 1 on both problems. (A theoretical analysis of this linear convergence is left to future work.) Figure 4 shows the performance of our strategy in a single-threaded mode; as discussed in Section VI-C, our strategy can be much more efficient under scenarios where parallel computational resources are available.

VII-C Comparison to Other Algorithms

We combined all of our acceleration strategies in Section VI and compared the resulting algorithm with three benchmarking algorithms, including SVRG in [5] (with mini-batch sampling of $\mathcal{B}_{s,t}$) and two state-of-the-art stochastic L-BFGS algorithms in [16] and [17]. For the algorithm in [17], we focused on its variant (b), since it consistently outperformed other variants in experiments. We tested all the three benchmarking algorithms on both ERM problems with different step sizes $\eta\in\{10^{-4},10^{-3},\ldots,1\}$ and selected the best $\eta$ for each algorithm. The outer iterate $\mathbf{x}^{s+1}$ in all these algorithms were selected via “the last inner iterate” scheme. From Figure 5, we observe that on both ERM problems, our algorithm yields faster convergence compared to all the benchmarking algorithms. This is due to the incorporation of the acceleration strategies in Sections VI-B and VI-C. The improvement of convergence is particularly significant on the ridge regression problem. This observation is consistent with our observations in Section VII-B. In addition, we indeed observe that with the aid of curvature information, all the stochastic L-BFGS methods outperform the stochastic first-order method SVRG.

VIII Future work and an open problem

We propose to pursue future work in the following two directions. First, we aim to develop and analyze proximal and momentum-based accelerated stochastic L-BFGS algorithms, based on our coordinate transformation framework. The proximal variant enables our algorithm to be applied to composite nonsmooth (convex) objective function. The accelerated variant can potentially improve the linear convergence rate of our algorithm, and thus reduce the total computational complexity. Second, we aim to analyze the convergence of the strategy in Section VI-C. In particular, Figure 4 suggests that Algorithm 1 may still converge linearly under this strategy.

Besides future work, there is also an open problem we hope to resolve. Although we have improved the linear convergence rate and computational complexity of Algorithm 1 as compared to those in [16] and [17], it seems our improved complexity in (49) is still inferior to that of SVRG. In SVRG, the complexity is $O((n+\kappa)d\log(1/\epsilon))$ [27], so our complexity (49) has an additional multiplicative factor $\kappa_{\mathrm{H}}$. This contradicts the experimental results in [16], [17] and Section VII-C, where stochastic L-BFGS-type algorithms have been repeatedly shown to outperform their first-order counterparts. *Therefore, an interesting problem consists in obtaining a (computational) complexity bound of the stochastic L-BFGS algorithm that is better (or at least as good as) that of SVRG. * Indeed, a careful analysis reveals that the additional $\kappa_{\mathrm{H}}$ arises from the uniform spectral bound of the metric matrices $\{\mathbf{H}_{r}\}_{r\geq 0}$. This uniform bound is effectively a worst-case bound, and does not reflect the local curvature information contained in recent iterates at any time $(s,t)$. Since the judicious use of curvature information serves as a very important reason for the fast convergence of the stochastic quasi-Newton algorithms, such information should also be reflected in theoretical analysis as well (possibly in an adaptive spectral bound for $\{\mathbf{H}_{r}\}_{r\geq 0}$). We believe an effective adaptive bound is critical for improving the complexity result in (49).

Interestingly, an incremental quasi-Newton (IQN) method was proposed by Mokhtari et al. [14] recently. The proposed algorithm makes use of the aggregated optimization variables, as well as the aggregated gradients and approximate Hessians of all the component functions to reduce the noise of gradient and Hessian approximations. As a result, it achieves the local superlinear convergence rate, but requires $\Theta(nd^{2})$ storage space. (Note that most of the stochastic L-BFGS methods, such as [15, 16] and Algorithm 1, only require $\Theta(d)$ memory.) The key idea in [14] is to show the descent direction in the IQN algorithm asympotically converges to that of the Newton’s method. If this condition holds, then the additional factor $\kappa_{\mathrm{H}}$ in the complexity estimate (49) may be removed. Therefore, how to design a stochastic quasi-Newton algorithm that satisfies this condition in the momeory-limited setting would also be an interesting problem to pursue in the future.

Appendix A Proof of Lemma 3

The proof of Lemma 3 is shown in (65) to (71) on page A.

In (67), we use the independence of $i_{j}$ and $i_{j^{\prime}}$ for $j\neq j^{\prime}$ and (95) in Lemma T-2. In (68), we use (95) and the fact that \mathbb{E}\big{[}\left\|\mathbf{a}-\mathbb{E}[\mathbf{a}]\right\|^{2}\big{|}\mathcal{G}\big{]}\leq\mathbb{E}\big{[}\left\|\mathbf{a}\right\|^{2}\big{|}\mathcal{G}\big{]} almost surely, for any $\sigma$-algebra $\mathcal{G}$. In addition, the inequalities in (70) and (71) follow from $\left\|\mathbf{a}+\mathbf{b}\right\|^{2}\leq 2\left\|\mathbf{a}\right\|^{2}+2\left\|\mathbf{b}\right\|^{2}$ and (96) respectively.

Appendix B Proof of Lemma 4

Our proof is inspired by [17]. For any $\mathbf{x}\in\mathbb{R}^{d}$ and $\mathcal{T}\subseteq[n]$ with cardinality $b_{\mathrm{H}}$, we have

[TABLE]

Define $\mathbf{V}_{k}\triangleq\mathbf{I}-{\mathbf{y}_{k}{\mathbf{s}_{k}}^{T}}/({\mathbf{y}_{k}}^{T}\mathbf{s}_{k})$ and $\mathbf{Q}_{k}\triangleq{\mathbf{s}_{k}{\mathbf{s}_{k}}^{T}}/({\mathbf{y}_{k}}^{T}\mathbf{s}_{k})$, then (6) becomes

[TABLE]

Fix any $r\in\mathbb{N}$. Since $\mathbf{y}_{k}=\nabla^{2}f_{\mathcal{T}_{r}}(\overline{\bf x}_{r})\mathbf{s}_{k}$, we have

[TABLE]

where $\widetilde{\mathbf{s}}_{k}\triangleq\nabla^{2}f_{\mathcal{T}_{r}}(\overline{\bf x}_{r})^{1/2}\mathbf{s}_{k}$. Hence

[TABLE]

Similarly, $\left\|\mathbf{Q}_{k}\right\|={\left\|\mathbf{s}_{k}\right\|^{2}}/({\mathbf{s}_{k}}^{T}\nabla^{2}f_{\mathcal{T}_{r}}(\overline{\bf x}_{r})\mathbf{s}_{k})\leq{1}/\overline{\mu}_{b_{\mathrm{H}}}$. By (73),

[TABLE]

We apply (74) repeatedly over $k\!=\!r-M^{\prime}+1,\ldots,r$, then

[TABLE]

where (a) follows from the definition of $\mathbf{H}_{r}^{(r-M^{\prime})}$ in (5) and

[TABLE]

and (b) follows from $\kappa_{b_{\mathrm{H}}}\geq 1$ and $M^{\prime}\leq M$. To show $\gamma={1}/{(M+1)\overline{L}_{b_{\mathrm{H}}}}$, it suffices to show $\left\|\mathbf{B}_{r}\right\|\leq(M+1)\overline{L}_{b_{\mathrm{H}}}$. We derive this bound using (18). Since

[TABLE]

and

[TABLE]

we have from (18) that

[TABLE]

Therefore,

[TABLE]

where (a) follows from (76).

Appendix C Proof of Proposition 2

Our proof is modified from that of Theorem 1. Specifically, our proof leverages a refined telescoping of (34). (The steps up to (34) are unchanged.) For each $t\in[m]$, we multiply both sides of (34) by $(1-\eta\gamma\overline{\mu})^{m-t}$ and obtain

[TABLE]

Telescope (77) over $t=1,\ldots,m$ and we have

[TABLE]

Now we consider using option IV to choose $\mathbf{x}^{s+1}$. Since $0<\beta\leq 1-\eta\gamma\overline{\mu}$, using (11) and Jensen’s inequality, we have

[TABLE]

Alternatively, if $\mathbf{x}^{s+1}$ is determined using option III, the definition of distribution $Q$ still yields (79). Finally, using (93) in Lemma T-1 to bound $\|\widetilde{\mathbf{x}}^{s,r^{\prime}}-\widetilde{\mathbf{x}}_{r^{\prime}}^{*}\|^{2}$ in (78), we have

[TABLE]

Taking expectation on both sides and using (11), we arrive at (52).

Appendix D Proof of Proposition 3

The subsampled gradient strategy essentially introduces errors in $\{\mathbf{g}_{s}\}_{s\geq 0}$. Therefore, we explicit model this error by $\mathbf{e}_{s}\triangleq\widetilde{\mathbf{g}}_{s}-\mathbf{g}_{s}$. We first bound the second moment of $\mathbf{e}_{s}$. Since $\mathbf{g}_{s}=\nabla f(\mathbf{x}^{s})$, by Lemma T-3 and $\widetilde{b}_{s}\geq nS^{2}\alpha_{s}/(S^{2}\alpha_{s}+(n-1)\xi^{2}\overline{\rho}^{2s})$, we have for any $s\in(S]$,

[TABLE]

As a result,

[TABLE]

The introduction of random sets $\{\widetilde{\mathcal{B}}_{s}\}_{s\geq 0}$ requires us to redefine the filtration $\{\mathcal{F}_{s,t}\}_{s\geq 0,t\in(m-1]}$ as

[TABLE]

As usual, we define $\mathcal{F}_{s}\triangleq\mathcal{F}_{s,0}$. In the sequel, we follow the naming convention in coordinate transformation framework as in Section IV-A. In particular, we define $\widetilde{\mathbf{e}}_{s,r^{\prime}}\triangleq\mathbf{H}_{r^{\prime}}^{1/2}\mathbf{e}_{s}$. Now, define $\widetilde{\mathbf{v}}^{\prime}_{s,t,r}\triangleq\widetilde{\mathbf{v}}_{s,t,r}+\widetilde{\mathbf{e}}_{s,r^{\prime}}$, then from (24) and (25), we have

[TABLE]

Define $\widetilde{\bm{\delta}}^{\prime}_{s,t,r}\triangleq\widetilde{\mathbf{v}}^{\prime}_{s,t,r}-\nabla\widetilde{f}_{r}(\widetilde{\mathbf{x}}_{s,t,r})$. We can derive an inequality similar to (33) in the proof of Theorem 1, i.e.,

[TABLE]

Taking expectation w.r.t. $\mathcal{B}_{s,t}$ and using (84) and (85), we have

[TABLE]

Using Cauchy-Schwartz inequality, we have

[TABLE]

Substituting (87) into (86), and using the telescoping techniques in Appendix C, we have

[TABLE]

With either option III or IV and (93) in Lemma T-1, we have

[TABLE]

Taking expectation on both sides and using (11), we have

[TABLE]

From (81) and (82), we have

[TABLE]

Substituting (90) and (91) into (89), we have

[TABLE]

Applying (92) recursively and we reach (54).

Appendix E Technical Lemmas

Lemmas T-1. T-2, and T-3 can be found in [45, Chapter 9], [27], and [46] respectively.

Lemma T-1**.**

If a function $f:\mathbb{R}^{d}\to\mathbb{R}$ is $\mu$-strongly convex and $L$-smooth, then for any $\mathbf{x}\in\mathbb{R}^{d}$,

[TABLE]

where $\mathbf{x}^{*}$ denotes the unique minimizer of $f$ on $\mathbb{R}^{d}$.

Lemma T-2**.**

Let $f$ be defined as in $\eqref{eq:problem}$ and satisfies Assumption 2. Define a distribution $p$ with support $[n]$ such that $p_{i}=L_{i}/(n\overline{L})$, for any $i\in[n]$. Then for any $\mathbf{x}\in\mathbb{R}^{d}$,

[TABLE]

where $\mathbf{x}^{*}$ denotes the unique minimizer of $f$ on $\mathbb{R}^{d}$.

Lemma T-3**.**

Let $\{\mathbf{z}_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{d}$ and define $\overline{\bf z}\triangleq 1/n\sum_{i=1}^{n}\mathbf{z}_{i}$. Uniformly sample a random subset $\mathcal{S}$ of $[n]$ with size $b$ without replacement. Then

[TABLE]